Kmer是從測序數(shù)據(jù)中滑窗提取出的長度為k的寡聚核苷酸序列，可以評估基因組大小、雜合度、重復(fù)序列比例等。
在測序reads均勻分布的前提下，有以下公式：
基因組長度 = 總堿基數(shù) / 平均測序深度 = 總kmer數(shù) / 平均kmer深度
根據(jù)該公式，可以使用總kmer數(shù)和平均kmer深度估算基因組長度。K值使用滿足以下公式計算的最大奇數(shù)：4 ^ K / genome > 200
做kmer有幾種軟件：SOAPdenovo2的Kmerfreq模塊，jellyfish和KAT（K-mer Analysis Toolkit）
jellyfish：簡書1、簡書2和GitHub
產(chǎn)生的kmer頻數(shù)分布數(shù)據(jù)需要用R包GenomeScope和findGSE來進(jìn)行統(tǒng)計估計。(http://qb.cshl.edu/genomescope/)
這里我還是用一種后來出現(xiàn)但整合了jellyfish和其他分析的軟件KAT( is accomplished through an integrated and modified version of Jellyfish2's counting method)
https://doi.org/10.1093/bioinformatics/btw663
Github: https://github.com/TGAC/KAT
Document：https://kat.readthedocs.io/en/latest/using.html
安裝
- 從bioconda里安裝最新版本的kat：
bioconda install kat
#這里我用以下這種簡單快捷的方法：
conda install -c bioconda kat
conda activate
kat -h
- 從GitHub安裝：
需要先安裝很多依賴包，不然可能很多功能用不了：
-GCC V4.8+
-make
-autoconf V2.53+
-automake V1.11+
-libtool V2.4.2+
-pthreads (probably already installed)
-zlib
-Python V3.5+ with the tabulate, scipy, numpy and matplotlib packages and C API installed. This is optional but highly recommended, without python KAT functionality is limited: no plots, no distribution analysis, and no documentation.
-Sphinx-doc V1.3+ (Optional: only required for building the documentation.
然后是一系列的安裝編譯過程：
git clone https://github.com/TGAC/KAT.git
cd KAT
./build_boost.sh
./autogen.sh
./configure
make
二、運行
- 給定一個k值（kmer = 17，小基因組一般用17，大基因組才用默認(rèn)的27），產(chǎn)生不同kmers數(shù)量的直方圖
usage：kat hist [options] (<input>)+
options: -o 輸出文件 [kat.hist]；-t [1]線程數(shù)； -l [1] histogram的最低值；-h [10000] histogram的最高值；-m [27] kmer的長度
nohup kat hist -o KAT_result/kat.hist -t 16 <(gzip -d -c NGS_10G_?.fq.gz) 2> kat_hist.log & #不支持壓縮格式，需先解壓
nohup kat hist -m 17 -o KAT_result/kat.hist -t 16 NGS_10G_1.clean.fq NGS_10G_2.clean.fq & #kat1.hist.png
nohup kat hist -m 17 -o KAT_result/kat.hist -h 400 -t 16 NGS_10G_1.clean.fq NGS_10G_2.clean.fq & #kat2.hist.png，縮短x軸看得更清楚測序深度
- Kmer的GC覆蓋圖
usage: kat gcp (<input>)+
options: -o 輸出文件 [kat.hist]；-t [1]線程數(shù); -x [1] 當(dāng)創(chuàng)建污染矩陣時，用于gc數(shù)據(jù)的bins的數(shù)量; -y [1000] 當(dāng)創(chuàng)建污染矩陣時，用于coverage數(shù)據(jù)的bins的數(shù)量; -m [27] kmer的長度
nohup kat gcp -m 17 -o KAT_result/kat.gcp -t 12 NGS_10G_1.clean.fq NGS_10G_2.clean.fq &
SOAPdenovo2的Kmerfreq模塊
ls NGS_10G_?.clean.fq > NGS10.clean.lst
KmerFreq_AR -k 17 -t 12 -m 1 -p KmerFreq.10G -q 33 NGS10.clean.lst > NGS.10G.kmer.count 2>NGS.10G.kmerfreq.cerr
freq.cz和.freq.cz.len文件是用于Corrector_AR對reads進(jìn)行校正分析
Corrector_AR -k 17 -t 12 -l 3 -Q 33 KmerFreq.10G.freq.cz KmerFreq.10G.freq.cz.len NGS10.clean.lst >corr.cout 2>corr.cerr ##沒有校正的reads就不需要再跑下一步
ls NGS_10G_1.clean.fq.cor.pair_1.fq.gz NGS_10G_2.clean.fq.cor.pair_2.fq.gz > kmer.input.fil.data
KmerFreq_AR -k 17 -t 12 -m 1 -p KmerFreq.10G.cor -q 33 kmer.input.fil.data > NGS.10G.kmer.cor.count 2>NGS.10G.kmerfreq.cor.cerr
jellyfish和genomescope
#/usr/bin/jellyfish jellyfish 1.1.10
nohup /usr/bin/jellyfish count -m 17 -s 1G -t 12 -C -o jellyfish.out <(gzip -d -c NGS_10G_?.fq.gz) jellyfish.out.log &
nohup /usr/bin/jellyfish count -m 17 -s 1G -t 12 -C -o jellyfish.out NGS_10G_1.fq NGS_10G_2.fq 2> jellyfish.out.log &
/usr/bin/jellyfish histo -t 12 jellyfish.out_0 > jellyfish.histo
~/Software/genomescope/genomescope.R jellyfish.histo 17 149 genomescope.result
