基因组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
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 158,117评论 4 360
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 66,963评论 1 290
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 107,897评论 0 240
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,805评论 0 203
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,208评论 3 286
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,535评论 1 216
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,797评论 2 311
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,493评论 0 197
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,215评论 1 241
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,477评论 2 244
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,988评论 1 258
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,325评论 2 252
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,971评论 3 235
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,055评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,807评论 0 194
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,544评论 2 271
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,455评论 2 266

推荐阅读更多精彩内容