基因组survey——K-mer频谱

Kmer是从测序数据中滑窗提取出的长度为k的寡聚核苷酸序列,可以评估基因组大小、杂合度、重复序列比例等。
在测序reads均匀分布的前提下,有以下公式:
基因组长度 = 总碱基数 / 平均测序深度 = 总kmer数 / 平均kmer深度
根据该公式,可以使用总kmer数和平均kmer深度估算基因组长度。K值使用满足以下公式计算的最大奇数:4 ^ K / genome > 200

做kmer有几种软件:SOAPdenovo2的Kmerfreq模块,jellyfish和KAT(K-mer Analysis Toolkit)

jellyfish:简书1简书2GitHub
产生的kmer频数分布数据需要用R包GenomeScopefindGSE来进行统计估计。(http://qb.cshl.edu/genomescope/)

这里我还是用一种后来出现但整合了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

安装

  1. 从bioconda里安装最新版本的kat:
bioconda install kat
#这里我用以下这种简单快捷的方法:
conda install -c bioconda kat
conda activate
kat -h
  1. 从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

二、运行

  1. 给定一个k值(kmer = 17,小基因组一般用17,大基因组才用默认的27),产生不同kmers数量的直方图
    usage:kat hist [options] (<input>)+
    options: -o 输出文件 [kat.hist];-t [1]线程数; -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轴看得更清楚测序深度
kat1.hist.png

kat2.hist.png
  1. Kmer的GC覆盖图
    usage: kat gcp (<input>)+
    options: -o 输出文件 [kat.hist];-t [1]线程数; -x [1] 当创建污染矩阵时,用于gc数据的bins的数量; -y [1000] 当创建污染矩阵时,用于coverage数据的bins的数量; -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 &
kat.gcp.mx.png

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进行校正分析

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
plot.png

推荐阅读更多精彩内容