在基因组分析 K-mer 第4回里我们成功估计到了和真实值很接近的长度。在这一回里想要介绍一下Coverage对预测精度的影响。有关Coverage的概念可以回过去看第2回的讲解。
在第4回中，我们首先虚拟生成了1000bp长度的长序列，然后根据这段长序列随机生成了200段长度为20bp的短序列，也就是Coverage为4的一组数据。最后得到的预测结果是原序列长度为842。
在这一期里我们通过增加coverage的数值，也就是增加短序列数量来看一下coverage对最后的结果会造成什么样的影响。

4.1 生成序列

关于代码有无法理解的地方可以参考第0回 随机生成fasta。
长序列生成和基因组分析 K-mer 第4回一样，短序列生成数量增加到500。

out_f1 <- "sample34_ref.fasta"       
out_f2 <- "sample34_ngs.fasta"        
param_len_ref <- 1000                  
narabi <- c("A","C","G","T")          
param_composition <- c(22, 28, 28, 22) 
param_len_ngs <- 20                   
param_num_ngs <- 500              # 短序列数量     
param_desc <- "contig"                   

library(Biostrings)                 

set.seed(1010)                       
ACGTset <- rep(narabi, param_composition)
param_composition
reference <- paste(sample(ACGTset, param_len_ref, replace=T), collapse="")
reference <- DNAStringSet(reference) 
names(reference) <- param_desc                                       

s_posi <- sample(1:(param_len_ref-param_len_ngs), param_num_ngs, replace=T)
hoge <- NULL                           
for(i in 1:length(s_posi)){           
  hoge <- append(hoge, subseq(reference, start=s_posi[i], width=param_len_ngs))
}
fasta <- hoge                       
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

4.2 设置k=10

in_f <- "sample34_ngs.fasta"           
out_f <- "hoge12.txt"                 
param_kmer <- 10                       #k-mer的k值为10

library(Biostrings)                    

fasta <- readDNAStringSet(in_f, format="fasta")

hoge <- oligonucleotideFrequency(fasta, width=param_kmer)
out <- colSums(hoge)                  

write.table(out, out_f, sep="\t", append=F, quote=F, row.names=T, col.names=F)
sum(out > 0)                           #计算出现频度为0次以上的k-mer种类

结果是988

> sum(out > 0)                          
[1] 988

真实数值是1000，已经非常接近了。怎么样，是不是体会到了coverage的重要性。