0605 Cloudy

K-mer是什么，举两个例子就知道了。

• 例1：在碱基长度为20的序列里，设置k=19, 那么就会产生两个k-mer

  
  

• 例2：在碱基长度为20的序列里，设置k=17，那么就会产生4个k-mer

  
  
• 同理，长度为L的序列中，根据k-mer的k对序列进行分割，就会产生(L-k+1)个k-mer

知道这个原理就不难根据k-mer来对全序列长度来进行预测了。但是在此之前我们得了解一个概念，就是Coverage。我们根据基因组分析 K-mer 第0回 随机生成fasta文件的内容生成一些序列来进行简单的验证和计算。

1.1 首先生成假想序列

• 首先生成一个长度为50的参照序列

out_f1 <- "sample32_ref.fasta"         
out_f2 <- "sample32_ngs.fasta"         
param_len_ref <- 50                    #指定参照序列长度
narabi <- c("A","C","G","T")           
param_composition <- c(22, 28, 28, 22) #指定A,C,G,T的比例
param_len_ngs <- 20                    #指定NGS短序列read序列的长度
param_num_ngs <- 10                    #指定NGS短序列read序列的断片数
param_desc <- "contig"                 #给fasta文件指定description的内容

• 运行包

library(Biostrings)                   

• 生成reference

set.seed(1010)                         
ACGTset <- rep(narabi, param_composition)#根据之前指定的比例进行生成
reference <- paste(sample(ACGTset, param_len_ref, replace=T), collapse="") 
reference <- DNAStringSet(reference)  
names(reference) <- param_desc     
reference    

• 根据reference生成短序列

s_posi <- sample(1:(param_len_ref-param_len_ngs), param_num_ngs, replace=T)#通过长序列的长度剪掉NGS短序列的长度得到一个值，从1到这个值为止随机取样N个数字，N等同于短序列Read的序列数。其实就是等于随机在长序列上寻找N歌段序列的切入点           
hoge <- NULL                         
for(i in 1:length(s_posi)){            
  hoge <- append(hoge, subseq(reference, start=s_posi[i], width=param_len_ngs))#把N段段序列保存到hoge里
}
fasta <- hoge                          #这个hoge就是我们需要的fasta文件

• 添加description，并且确认一下短序列

description <- paste(param_desc, s_posi, (s_posi+param_len_ngs-1), sep="_")
names(fasta) <- description           
fasta                                 

• 最后保存文件

writeXStringSet(reference, file=out_f1, format="fasta", width=50)
writeXStringSet(fasta, file=out_f2, format="fasta", width=50)

1.2 计算Coverage

在用k-mer预测全长之前要理解一个概念，就是coverage。什么？你不知道coverage是什么？！在这里短序列一共有10段，每段20bp，这个就相当于NGS的下机数据已拥有10 x 20 = 200bp,另外一方面由于长序列是50bp, 这就相当于把长序列覆盖了 200 / 50 =4 遍，也就是coverage等于4。这样说有没有理解coverage的意思了呢？