0607 Cloudy
先定義一下標(biāo)題里的“估計(jì)”兩個(gè)字在這里是什么意思街图。根據(jù)什么估計(jì)什么。
- 1.根據(jù)NGS短序列數(shù)據(jù)進(jìn)行K-mer分析
- 2.對(duì)原完整的長(zhǎng)序列的總堿基數(shù)量(總長(zhǎng)度）進(jìn)行估計(jì)
這個(gè)步驟一般會(huì)用在de novo拼接前懒构。
2.1 估計(jì)序列長(zhǎng)度(k=2)
- 讀取根據(jù)之前的文章內(nèi)容生成的fasta文件
注意：長(zhǎng)序列reference的真實(shí)堿基長(zhǎng)度為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")
- 確認(rèn)一下fasta文件
fasta
out <- oligonucleotideFrequency(fasta, width=param_kmer) # 計(jì)算一下k=2也就是連續(xù)兩個(gè)堿基的出現(xiàn)頻度耘擂。
tmp <- cbind(names(fasta), out) #第一列為ID胆剧、與out結(jié)合形成新的數(shù)據(jù)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#保存文件
- 這樣我們就獲得了每段序列的連續(xù)兩個(gè)堿基出現(xiàn)的頻度數(shù)據(jù)，查看一下tmp
- 我們?cè)俳y(tǒng)計(jì)一下每個(gè)組合出現(xiàn)的頻度醉冤，如下圖秩霍。
colSums(out)
在所有的數(shù)據(jù)里，AA出現(xiàn)了1次蚁阳，AC出現(xiàn)了11次铃绒，AG出現(xiàn)了24次......以此類推。
按照基礎(chǔ)的數(shù)學(xué)排列組合邏輯來(lái)考慮螺捐，k=2的時(shí)候颠悬，ATGC四種堿基一共可以產(chǎn)生4的2次方也就是16個(gè)組合的k-mer，也就是上圖中的列數(shù)定血。其中出現(xiàn)了一次以上的k-mer是14種赔癌。也就是說(shuō)當(dāng)k=2的時(shí)候產(chǎn)生的16種kmer里，14/16也就是說(shuō)幾乎每一種組合都在一次以上澜沟。
其實(shí)通過(guò)k-mer估計(jì)基因序列的原理其實(shí)就是把k值設(shè)定成比較大的數(shù)字灾票，然后統(tǒng)計(jì)出現(xiàn)了1次以上k-mer的種類。我們的reference的真實(shí)長(zhǎng)度是50茫虽，所以k=2并不能滿足我們的需求去估計(jì)16以上的序列刊苍。估計(jì)結(jié)果和真實(shí)結(jié)果相差比較遠(yuǎn)。
2.2 估計(jì)序列長(zhǎng)度(k=3)
- 我們把k值設(shè)定成3濒析，繼續(xù)嘗試正什。
當(dāng)k=3的時(shí)候，理論上一共會(huì)出現(xiàn)4的3次方一共就是64種k-mer,這個(gè)數(shù)值已經(jīng)大于reference的真實(shí)長(zhǎng)度50,可以比 k=2更加貼切的估計(jì)出真實(shí)長(zhǎng)度悼枢。
param_kmer=3
hoge <- oligonucleotideFrequency(fasta, width=param_kmer)
out <- colSums(hoge)
sum(out>0)
- 返回結(jié)果是29, 雖然離50還差了點(diǎn)距離埠忘，但是比起K=2的14，已經(jīng)有了明顯的改善馒索。當(dāng)然在真實(shí)的世界里莹妒，都是Gb單位的長(zhǎng)度，這種幾十幾百的bp都只能算是誤差范圍绰上。
> sum(out>0)
[1] 29
所以之后我們會(huì)用稍微長(zhǎng)一點(diǎn)的序列來(lái)進(jìn)行測(cè)試旨怠。
