CHAVEZ_CAROLINE_hw1

IMPORTANT: Name your file using your name: “LAST_FIRST_hw1”

Whenever you write your own function, it has to be annotated, describing what the function does and what input arguments it takes.

For the homework to be considered complete, submit it as both an .rmd file, as well as an .html or .pfd version using the ‘Knit’ option.

  1. create an array containing numbers from 1 to 100 and name it n100. TIP: use the ":" annotation to create the 100 numbers. (1p)
n100=c(1:100)
n100
##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##  [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
##  [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
##  [55]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
##  [73]  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
##  [91]  91  92  93  94  95  96  97  98  99 100
  1. Create a new variable called v50 by sampling 50 numbers at random from your ‘n100’ array. Sample with replacement. Then create v50b by sampling 50 numbers from n100 without replacement. (3p)
v50=sample(n100,50,replace=T)
v50
##  [1] 67 97 12 81 70 17 88 37 51  9 97 26 76 70 34 39 97 90 71 26 86 65 71 81 59
## [26] 61 89 68 68 50 48 95 62 68 79 45 10 41 45 58  9 80 52 41 74 61 74 21 36  6
v50b=sample(n100,50,replace=F)
v50b
##  [1]   4  35  47 100  16  95  29  62  92  79  55  27   2  57  26  89  85  78  82
## [20]  52  74  51  42  18  25  34  22  64  23  13  44  37  63  98  32  67  66  39
## [39]  90  73  61  58  43   3   6  91  49  14  10  31
  1. How many numbers occur more than once in v50? How many in v50b? Tip: you can do this by using the necessary function twice. (2p)
sum(table(v50) > 1)
## [1] 11
sum(table(v50b) > 1)
## [1] 0
  1. Create a function named my_mean that calculates the mean, without using the mean() base function. You are allowed to use other base functions and if you forgot how to compute the mean manually, feel free to look it up. Compute the mean of v1 using your my_mean function and compare it to R’s base function mean(). Are they the same? (3p)
my_mean=function(n){(sum(n))/length(n)}
my_mean(v50)
## [1] 57.16
mean(v50)
## [1] 57.16
  1. What is the standard deviation of v50? Compute it manually by relying on your my_mean function. Name the function my_stdev. Note, my_stdev should take a vector as input, in this case your v50 vector, and it should call my_mean from within the my_stdev function. Do you get the same result as R’s sd() base function ? (3p)
my_stdev=function(n){sqrt(sum((n-my_mean(n))^2)/(length(n)-1))}
my_stdev(v50)
## [1] 26.15758
sd(v50)
## [1] 26.15758
  1. Make a data.frame called df10 with 10 rows and 2 columns named col1 and col2 by sub-sampling from v50 without replacement twice (once for col1 and once for col2). Display the first 5 rows of this data.frame df10? (3p)
df10=data.frame(col1=c(sample(v50,10,replace=F)), col2=c(sample(v50,10,replace=F)))
df10[1:5,1:2]
##   col1 col2
## 1   61   26
## 2    9   67
## 3   74   21
## 4   12   80
## 5   68   68
  1. Sort the numbers in v50 in a descending manner and print the last 4 elements of the sorted vector. (2p)
sort(v50, decreasing=T)
##  [1] 97 97 97 95 90 89 88 86 81 81 80 79 76 74 74 71 71 70 70 68 68 68 67 65 62
## [26] 61 61 59 58 52 51 50 48 45 45 41 41 39 37 36 34 26 26 21 17 12 10  9  9  6
tail(sort(v50,decreasing=T), n=4)
## [1] 10  9  9  6
  1. How else could you have found the last number of this sequence? Hint: Use a function from the lecture slides. (1p)
rev(sort(v50,decreasing=T))
##  [1]  6  9  9 10 12 17 21 26 26 34 36 37 39 41 41 45 45 48 50 51 52 58 59 61 61
## [26] 62 65 67 68 68 68 70 70 71 71 74 74 76 79 80 81 81 86 88 89 90 95 97 97 97
rev(sort(v50,decreasing=T))[1]
## [1] 6
  1. What class is v50? (1p)
class(v50)
## [1] "integer"
  1. Turn v50 from its current class into a numeric vector. Store this new vector as v50num and verify that the class has changed. (2p)
v50num=as.numeric(v50)
class(v50num)
## [1] "numeric"
  1. Turn it into a factor vector called v50fac. Do you have more, less, or equal levels as there are entries in v50fac? EXPLAIN the result! (2p)
v50fac=factor(v50)
class(v50fac)
## [1] "factor"
v50fac
##  [1] 67 97 12 81 70 17 88 37 51 9  97 26 76 70 34 39 97 90 71 26 86 65 71 81 59
## [26] 61 89 68 68 50 48 95 62 68 79 45 10 41 45 58 9  80 52 41 74 61 74 21 36 6 
## 37 Levels: 6 9 10 12 17 21 26 34 36 37 39 41 45 48 50 51 52 58 59 61 62 ... 97
#There are less entries in v50fac than v50num because v50fac treats each number as its own level, so if there are repeat numbers (which there are because it is sampled with replacement), it would just be counted as one level.
  1. Create my_fun_function that takes v50 and another vector of your own design as input. Write a function of your choosing, following these instructions: (3p)
my_fun_function=function(v50,rnum){abs(max(v50)-min(rnum))}
my_fun_function(v50,sample(1:10,5,replace=F))
## [1] 95
  1. Bonus: You work in a protein engineering lab and want to optimize your protein’s binding pocket. You decide to start by generating flexible amino acid linker regions that are 16 amino acids long and contain only aliphatic amino acids (glycine - G, alanine - A, valine - V, leucine - L , and isoleucine - I) and the positively charged amino acids Arginine-R and Lysine-K. You know from literature that the optimal R and K frequencies should be around 25% each (~4 amino acids each). Use sample() and paste() to generate 16-amino acid linkers that fulfill the above criteria. Use ?sample to learn how you can specify sampling probabilites. Create 5 different versions and ensure your versions are reproducible. Do all versions have exactly 25% Rs and Ks? Why do you think that is? (3p)
aamino=c('G','A','V','L','I','R','K')
set.seed(1)

p1=sample(aamino,16,T,prob=c(.10,.10,.10,.10,.10,.25,.25))
p2=sample(aamino,16,T,prob=c(.10,.10,.10,.10,.10,.25,.25))
p3=sample(aamino,16,T,prob=c(.10,.10,.10,.10,.10,.25,.25))
p4=sample(aamino,16,T,prob=c(.10,.10,.10,.10,.10,.25,.25))
p5=sample(aamino,16,T,prob=c(.10,.10,.10,.10,.10,.25,.25))

paste(p1,collapse="")
## [1] "KKVGRAGIIRRRIKLK"
paste(p2,collapse="")
## [1] "LGKLGRIRKKRKAKKV"
paste(p3,collapse="")
## [1] "KRAILRLKAILVVLRK"
paste(p4,collapse="")
## [1] "LIKAKRRRKVIKGKKK"
paste(p5,collapse="")
## [1] "IKKLRAKAKKKAAKLG"
table(p1)/16*100
## p1
##     A     G     I     K     L     R     V 
##  6.25 12.50 18.75 25.00  6.25 25.00  6.25
table(p2)/16*100
## p2
##     A     G     I     K     L     R     V 
##  6.25 12.50  6.25 37.50 12.50 18.75  6.25
table(p3)/16*100
## p3
##     A     I     K     L     R     V 
## 12.50 12.50 18.75 25.00 18.75 12.50
table(p4)/16*100
## p4
##     A     G     I     K     L     R     V 
##  6.25  6.25 12.50 43.75  6.25 18.75  6.25
table(p5)/16*100
## p5
##     A     G     I     K     L     R 
## 25.00  6.25  6.25 43.75 12.50  6.25
#Not all of the versions have exactly 25% R and K because when we set the probability to 0.25, it just means 25% for each individual protein to be chosen during the sampl, and does not guarantee that 4 of 16 will be R and 4 will be K. This is because of random sampling which does not guarantee an exact split in results.