如何读取DNA序列，并赋值计算

慢的原因是程序里每一轮循环里的赋值都会重新copy整个数组。
所以这也是为什么

• 好多有其他语言编程经验的小伙伴学R的时候会感觉循环特别慢。
• R里需要尽量避免使用循环。
• R里大部分函数都为向量化(vectorized)的函数，这是正确的姿势。

你的代码完全不需要循环,使用相同数据来比较一下：

library(seqinr)
library(Biostrings)
#devtools::install_github("collectivemedia/tictoc")
library(tictoc)

sampleseq="ACAAGATGCCATTGTCCCCCGGCCTCCTGCTGCTGCTGCTCTCCGGGGCCACGGCCACCGCTGCCCTGCCCCTGGAGGGTGGCCCCACCGGCCGAGACAGCGAGCATATGCAGGAAGCGGCAGGAATAAGGAAAAGCAGCCTCCTGACTTTCCTCGCTTGGTGGTTTGAGTGGACCTCCCAGGCCAGTGCCGGGCCCCTCATAGGAGAGGAAGCTCGGGAGGTGGCCAGGCGGCAGGAAGGCGCACCCCCCCAGCAATCCGCGCGCCGGGACAGAATGCCCTGCAGGAACTTCTTCTGGAAGACCTTCTCCTCCTGCAAATAAAACCTCACCCATGAATGCTCACGCAAGTTTAATTACAGACCTGAA"
seq_b <- Biostrings::DNAString(sampleseq)
seq<-seq_b

你的代码

tic("forloop used:")

l<-length(seq)
c<-rep(NA,l)
g<-rep(NA,l)
a<-rep(NA,l)
t<-rep(NA,l)
x<-rep(NA,l)
y<-rep(NA,l)
z<-rep(NA,l)
for(i in 1:l){
   cc<-letterFrequencyInSlidingView(seq,"C",view.width = 1)
   c [ i ]<-sum(cc[1:i]) # 计算前 i 个字母中c的个数，并赋值给c [ i ]
   gg<-letterFrequencyInSlidingView(seq,"G",view.width = 1)
   g[ i]<-sum(gg[1:i])
   aa<-letterFrequencyInSlidingView(seq,"A",view.width = 1)
   a[ i]<-sum(aa[1:i])
   tt<-letterFrequencyInSlidingView(seq,"T",view.width = 1)
   t[ i ]<-sum(tt[1:i])
   x [ i ]<-(a[ i]+g[ i])-(c [i]+t[ i])
   y [ i ]<-(a[ i]+c[ i ])-(g[ i]+t[ i])
   z [ i ]<-(a[ i]+t[ i])-(c[ i]+g[ i])
}

toc()


forloop used:: 0.67 sec elapsed

使用cumsum

library(tictoc)

tic("cumsum used: ")
seq_v <- seqinr::getSequence(sampleseq)
a2 <- cumsum(seq_v=="A")
c2 <- cumsum(seq_v=="C")
t2 <- cumsum(seq_v=="T")
g2 <- cumsum(seq_v=="G")

x2 <- a2 + g2 - c2 - t2
y2 <- a2 + c2 - g2 - t2
z2 <- a2 + t2 - c2 - g2

toc()

cumsum used: : 0.16 sec elapsed

结果for循环代码0.67s, 向量化代码0.16s

你可以试试使用更长的测试数据，差距会更大。测了个十倍，3000nt的序列，for循环变为3.12s，向量化代码0.1s

结果完全一样：

df <- data.frame(a=a,a2=a2,t=t,t2=t2,c=c,c2=c2,g=g,g2=g2,x=x,x2=x2,y=y,y2=y2,z=z,z2=z2)
dplyr::as_tibble(df)

# A tibble: 368 x 14
       a    a2     t    t2     c    c2     g    g2     x    x2     y    y2     z    z2
   <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
 1     1     1     0     0     0     0     0     0     1     1     1     1     1     1
 2     1     1     0     0     1     1     0     0     0     0     2     2     0     0
 3     2     2     0     0     1     1     0     0     1     1     3     3     1     1
 4     3     3     0     0     1     1     0     0     2     2     4     4     2     2
 5     3     3     0     0     1     1     1     1     3     3     3     3     1     1
 6     4     4     0     0     1     1     1     1     4     4     4     4     2     2
 7     4     4     1     1     1     1     1     1     3     3     3     3     3     3
 8     4     4     1     1     1     1     2     2     4     4     2     2     2     2
 9     4     4     1     1     2     2     2     2     3     3     3     3     1     1
10     4     4     1     1     3     3     2     2     2     2     4     4     0     0
# ... with 358 more rows