r simulation代写 | Programming Projects | DNA/protein alignment biology

首页/ r simulation代写 | Programming Projects | DNA/protein alignment biology

Programming Projects

You may use any language: python, perl, matlab, R, C, C++, Java, etc.

Your solution should include

 

  • Explanation of your approach, implementation, and results in well organized and grammatical text.
  • source code
  • instructions to install and run the code
  • test data and expected results
  • Answers to any specific questions

 

You may use existing packages/libraries to read and write sequences, but, unless otherwise instructed, generate your own code. Do not use solutions you look up on the internet.

 

Alignments

 

  1. DNA/protein alignment

 

Align a DNA sequence and a protein sequence based on the translated sequences of the DNA sequences.

Only the forward reading frames need be considered. If desired, you may use existing packages for

reading/writing sequences, scoring tables, and formatted results. Display the results in conventional pairs of

lines format showing both DNA and translated protein sequences.

  • Use a local alignment algorithm
  • Implement both global and local alignments
  • Assume the sequences are in fasta format
  • Use a scoring table such as Blosum or PAM
  • Allow affine gaps within reading frames
  • Allow length independent gaps between reading frames.
  • allow affine gaps between reading frames
  • Align two DNA sequences, same requirements as above, based on their protein translations

 

  1. Suboptimal alignment

 

For DNA or protein sequences, report n non-intersecting suboptimal alignments. As discussed in class, sub-optimal alignments use no pair of letters used in earlier alignments. This can be achieved by recalculating the score matrix, setting the score for all letter-pairs in previous alignments to be zero.

  • allow selection of an appropriate scoring table for DNA or protein
  • use a local algorithm with affine gap penalties
  • display aligned sequences in conventional pairs of lines format

 

Simulations

  1. Sequence randomization

 

Randomize a sequence preserving its n-mer word content (i.e, 1 letter, 2-letter, 3-l;etter, …). This may be solved simply using sampling with or without replacement, but the n-mer content of the randomized sequence will not exactly match the original sequence.

  • input and output in Fasta format. Your program should work for any alphanumeric letter sequence, but especially DNA and protein alphabets.
  • n-mer should be selectable in range 1 to 5 (or more)
  • Additional output should show original and randomized n-mer word content
  • Exact solution (all n-mer counts exactly the same) using an implementation of Euler’s algorithm.

 

  1. Sequence evolution – Is a simple i.i.d. randomization a good model for proteins? Choose a random set of 50 protein sequences (from ncbi or uniprot). This problem can be combined with Simulation #1
  • For the test data set, calculate the distribution of word scores using words of various length (minimum: 8, 12, 16) using a Blosum scoring table. Display as a frequency table (score vs count).
  • Display as a histogram plot
  • Compare to the scores distributions generated for the random sequences with different n-mer word lengths (n=1..5) preserved. Use the code in problem 1, above, or the program Ushuffle.
  • Mean and standard deviation are sufficient for the comparison, although you may include more sophisticated tests.
  • repeat for DNA using word lengths 4, 6, and 8.
  • What n-mer model best fits the unrelated protein data?

 

Evolution/Trees

 

  1. Starting with a protein sequence (at least 100 residues), apply a random mutation process based on the PAM model. At random times, duplicate the sequences creating speciation events. At random times, remove some species (extinctions).
  • Input parameters should include total evolution time (% true mutations), and number of leaf sequences to be generated.
  • speciation and extinctions can follow simple a simple probability, i.e., probability of speciation or extinction at each step is constant.
  • Implement a more sophisticated model where these events follow a distribution such as normal or gamma
  • report the final set of sequences as a multiple alignment. No alignment algorithm is necessary, since the sequences are generated by random mutation, they are already aligned.
  • report the tree in Newick format.

 

  • report the true number of mutations between each pair of sequences (from the simulation history) and the observed number of differences (n x n table, where n is number of leaf nodes)

 

  1. parsimony – requires solving Evolution #1, above.
  • Implement the parsimony algorithm.
  • for 50 random trees with 20 leaf nodes, calculate the tree lengths and report the true tree lengths (known from the simulation), estimated tree length (from parsimony). Calculate the mean and standard deviation of the difference between the known and estimated tree lengths.
  • Repeat for evolutionary distances of 50%, 100%, 200%. Are there significant differences between the known and estimated trees?

 

  1. Same as above for UPGMA trees. – requires solving Evolution #1, above.
  • use 50 random trees with 20 leaf nodes
  • Implement UPGMA
  • Use the Fitch-Margoliash method to estimate branch lengths
  • Use the observed mutational distances to construct trees, compare the estimated distances from the tree to the known true distances from the simulation.
  • Calculate the mean and standard deviation of the RMS difference between the known and estimated numbers of mutations.
  • Repeat for evolutionary distances of 50%, 100%, 200%. Are there significant differences between the known and estimated trees?

 

成果展示,部分。如需需要成品，左边二维码联系我们的客服。

 

r simulation代写
Simulation

1

The function seq.Rand shuffles a string sequence. The argument seq for a sequence, k for n-mer, and seed for setting set to make the result repeatable. The output will include the original one, the randamized one, and the Eulerian path.

fasta.form() function combines seq.Rand and the function for dealing with fasta format.

Codes:

# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("Biostrings", version = "3.8")
# BiocManager::install("msa")
library(Biostrings)

library(msa)
library(knitr)
data("BLOSUM62")

seed <- 2018 # an example
k <- 5 # an example
seq <- "ACTGTGCAGTCGATGCTAGGTTCCGATATCTAGCTCGCTAGATC" # an example

seq.Rand <- function(seq, k, seed = NULL){
  if(!is.null(seed)) set.seed(seed) # set seed, if any
  if(k==1){ # if k = 1, sampling the vector without replacement
    seq.v <- unlist(strsplit(seq,"")) # break a sequence to a vector, each entry only has one char
    seq.v[-1] <- seq.v[sample(2:length(seq.v))] # first entry remains the same
    seq.new <- Eulerian <- paste(seq.v,collapse = "") 
  } else{
    seq.v <- unlist(strsplit(seq,"")) # break a sequence to a vector, each entry only has one char
    letter.dict <- unique(seq)
    
    len <- length(seq.v)
    kmer.o <- sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")})
    
    seq.new <- Eulerian <- NULL
    seq.new[1] <- Eulerian[1] <- kmer.o[1]
    kmer <- kmer.o[-1]
    for(i in 2:(length(kmer)+1)){ # if k>=2
      # replace without relacement
      temp0 <- seq.new[i-1] 
      overlap <- substr(temp0, 2, k)
      dict.temp <- which(substr(kmer, 1, k-1) == overlap)
      if (length(dict.temp) == 0) break
      ind.temp <- sample(dict.temp,1)
      seq.new[i] <- kmer[ind.temp]
      kmer <- kmer[-ind.temp]
      # print(kmer)
    }
    
    #Eulerian
    kmer <- kmer.o[-1]
    for(i in 2:(length(kmer)+1)){
      temp0 <- Eulerian[i-1]
      overlap <- substr(temp0, 2, k)
      dict.temp <- which(substr(kmer, 1, k-1) == overlap)
      if (length(dict.temp) == 0) { # to make sure every node is visited
        dict.temp <- which(substr(kmer.o, 1, k-1) == overlap)
        if (length(dict.temp) == 0) break
        ind.temp <- sample(dict.temp,1)
        Eulerian[i] <- kmer.o[ind.temp]
        # print(kmer)
      } else {
        ind.temp <- sample(dict.temp,1)
        Eulerian[i] <- kmer[ind.temp]
        kmer <- kmer[-ind.temp]
        # print(kmer)
      }
    }
    
    # format the sequence
    seq.new <- paste(seq.new[1], 
                     paste(substr(seq.new[2:length(seq.new)],k-1,k),collapse = ""),
                     sep = "")
    Eulerian <- paste(Eulerian[1], 
                      paste(substr(Eulerian[2:length(Eulerian)],k-1,k),collapse = ""),
                      sep = "")
  }
  
  # output:
  return(list(randseq = seq.new,
              Eulerian = Eulerian,
              original = seq))
  
}

seq.Rand(seq, k,2018)

## $randseq
## [1] "ACTGTTGGCCAAGGTTCCGGAATTAGGTAAGGACGTCAGGACCTC"
## 
## $Eulerian
## [1] "ACTGTTGGCCAAGGTTCCGGAATTGGCCTTAAGGCCTGGGTTTTCCCCGGAATTAATTCTGGTTGTGGTGTTGTGGCTGTGGT"
## 
## $original
## [1] "ACTGTGCAGTCGATGCTAGGTTCCGATATCTAGCTCGCTAGATC"

#fasta
fasta.form <- function(fasta, k, seed=NULL){
  seq.name <-  names(fasta) # deal with fasta format
  seq <- paste(fasta)
  result <- seq.Rand(seq, k,seed)
  return(list(name = seq.name,
              randseq = result$randseq,
              Eulerian = result$Eulerian,
              original = result$original))
}

data <- readDNAStringSet("./sequence.fasta-3.txt")

## Warning in .Call2("fasta_index", filexp_list, nrec, skip, seek.first.rec, :
## reading FASTA file ./sequence.fasta-3.txt: ignored 12146 invalid one-letter
## sequence codes

fasta.form(data[1], k,2018)

## $name
## [1] "BAG70982.1 phenylalanine ammonia-lyase [Musa balbisiana]"
## 
## $randseq
## [1] "MAKAVVVVNNGGAACCKKAADDNNWWKKAAAASSTTGGSSHHDDVVKKRRMMVVRRKKVVRRGGAATTTTSSVVAAAAVVAAAAAARRAAAAVANRSSVVRRVVSSAARRDDGGVVRRAASSSASTATKHRVMWV"
## 
## $Eulerian
## [1] "MAKAVVVVNNGGAACCKKAADDNNWWKKAAAASSTTGGSSHHDDVVKKRRMMVVRRKKVVRRGGAATTTTSSVVAAAAVVAAAAGSGSAARRSSVVRAAKAACCKNGNGNGVVVVVNVVVNVVVVVVAVNTATYGDGGVVRRDKVVRRGNWVVVVVVVVVNVVVVVVAVRASSNAMNSSRVVSAVTTSTVVVNAVSMRGRGRKCKACACNGNGGAAVNGGTTDSWWVARAAAGNKRAASTDDHSAGAVRSTVVAVYSAATSSGKRKGDSGTASKAASHRGGVRKAVVAVTTGGGARTTTVMMSVRGSASKAACAVSGTKVVVNVNAVASAASHHDKANGAVAKVRAAASGANGNGVNVVVVAVAATGNWDNVVVVAVCRRRCKVNNGGAGAVNVVVVVVAVSGAKTTVTDVNGGAGAAVTANATGSMKRVNAVVRMVATSYDGVRDNACAVGTWKKAACACVNVNVVVNVVVVVNVNVVVVVVVNVVVNAVVSGSATACCKGANGVVVNVVAVVNKRCKCKNGGAVVVNVNAVAAGVAATTADKAADVVAVSYTCAARAMGNTVRVNNGNGVVVVAVGRDVAANWVVAVAGGNKAKAADNGNGVNVNVVVVVNVNVNNGAVRNVTAAVNNGGAVNVVVNAVSGRVVNGSDNKAVNNGNGVNVVVVAVAGNAGGTDARAAVNNGGAVVVNVNVNAVVGGYKAKACKVVVVAVVVVNVVVNAVVGTGKGAVAADNVNNGGAACAVVGYAATTSRDSHKACKAVVVVNVNVNNGGAAVVGNNVMKVWKWKDNGAVNNGGAVVAVVGTGWKACCKKAACNGAVVGGRGDKANTNKKAWKVNVVAVVVAVVGGDGSRDTSVVVNNGNGVVVNNGNGAVVGTGRGRVAAAAAVVVVNNGVVVNVVVVVVVVAVVVVVVN"
## 
## $original
## [1] "MAKAVVNGACKADNWKAASTGSHDVKRMVRKVRGATTSVAAVAAARSVRVSARDGVRASSWVMSMNKGTDSYGVTTGGATSHRRTKGGAKRNAGGSGSGNTSSAAKAAMVRVNTGYSGRAASNNGTCRGTTASGDVSYAGTGRNAKAVGDGKVGAAARASADGKGAVNGTAVGSGASMVANVAVSVSAVCVMGKTDHTHKKHHGAAAMHGSSYMKMAKKHDKKDRYARTSWGVRASTKSRNSVNDNDVSRSKAHGGNGTGVSMDNTRAVAAGKMASVNDYNNGSNSGGRNSDYGKGAAMAAYCSANVTNHVSAHNDVNSGSARKTAAVDKMSTTYVACAVDRHNKSAVKSTVSVAKRVTMGANGHARCKKVVDRHVTYVDDCSATYMKRVVAHANGDKKDAGSSKATTAKVAARAAVGGKAANRCRSYYRVRKTGTGKVRSGDKVDACGKVDCKWNGAC"

2

the distribution of word scores

Please refer to the following code and out put for the frequency tables and histgrams for the scores derived by different lengths of word.

The mean and Variance are summarized as followings.
length of word
 
Mean
 
Variance
 
8
 
32.31412
 
41.9872
 
12
 
47.30298
 
104.3903
 
16
 
62.61114
 
190.2871
 
The larger the length of word, the higher mean and variance since the maginitude of the scores increaase.

shuffle

Please refer to the following code and out put for the shuffled word, which derives the the frequency tables and histgrams for shuffled words with different n-mer.

The mean and Variance of scores for different n-mer are summarized as followings.
n
 
Mean
 
Variance
 
1
 
25.559
 
7.6516
 
2
 
25.761
 
9.493
 
3
 
25.66886
 
8.859
 
4
 
24.917
 
4.986
 
5
 
26.48
 
8.758
 
What n-mer model best fits the unrelated protein data?

The 5-mer obtains larger mean, indicating 5-mer model is the best one among the above five mdoels to fit the unrelated protein data.

Redo for DNA

Please refer to the code and out put for DNA for details.

The mean and Variance are summarized as followings.
length of word
 
Mean
 
Variance
 
4
 
32.31412
 
41.9872
 
6
 
47.30298
 
104.3903
 
8
 
62.61114
 
190.2871
 
code and out of the score.

lengthof word: 8

data <- readDNAStringSet("./sequence.fasta-3.txt")

## Warning in .Call2("fasta_index", filexp_list, nrec, skip, seek.first.rec, :
## reading FASTA file ./sequence.fasta-3.txt: ignored 12146 invalid one-letter
## sequence codes

seq.name <-  names(data)
seq <- paste(data)
data <- data.frame(seq.name, seq)
data$seq <- as.character(data$seq)
data$seq.name <- as.character(data$seq.name)

k <- 8
word.list <- NULL
for(i in 1:50){
  seq.v <- unlist(strsplit(seq[i],""))
  len <- length(seq.v)
  word.list <- c(word.list, sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")}))
}

word.list <- unique(word.list)

# seq.v <- unlist(strsplit(seq,""))
# letter.dict <- unique(seq)
# 
# len <- length(seq.v)
# kmer.o <- sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")})

scores.list <- list()
i <- 1
  seq.v <- unlist(strsplit(seq[i],""))
  len <- length(seq.v)
  scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
  for(j in 1:length(word.list)){
    scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
  }
  scores.temp <- as.vector(scores.temp)
  scores.list[[i]] <- scores.temp

socres.8 <- unlist(scores.list)
T.8 <- quantile(socres.8,.999) #set a T
socres.8 <- socres.8[socres.8>=T.8]
table(socres.8)

## socres.8
##  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41 
## 238 223 191 152 163 129 120 131 114  95 111 103 115 107 105  89  82  88 
##  42  43  44  45  46  47  48  49  50  51  52  53  55  56 
##  63  57  45  34  20  16  10   4   5   4   2   2   1   1

hist(socres.8)




mean(socres.8)

## [1] 32.31412

var(socres.8)

## [1] 41.9872

lengthof word: 12

# 12
k <- 12
word.list <- NULL
for(i in 1:50){
  seq.v <- unlist(strsplit(seq[i],""))
  len <- length(seq.v)
  word.list <- c(word.list, sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")}))
}

word.list <- unique(word.list)

# seq.v <- unlist(strsplit(seq,""))
# letter.dict <- unique(seq)
# 
# len <- length(seq.v)
# kmer.o <- sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")})

scores.list <- list()
i <- 1
seq.v <- unlist(strsplit(seq[i],""))
len <- length(seq.v)
scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
for(j in 1:length(word.list)){
  scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
}
scores.temp <- as.vector(scores.temp)
scores.list[[i]] <- scores.temp

socres.12 <- unlist(scores.list)
T.12 <- quantile(socres.12,.999)
socres.12<- socres.12[socres.12>=T.12]
table(socres.12)

## socres.12
##  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49 
## 135 125 118 115 105 108 114  90 102 109  93  88 114  90  97  72  92  81 
##  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67 
##  97  97  73  94  87  88  90  74  86  72  81  70  75  51  50  46  25  32 
##  68  69  70  71  72  73  74  75  77  80 
##  15  11   9   4   2   1   2   3   1   1

hist(socres.12)




mean(socres.12)

## [1] 47.30298

var(socres.12)

## [1] 104.3903

lengthof word: 16

#16
k <- 16
word.list <- NULL
for(i in 1:50){
  seq.v <- unlist(strsplit(seq[i],""))
  len <- length(seq.v)
  word.list <- c(word.list, sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")}))
}

word.list <- unique(word.list)

# seq.v <- unlist(strsplit(seq,""))
# letter.dict <- unique(seq)
# 
# len <- length(seq.v)
# kmer.o <- sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")})

scores.list <- list()
i <- 1
seq.v <- unlist(strsplit(seq[i],""))
len <- length(seq.v)
scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
for(j in 1:length(word.list)){
  scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
}
scores.temp <- as.vector(scores.temp)
scores.list[[i]] <- scores.temp

socres.16 <- unlist(scores.list)
T.16 <- quantile(socres.16,.999)
socres.16 <- socres.16[socres.16>=T.16]
table(socres.16)

## socres.16
##  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 
## 131 113  91 100  99  85  93 100  90  99  73  75  74  75  77  86  72  72 
##  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77 
##  62  72  64  60  66  70  56  67  66  60  71  75  72  75  61  92  71  72 
##  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  94  95  96 
##  67  88  71  63  57  52  51  49  44  28  25  13   5  15   6   2   3   2 
##  97  98  99 
##   2   1   1

hist(socres.16)




mean(socres.16)

## [1] 62.61114

var(socres.16)

## [1] 190.2871

code and out put for shuffle

## shuffle
shuffle.1 <- seq.Rand(seq[1],1,2018)$randseq
shuffle.2 <- seq.Rand(seq[1],2,2018)$randseq
shuffle.3 <- seq.Rand(seq[1],3,2018)$randseq
shuffle.4 <- seq.Rand(seq[1],4,2018)$randseq
shuffle.5 <- seq.Rand(seq[1],5,2018)$randseq

# 1mer
seq.v <- unlist(strsplit(shuffle.1,""))
len <- length(seq.v)
scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
for(j in 1:length(word.list)){
  scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
}
scores.temp <- as.vector(scores.temp)
scores.list[[i]] <- scores.temp

socres.16 <- unlist(scores.list)
T.16 <- quantile(socres.16,.999)
socres.16 <- socres.16[socres.16>=T.16]
hist(socres.16)




mean(socres.16)

## [1] 25.559

var(socres.16)

## [1] 7.651598

kable(table(socres.16))
socres.16
 
Freq
 
23
 
965
 
24
 
724
 
25
 
564
 
26
 
420
 
27
 
322
 
28
 
237
 
29
 
149
 
30
 
103
 
31
 
61
 
32
 
53
 
33
 
49
 
34
 
26
 
35
 
11
 
36
 
8
 
37
 
8
 
38
 
2
 
39
 
2
 
40
 
1
 
41
 
4
 
42
 
1
 
43
 
2
 
# 2mer
seq.v <- unlist(strsplit(shuffle.2,""))
len <- length(seq.v)
scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
for(j in 1:length(word.list)){
  scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
}
scores.temp <- as.vector(scores.temp)
scores.list[[i]] <- scores.temp

socres.16 <- unlist(scores.list)
T.16 <- quantile(socres.16,.999)
socres.16 <- socres.16[socres.16>=T.16]
hist(socres.16)




mean(socres.16)

## [1] 25.76066

var(socres.16)

## [1] 9.493101

kable(table(socres.16))
socres.16
 
Freq
 
23
 
2026
 
24
 
1541
 
25
 
1055
 
26
 
843
 
27
 
648
 
28
 
467
 
29
 
363
 
30
 
251
 
31
 
184
 
32
 
134
 
33
 
86
 
34
 
89
 
35
 
49
 
36
 
32
 
37
 
19
 
38
 
10
 
39
 
9
 
40
 
4
 
41
 
4
 
42
 
5
 
43
 
3
 
44
 
2
 
45
 
3
 
51
 
1
 
54
 
2
 
# 3mer
seq.v <- unlist(strsplit(shuffle.3,""))
len <- length(seq.v)
scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
for(j in 1:length(word.list)){
  scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
}
scores.temp <- as.vector(scores.temp)
scores.list[[i]] <- scores.temp

socres.16 <- unlist(scores.list)
T.16 <- quantile(socres.16,.999)
socres.16 <- socres.16[socres.16>=T.16]
hist(socres.16)




mean(socres.16)

## [1] 25.66886

var(socres.16)

## [1] 8.858734

kable(table(socres.16))
socres.16
 
Freq
 
23
 
930
 
24
 
732
 
25
 
553
 
26
 
378
 
27
 
288
 
28
 
220
 
29
 
168
 
30
 
118
 
31
 
77
 
32
 
56
 
33
 
34
 
34
 
21
 
35
 
19
 
36
 
10
 
37
 
12
 
38
 
6
 
39
 
8
 
40
 
4
 
41
 
2
 
42
 
1
 
43
 
2
 
44
 
2
 
45
 
1
 
# 4mer
seq.v <- unlist(strsplit(shuffle.4,""))
len <- length(seq.v)
scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
for(j in 1:length(word.list)){
  scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
}
scores.temp <- as.vector(scores.temp)
scores.list[[i]] <- scores.temp

socres.16 <- unlist(scores.list)
T.16 <- quantile(socres.16,.999)
socres.16 <- socres.16[socres.16>=T.16]
hist(socres.16)




mean(socres.16)

## [1] 24.91729

var(socres.16)

## [1] 4.985532

kable(table(socres.16))
socres.16
 
Freq
 
23
 
48
 
24
 
25
 
25
 
19
 
26
 
15
 
27
 
9
 
28
 
7
 
29
 
2
 
30
 
3
 
31
 
3
 
32
 
1
 
33
 
1
 
# 5mer
seq.v <- unlist(strsplit(shuffle.5,""))
len <- length(seq.v)
scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
for(j in 1:length(word.list)){
  scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
}
scores.temp <- as.vector(scores.temp)
scores.list[[i]] <- scores.temp

socres.16 <- unlist(scores.list)
T.16 <- quantile(socres.16,.999)
socres.16 <- socres.16[socres.16>=T.16]
#table(socres.16)
hist(socres.16)




mean(socres.16)

## [1] 26.48

var(socres.16)

## [1] 8.757766

kable(table(socres.16))
socres.16
 
Freq
 
24
 
262
 
25
 
204
 
26
 
144
 
27
 
98
 
28
 
79
 
29
 
39
 
30
 
32
 
31
 
28
 
32
 
20
 
33
 
14
 
34
 
9
 
35
 
4
 
36
 
4
 
37
 
4
 
38
 
3
 
39
 
2
 
40
 
2
 
44
 
1
 
51
 
1
 
For DNA data:

lengthof word: 4

data <- readDNAStringSet("./sequence.dna.txt")
seq.name <-  names(data)
seq <- paste(data)
seq <- seq[1:50]

k <- 4
word.list <- NULL
for(i in 1:50){
  seq.v <- unlist(strsplit(seq[i],""))
  len <- length(seq.v)
  word.list <- c(word.list, sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")}))
}

word.list <- unique(word.list)

# seq.v <- unlist(strsplit(seq,""))
# letter.dict <- unique(seq)
#
# len <- length(seq.v)
# kmer.o <- sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")})

scores.list <- list()
i <- 1
  seq.v <- unlist(strsplit(seq[i],""))
  len <- length(seq.v)
  scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
  for(j in 1:length(word.list)){
    scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
  }
  scores.temp <- as.vector(scores.temp)
  scores.list[[i]] <- scores.temp

socres.4 <- unlist(scores.list)
T.4 <- quantile(socres.4,.999) #set a T
socres.4 <- socres.4[socres.4>=T.4]
table(socres.4)

## socres.4
## 26 27 28 29 30 31 32 33 36 
## 71 65 59 30  7  3  9  6  1

hist(socres.4)




mean(socres.4)

## [1] 27.68127

var(socres.4)

## [1] 3.066008

lengthof word: 6

# 12
k <- 6
word.list <- NULL
for(i in 1:50){
  seq.v <- unlist(strsplit(seq[i],""))
  len <- length(seq.v)
  word.list <- c(word.list, sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")}))
}

word.list <- unique(word.list)

# seq.v <- unlist(strsplit(seq,""))
# letter.dict <- unique(seq)
#
# len <- length(seq.v)
# kmer.o <- sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")})

scores.list <- list()
i <- 1
seq.v <- unlist(strsplit(seq[i],""))
len <- length(seq.v)
scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
for(j in 1:length(word.list)){
  scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
}
scores.temp <- as.vector(scores.temp)
scores.list[[i]] <- scores.temp

socres.6 <- unlist(scores.list)
T.6 <- quantile(socres.6,.999)
socres.6<- socres.6[socres.12>=T.6]
table(socres.6)

## socres.6
##    -18    -17    -16    -15    -14    -13    -12    -11    -10     -9 
##      1     19    122    502   1616   4148   9135  17839  30029  46708 
##     -8     -7     -6     -5     -4     -3     -2     -1      0      1 
##  65803  82215 101207 113498 121924 137414 146644 158432 172517 178802 
##      2      3      4      5      6      7      8      9     10     11 
## 179153 177808 173615 162771 163037 151507 145695 133913 121480 109037 
##     12     13     14     15     16     17     18     19     20     21 
##  98724  88369  76383  67296  55174  49082  42024  34643  27584  22694 
##     22     23     24     25     26     27     28     29     30     31 
##  18865  15985  12269   9021   6915   5965   4384   3356   2229   1766 
##     32     33     34     35     36     37     38     39     40     41 
##   1335   1054    637    405    278    238    173     94     47     37 
##     42     43     44     45     46     47     48 
##     38     13      3      2      6      4      1

hist(socres.6)




mean(socres.6)

## [1] 4.324427

var(socres.6)

## [1] 62.15408

lengthof word: 8

# 
# #16
# k <- 8
# word.list <- NULL
# for(i in 1:50){
#   seq.v <- unlist(strsplit(seq[i],""))
#   len <- length(seq.v)
#   word.list <- c(word.list, sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")}))
# }
# 
# word.list <- unique(word.list)
# 
# # seq.v <- unlist(strsplit(seq,""))
# # letter.dict <- unique(seq)
# # 
# # len <- length(seq.v)
# # kmer.o <- sapply(1:(len-k+1), function(x){paste(seq.v[x:(x+k-1)],collapse="")})
# 
# scores.list <- list()
# i <- 1
# seq.v <- unlist(strsplit(seq[i],""))
# len <- length(seq.v)
# scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
# for(j in 1:length(word.list)){
#   scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
# }
# scores.temp <- as.vector(scores.temp)
# scores.list[[i]] <- scores.temp
# 
# socres.8 <- unlist(scores.list)
# T.8 <- quantile(socres.8,.999)
# socres.16 <- socres.8[socres.8>=T.8]
# table(socres.8)
# hist(socres.8)
# mean(socres.8)
# var(socres.8)

shuffle

# ## shuffle
# shuffle.1 <- seq.Rand(seq[1],1,2018)$randseq
# shuffle.2 <- seq.Rand(seq[1],2,2018)$randseq
# shuffle.3 <- seq.Rand(seq[1],3,2018)$randseq
# shuffle.4 <- seq.Rand(seq[1],4,2018)$randseq
# shuffle.5 <- seq.Rand(seq[1],5,2018)$randseq
# 
# # 1mer
# seq.v <- unlist(strsplit(shuffle.1,""))
# len <- length(seq.v)
# scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
# for(j in 1:length(word.list)){
#   scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
# }
# scores.temp <- as.vector(scores.temp)
# scores.list[[i]] <- scores.temp
# 
# socres.8 <- unlist(scores.list)
# T.8 <- quantile(socres.8,.999)
# socres.8 <- socres.8[socres.8>=T.8]
# hist(socres.8)
# mean(socres.8)
# var(socres.8)
# kable(table(socres.8))
# 
# # 2mer
# seq.v <- unlist(strsplit(shuffle.2,""))
# len <- length(seq.v)
# scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
# for(j in 1:length(word.list)){
#   scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
# }
# scores.temp <- as.vector(scores.temp)
# scores.list[[i]] <- scores.temp
# 
# socres.8 <- unlist(scores.list)
# T.8 <- quantile(socres.8,.999)
# socres.8 <- socres.8[socres.8>=T.8]
# hist(socres.8)
# mean(socres.8)
# var(socres.8)
# kable(table(socres.8))
# 
# # 3mer
# seq.v <- unlist(strsplit(shuffle.3,""))
# len <- length(seq.v)
# scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
# for(j in 1:length(word.list)){
#   scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
# }
# scores.temp <- as.vector(scores.temp)
# scores.list[[i]] <- scores.temp
# 
# socres.8 <- unlist(scores.list)
# T.8 <- quantile(socres.8,.999)
# socres.8 <- socres.8[socres.8>=T.8]
# hist(socres.8)
# mean(socres.8)
# var(socres.8)
# kable(table(socres.8))
# 
# # 4mer
# seq.v <- unlist(strsplit(shuffle.4,""))
# len <- length(seq.v)
# scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
# for(j in 1:length(word.list)){
#   scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
# }
# scores.temp <- as.vector(scores.temp)
# scores.list[[i]] <- scores.temp
# 
# socres.8 <- unlist(scores.list)
# T.8 <- quantile(socres.8,.999)
# socres.8 <- socres.8[socres.8>=T.8]
# hist(socres.8)
# mean(socres.8)
# var(socres.8)
# kable(table(socres.8))
# 
# # 5mer
# seq.v <- unlist(strsplit(shuffle.5,""))
# len <- length(seq.v)
# scores.temp <- matrix(rep(0,length(word.list)*(len+1-k)),length(word.list))
# for(j in 1:length(word.list)){
#   scores.temp[j, ] <- sapply(1:(len+1-k),function(x){sum(diag(BLOSUM62[unlist(strsplit(word.list[j],"")), seq.v[x:(x+k-1)]]))})
# }
# scores.temp <- as.vector(scores.temp)
# scores.list[[i]] <- scores.temp
# 
# socres.8 <- unlist(scores.list)
# T.8 <- quantile(socres.8,.999)
# socres.8 <- socres.8[socres.8>=T.8]
# #table(socres.8)
# hist(socres.8)
# mean(socres.8)
# var(socres.8)
# kable(table(socres.8))

联系客服提交作业获取报价与时间?

最快2~12小时即可完成，用技术和耐心帮助客户高效高质量完成作业.

联系客服