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的意思了呢？
