宏基因组之基因丰度计算

目前有两种方式可计算宏基因组基因的丰度,一种是基于比对,比如bwa,bowtie,soapaligner等主流的比对软件,另一种是不比对快速估计基因丰度,可以用近俩年来流行的salmon软件,这个软件在转录组的数据比对中也经常用到,可以直接计算出原始的Counts值和标准化的TPM值,此外由于是基于非比对,计算的速度得到很大的提升,同时也节省了部分的内存(少了庞大的sam/bam文件),可以说是多快好省,但是目前的高分文章中也还是不少用基于比对的方法去计算宏基因组的基因丰度的,尤其是Soapaligner比对软件,是由华大公司开发的Soap系列的工具,首先在速度上具有较大的优势,同时也很好的解决了在小内存情况下将短序列比对到参考基因组上的问题,作用也不容小觑。当然也有不少研究人员在前期质控环节进行去宿主也是不错的选择~下面我就分别简单介绍一下基于比对的soapaligner和不比对快速估计的samlon俩个软件的操作流程!!

Option1 SOAPaligenr

图片.png

安装与比对


1. 安装  
conda install soap
2. 构建索引  
2bwt-builder contig.fa 
会产生13个索引文件
3.双端序列进行比对  
soap -p 6 -r 2 -m 200 -x 400 -M 4 -a ./cleandata/B1.clean.fq1.gz -b ./sample1/B1.clean.fq2.gz -D contig.index -o B1_PE.soap -2 B1_SE.soap -u B1_UN.fasta

这个软件的参数有点多,-a,-b分别输入双端序列, -D 索引文件的地址,-o,输出能比对上组装后的参考序列且reads成对存在的比对文件,-2,能比对上组装后的参考序列但reads不成对存在的比对文件;-u,不能比对上组装后的参考序列。至于-m,-x参数比较主观,看过很多文献都有不同的设定,这里贴出一篇文献供参考(DOI:10.3389/fmicb.2017.01546),文献中的参数为−m 4 −r 2 −m 100 −x 1000!!

标准化

看过十几篇文献在基因丰度标准化这一块都引用到了一篇2012年发表在Nature上的一篇文献(qin2012,doi:10.1038/nature11450),看过之后在该文章的补充材料里发现了计算基因丰度的步骤,如下图俩步骤,看过之后感觉其实和转录组中的RPKM/TPM等标准化的思路较为类似,先标准化基因的长度,然后再除以标准化的值的总和。还是挺简单的嘛。


图片.png

我们也按照这个思路根据soapaligner比对的结果来进行标准化,在结果文件中我们能得到三个文件,分别为:B1_PE.soap:能比对上组装后的参考序列且reads成对存在的比对文件;B1_SE.soap: 能比对上组装后的参考序列但reads不成对存在的比对文件;sample1_UN.fasta:不能比对上组装后的参考序列,通常为fasta格式;我们只需要比对上的文件B1_PE.soap和B1_SE.soap即可。这两个文件的第8列中比对上的基因,我们将其取出来sort|uniq -c 就可以计算出基因的原始copy数目,也就是比对上特定基因的数目。之后再计算出比对基因的长度,也就是你组装去冗余之后的非冗余基因集的各基因长度,这步可以写个小脚本或者用现成的轮子seqkit/seqtk都可以很快的解决,最后利用上面的公式进行标准化就行了~ 还需要注意的一点,没有比对上的基因我们就用0代替他的原始丰度就好了。

Option2 Samlon

这种方法就更加简单了,也是大致两步走,建立索引,序列比对,只不过时间将大大缩短,这种方法网上的教程有很多,我这里贴出我的代码大家可以比较一下使用~ 建议大家使用时候都安装最新版本的软件,github地址为:https://github.com/COMBINE-lab/salmon/releases/

#1.建立索引
salmon index -t B1_NR100nl.fasta -p 9 -k 31 -i ./index
#2.比对
salmon quant --validateMappings -i ./index -l A -p 3  --meta -1 ../clean1_${b}.fastq -2 ../clean2_${b}.fastq -o ${b}.quant #基因定量

注意的就是宏基因组数据加上-p meta这个参数,其他参数默认就好~,得到的结果有个quant.sf文件,丰度信息就在里面。为5列的信息,第4列就是我们想要的标准化后的基因丰度了

图片.png

OK,拿到基因丰度文件之后,我们就可以用来进行后续分析,比如计算eggnog数据库比对后的功能基因丰度~这个留到下次再说吧!!

下面是刚创建个人公众号,会定时更新R、linux、python,组学方面的学习内容,请多多支持呦~~


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