3-跟着science学习宏基因组-序列比对组装

写在前面

宏基因组相信会有很多人都会对他产生兴趣，但是却不像扩增子那样每个人都可以在自己的电脑上运行宏基因组数据。我们这份教程大部分都可以在帮你基本上运行了，只是有的数据库实在是太大了。例如：nr，所以大家还是尽量下载建好的索引。

温馨提示: 有的步骤我标记了step，代表的是原流程的step，大家可以对照原流程查看。例如：step10。

序列比对部分

使用megahit 组装

megahit的安装也挺简单，直接使用conda就可以。安装megahit；

conda instal megahit

这里使用去除宿主序列，采用双端模式进行序列组装。

#-使用去除宿主的序列组装
megahit -1 ./host_filtering/raw/SUBERR793599_R2_filtered.fastq -2 ./host_filtering/raw/SUBERR793599_R2_filtered.fastq -o ./megagta/SUBERR793599 -t 6

# 使用对其的序列
megahit -1 ./host_filtering/raw/SUBERR793599_R1_filtered.fastq_pairs_R1.fastq -2 ./host_filtering/raw/SUBERR793599_R2_filtered.fastq_pairs_R2.fastq -o ./megahit/SUBERR793599 -t 6

  • t 指定线程数量，推荐32线程以上；

本次测试数据仅仅用16s即可跑完一个；

批量组装

注意，如果单个样本运行过了，要删除一下输出文件夹。防止megahit报错。

mkdir megahit

for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
megahit -1 ./host_filtering/raw/${base}_R1_filtered.fastq_pairs_R1.fastq -2 ./host_filtering/raw/${base}_R2_filtered.fastq_pairs_R2.fastq -o ./megahit/${base} -t 6
done

# tail -n+2 ./metaGroup.txt |cut -f 1|sed 's/^/temp\/qc\//;s/$/_R1_kneaddata_paired_1.fastq/'| tr '\n' ','|sed 's/,$//'`

step10 使用Diamond比对nr库

使用Diamond进行nr库比对

52G的nr蛋白库，解压缩，压缩文件55.8G，解压后142G；这个数据库比较难下载，建议大家耐心。

# gunzip -c nr.gz > nr.fa
#尝试使用pigz多线程解压缩
time unpigz -p 32 nr.gz

数据库下载完成后建立索引：diamond使用有限制，只能用于比对蛋白数据库。它的优势是处理>1M 的query，量越大速度越快。构建索引过程花费了442秒，使用了32核心。本服务器地址：~/db/nr.dmnd

nr库diamond索引构建

# conda install diamond# 安装工具
diamond makedb --in nr.fa -d nr -p 6

运行前端数据，具体的diamond使用可以参照地址：http://www.chenlianfu.com/?p=2703;本次测试需要时间1.5h。

#建立比对文件夹，不见文件夹会--daa选项会提示错误
mkdir diamond
# 开始比对
diamond blastx -c 1 --db ~/db/nr/db/nr.dmnd -t /tmp -p 6 -q ./trimmomatic/SUBERR793599_forward_paired.fq.gz --daa ./diamond/SUBERR793599.1.daa

  • t 临时文件夹位置

  • p 设置线程数量

  • db 数据库

  • q 输入文件

后端运行

diamond blastx -c 1 --db ~/db/nr/db/nr.dmnd -t /tmp -p 34 -q ./trimmomatic/SUBERR793599_reverse_paired.fq.gz --daa ./diamond/SUBERR793599.2.daa

大家注意翻译软件

例如有道云等软件会调用crtl + C ，导致运行终止，这里我鼠标移动了一下，就终止了，手动补上来。

# diamond blastx -c 1 --db ~/db/nr/db/nr.dmnd -t /tmp -p 20 -q ./trimmomatic/SUBERR793607_reverse_paired.fq.gz --daa ./diamond/SUBERR793607.2.daa

# diamond blastx -c 1 --db ~/db/nr/db/nr.dmnd -t /tmp -p 20 -q ./trimmomatic/SUBERR793608_forward_paired.fq.gz --daa ./diamond/SUBERR793608.1.daa

# diamond blastx -c 1 --db ~/db/nr/db/nr.dmnd -t /tmp -p 20 -q ./trimmomatic/SUBERR793608_reverse_paired.fq.gz --daa ./diamond/SUBERR793608.2.daa

批量运行

mkdir diamond

for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
# 开始比对
diamond blastx -c 1 --db ~/db/nr/db/nr.dmnd -t /tmp -p 6 -q ./trimmomatic/${base}_forward_paired.fq.gz --daa ./diamond/${base}.1.daa

diamond blastx -c 1 --db ~/db/nr/db/nr.dmnd -t /tmp -p 34 -q ./trimmomatic/${base}_reverse_paired.fq.gz --daa ./diamond/${base}.2.daa
done

附

step9 使用megagta组装

得到了clean的数据，这下我们需要进行组装：MegaGTA在Xander算法上进行了改进，并实现了更高的灵敏度，准确性和速度。而且，它能够从超大型宏基因组学数据集中组装基因序列。本次我修改了，是因为这个软件已经超过5年没有维护了，所以我就不用了。

安装

conda install megagta

这里我们去除了megagta工具的使用，实在是很难配置，相信大家也一样。我们做一个替换。https://www.ncbi.nlm.nih.gov/pubmed/29072142

megahit

mkdir megagta
megagta.py --continue -1 ./host_filtering/raw/SUBERR793599_R1_filtered.fastq_pairs_R1.fastq -2 ./host_filtering/raw/SUBERR793599_R2_filtered.fastq_pairs_R2.fastq -o ./megagta/SUBERR793599 -g gene_list.txt -t 6 -m 0.5 --min-contig-len 300
2> ./megagta/SUBERR793599/megagta.log

~/install/megagta/bin/post_proc.sh -g /mnt/data/home/NIOO/mattiash/install/megagta/share/RDPTools/Xander_assembler/gene_resource -d {params.outdir}contigs -m 16G -c 0.01

这部分测试失败了：但是我并没有继续寻找 原因，因为这个拼装软件已经有五年没有更新过了，所以我想也没有必要继续学习了。文中我也做了个替换。

根际互作生物学研究室 简介

根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。

团队工作及其成果 （点击查看）

了解 交流 合作

  • 团队成员邮箱 袁军：

    junyuan@njau.edu.cn；

    文涛：

    2018203048@njau.edu.cn

  • 团队公众号：

赞 (0)

相关推荐