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

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


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


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
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.


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


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)


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?


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案例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 <-