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
所以之后我们会用稍微长一点的序列来进行测试。
网友评论