基因组分析 K-mer 第2回 估计50bp的全基因长度

0607 Cloudy
先定义一下标题里的“估计”两个字在这里是什么意思。根据什么估计什么。

• 1.根据NGS短序列数据进行K-mer分析
• 2.对原完整的长序列的总碱基数量(总长度）进行估计

这个步骤一般会用在de novo拼接前。

2.1 估计序列长度(k=2)

• 读取根据之前的文章内容生成的fasta文件
  注意：长序列reference的真实碱基长度为50。
  基因组分析 K-mer 第1回 理解K-mer和Coverage的基本概念

in_f <- "sample32_ngs.fasta"           
out_f <- "hoge7.txt"                   
param_kmer <- 2                        #指定k-mer的k值
library(Biostrings) 
fasta <- readDNAStringSet(in_f, format="fasta")

• 确认一下fasta文件

fasta

out <- oligonucleotideFrequency(fasta, width=param_kmer)      # 计算一下k=2也就是连续两个碱基的出现频度。
tmp <- cbind(names(fasta), out)        #第一列为ID、与out结合形成新的数据
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#保存文件

• 这样我们就获得了每段序列的连续两个碱基出现的频度数据，查看一下tmp

• 我们再统计一下每个组合出现的频度，如下图。

colSums(out)

在所有的数据里，AA出现了1次，AC出现了11次，AG出现了24次......以此类推。
按照基础的数学排列组合逻辑来考虑，k=2的时候，ATGC四种碱基一共可以产生4的2次方也就是16个组合的k-mer，也就是上图中的列数。其中出现了一次以上的k-mer是14种。也就是说当k=2的时候产生的16种kmer里，14/16也就是说几乎每一种组合都在一次以上。
其实通过k-mer估计基因序列的原理其实就是把k值设定成比较大的数字，然后统计出现了1次以上k-mer的种类。我们的reference的真实长度是50，所以k=2并不能满足我们的需求去估计16以上的序列。估计结果和真实结果相差比较远。

2.2 估计序列长度(k=3)

• 我们把k值设定成3，继续尝试。
  当k=3的时候，理论上一共会出现4的3次方一共就是64种k-mer,这个数值已经大于reference的真实长度50,可以比 k=2更加贴切的估计出真实长度。

param_kmer=3
hoge <- oligonucleotideFrequency(fasta, width=param_kmer)
out <- colSums(hoge)
sum(out>0)

• 返回结果是29, 虽然离50还差了点距离，但是比起K=2的14，已经有了明显的改善。当然在真实的世界里，都是Gb单位的长度，这种几十几百的bp都只能算是误差范围。

> sum(out>0)
[1] 29

所以之后我们会用稍微长一点的序列来进行测试。