0608 Cloudy
在第二回的文章里我们用k-mer法预测了比较短的序列长度。因为序列本身较短，再加上是预测上的误差导致结果差强人意，在这一回里我会尝试用一毛一样的方法来预测1000bp的序列。废话不多说，直接来操作吧。

3.1 生成1000bp的长序列

代码有不懂的地方请参照前两回的讲解。主要做了以下两件事

• 生成一段1000bp的长序列
• 根据长序列(以长序列为参照序列)生成200段20bp的短序列

• 复习
  根据基因组分析 K-mer 第1回 理解K-mer和Coverage的基本概念我们可以知道这组数据的coverage是4

out_f1 <- "sample33_ref.fasta"        
out_f2 <- "sample33_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 <- 200                  
param_desc <- "contig"                    
library(Biostrings)                    
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                    

• 长序列确认:

  
  

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)

3.2 用短序列来预测原长序列的长度

方法和第二回的内容非常相似。首先用K=3来尝试一下。

3.2.1 k=3的k-mer

in_f <- "sample33_ngs.fasta"           
out_f <- "hoge10.txt"                  
param_kmer <- 3                        

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)

打开out你可以看到各k-mer出现的频度，如下图

如同第二回所述，出现0次以上的k-mer的种类数就是原参照序列也就是长序列的长度。在此，所有的k-mer都在0次以上,一共有64个k-mer,所以原长度是64?
No, No, No,这个问题在第二回的讲解里也出现过，K=3的话撑死了也就64个k-mer, 远远不能满足1000bp长度的基因，所以，需要选择更大的k值。

3.2.2 k=10, start again

把k设置成10，再来试一次

in_f <- "sample33_ngs.fasta"          
out_f <- "hoge11.txt"                  
param_kmer <- 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)

通过产看k-mer的理论上的种类可以看出就是4的10次方。

length(out)        
[1] 1048576

其中大于0的k-mer是842,也就是推算出原基因长度就是842，怎么样，是不是已经很接近真实答案1000了。

sum(out > 0)   
 [1] 842