宏基因组学习笔记

分析策略

以下三个图来自微信公众号[希研未来]。

• 基于reads

• 基于组装

• 基于Bin

软件安装

常用的这些软件用Conda或者是mamba都很好处理，但是metaWRAP这个软件是真的难安装啊。踩坑多次终于安装上了：

mamba create -y --name metawrap-env --channel ursky metawrap-mg=1.3.2     

数据配置

CheckM

mkcd ~/database/checkm
axel -n 50 ftp://download.nmdc.cn/tools/meta/checkm/checkm_data_2015_01_16.tar.gz      
tar -xvf checkm_data_2015_01_16.tar.gz   

可以通过运行checkm data setRoot设置数据库路径。

kraken2

mkcd ~/database/kraken2
axel -n 50 ftp://download.nmdc.cn/tools/meta/kraken2/k2_pluspf_20230314.tar.gz
tar -xvf k2_pluspf_20230314.tar.gz

NCAI-nt

mkcd ~/database/ncbi.nt
~/.aspera/connect/bin/ascp -i ~/mambaforge/envs/tools4bioinf/etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l6000m \ 
    anonftp@ftp.ncbi.nlm.nih.gov:blast/db/v4/nt_v4.{00..85}.tar.gz ./

NCBI-taxonomy

mkcd ~/database/ncbi.taxonomy
axel -n 50 ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
tar -xvf taxdump.tar.gz

宿主基因组

通常第一步是去除宿主基因组，需要构建索引，（其他宿主类似），如下：

mkcd ~/database/human.genome
axel -n 50 http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz   
gunzip hg38.fa.gz
cat *fa > hg38.fa
rm chr*.fa
bmtool -d hg38.fa -o hg38.bitmask
srprism mkindex -i hg38.fa -o hg38.srprism -M 100000

为了识别不同的宿主，该软件参数-x可以指定宿主，如-x rice，默认的宿主是human，文件前缀是hg38.

修改配置文件

vi ~/mambaforge/envs/metawrap-env/bin/config-metawrap
# Paths to metaWRAP scripts (dont have to modify)
mw_path=$(which metawrap)
bin_path=${mw_path%/*}
SOFT=${bin_path}/metawrap-scripts
PIPES=${bin_path}/metawrap-modules

# CONFIGURABLE PATHS FOR DATABASES (see 'Databases' section of metaWRAP README for details)
# path to kraken standard database
# KRAKEN_DB=~/KRAKEN_DB
KRAKEN2_DB=~/database/kraken2.plus.plantandfungi 

# path to indexed human (or other host) genome (see metaWRAP website for guide). This includes .bitmask and .srprism files
BMTAGGER_DB=~/database/human.genome
#
# paths to BLAST databases
BLASTDB=~/database/ncbi.nt
TAXDUMP=~/database/ncbi.taxonomy

数据质控

如果宿主不是人类，可以使用参数--skip-bmtagger跳过：

metawrap read_qc --skip-bmtagger -1 data/C3_1.R1.raw.fastq -2 data/C3_1.R2.raw.fastq -t 60 -o readqc/C3_1

有多个文件时，批量运行：

for F in RAW_READS/*_1.fastq; do
R=${F%_*}_2.fastq
BASE=${F##*/}
SAMPLE=${BASE%_*}
metawrap read_qc -1 $F -2 $R -t 1 -o
READ_QC/$SAMPLE &
done

测试两个文件差不多35G，60线程跑了45分钟左右。

宏基因组组装

可以单个样品分别组装，也可以将群落所有文件一起组装。如果要全部一起组装，就先把文件合并：

cat CLEAN_READS/ERR*_1.fastq > CLEAN_READS/ALL_READS_1.fastq
cat CLEAN_READS/ERR*_2.fastq > CLEAN_READS/ALL_READS_2.fastq

开始组装：

metawrap assembly -1 readqc/C3_1/C3_1.R1.raw_val_1.fq -2 readqc/C3_1/C3_1.R2.raw_val_2.fq -m 100 -t 60 --megahit -o assembly  

有两种组装方式可以选择：

• megahit：默认的算法，速度更快。
• metaspades：速度更慢且需要更多的内存，但是结果更准确。

组装完成的结果文件是final_assembly.fasta，对应的报告文件是assembly_report.html.

测试两个文件差不多35G，60线程+100G内存跑了3小时20分钟左右。组装大小是108M，一共有67463个contigs.

分类和丰度

使用kraken2的，同时对原始的reads和组装后的contigs进行分析，但是目的不一样：

• 基于reads：是看样品中哪些群落被检测到了；
• 基于contigs：目的是看那些物种比其他的组装得更好，通常不适用contigs来推断整个群落组成结构。

metawrap kraken2  --no-preload -t 60 -s 1000000 -o kraken readqc/C3_1/C3_1.R* assembly/final_assembly.fasta

结果文件：

• .kraken是估计的分类特征表；
• .krona是KronaTools输出的描述统计信息；
• kronagram.html包含了所有的分类信息。

分箱Bin

使用三种算法对组装的结果进行分箱。这一步对文件的后缀貌似非常严格，好像必须是fastq. 耗时1小时10分钟。

metawrap binning -o init.binning -t 60 -a assembly/final_assembly.fasta --metabat2 --maxbin2 --concoct readqc/C3_1/C3_1.R1.raw_val_*fastq

合并分箱结果

这一步需要使用CheckM这个数据库，见数据配置环节。

metawrap bin_refinement -o bin.refinement -t 60 -A init.binning/metabat2_bins -B init.binning/maxbin2_bins -C init.binning/concoct_bins -c 50 -x 10  

运行时间差不多3小时。

分箱可视化

运行：

metawrap blobology -a assembly/final_assembly.fasta -t 60 -o blobology --bins bin.refinement/metawrap_50_10_bins readqc/C3_1/C3_1.R1.raw_val_*fastq  

遇到报错：

BLAST Database error: Error: Not a valid version 4 database.

解决方法：

ln -s /home/xxx/mambaforge/envs/blast+/bin/blastn /home/xxx/mambaforge/envs/metawrap-env/bin/blastn    

运行时间差不多10小时。

基因组草图丰度

metawrap quant_bins -b bin.refinement/metawrap_50_10_bins -o quant.bins -a assembly/final_assembly.fasta readqc/C3_1/C3_1.R1.raw_val_*fastq    

耗时1小时。

结果文件存放在bin_abundance_table.tab这个文件中。

重新组装bins

相当于对最佳的bin再次进行重新组装，只有让bins质量变得更好的reads才会被添加到bins中。这一步是用所有的reads进行分析的。

metawrap reassemble_bins -o bin.reassembly -1 readqc/C3_1/C3_1.R1.raw_val_1.fastq -2 readqc/C3_1/C3_1.R1.raw_val_2.fastq -t 60 -m 100 -c 50 -x 10 -b bin.refinement/metawrap_50_10_bins

耗时30分钟。

分箱物种丰度

耗时1小时13分钟。

metawrap classify_bins -b bin.reassembly/reassembled_bins -o bins.classification -t 60    

分箱功能注释

metawrap annotate_bins -o bins.anno -t 60 -b bin.reassembly/reassembled_bins  

遇到报错：

Please rename your contigs or use --centre XXX to generate clean contig names.

原因是bin中的序列名称超过20个字符了，因此需要修改序列名称。解决方法（见Prokka官网）：

写了个python代码解决这个问题：

import sys

fasta_file = sys.argv[1]
out_file = sys.argv[2]
index_file = sys.argv[3]

i = 0  # 从数字0开始

with open(fasta_file, "r") as f, open(out_file, "w") as out, open(index_file, "w") as index:
    for line in f:
        if line.startswith(">"):
            i += 1
            # 修改序列名称
            out.write(">seq{0}\n".format(i))
            # 输出新旧ID对应关系
            old_id = line.replace(">","")
            new_id = ">seq" + str(i)
            index.write(str(new_id) + "\t" + old_id)
        else:
            out.write(line)
    
print("Done.")

python3 ~/scripts/rename.fasta.id.py reassembled_bins/bin.1.orig.fa renamed.bins/bin1.fa bin1.id.txt

重新运行注释代码：

metawrap annotate_bins -o bins.anno -t 60 -b bin.reassembly/renamed.bins

运行很快，基本上就是一分钟左右。

• 功能注释：

既然拿到了gff文件，那么就可以计算不同基因在不同样品的中丰度了。

• 翻译的蛋白：

既然拿到了蛋白序列，那就可以用eggNOG-mapper进行功能注释了（参考博客）：

python3 ~/software/eggNOgmapper/emapper.py -i bin_translated_genes/bin1.faa --output eggnog -m diamond --cpu 60

基本上需要的信息都有了。

实战

数据质控

for F in 0.data/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc --skip-bmtagger -1 $F -2 $R -t 70 -o 1.read.qc/$SAMPLE & done

将质控后的数据移动到新的文件中：

mkdir -p 2.clean.reads/all
for i in 1.read.qc/*; do
b=${i#*/}
mv ${i}/final_pure_reads_1.fastq 2.clean.reads/${b}_1.fastq
mv ${i}/final_pure_reads_2.fastq 2.clean.reads/${b}_2.fastq
done

宏基因组组装

将F和R分别合并后再进行组装。

• 先合并

cat 2.clean.reads/*_1.fastq > 2.clean.reads/all/all_reads_1.fastq
cat 2.clean.reads/*_2.fastq > 2.clean.reads/all/all_reads_2.fastq

• 开始组装

nohup metawrap assembly -1 2.clean.reads/all/all_reads_1.fastq -2 2.clean.reads/all/all_reads_2.fastq -m 100 -t 60 --megahit -o 3.assembly &  

kraken2注释

nohup metawrap kraken2 -o 4.kraken2 -s 1000000 2.clean.reads/*.fastq 3.assembly/final_assembly.fasta &

分箱

metawrap binning -o 5.init.binning -t 60 -a 3.assembly/final_assembly.fasta --metabat2 --maxbin2 --concoct 2.clean.reads/*.fastq  

合并分箱结果

nohup metawrap bin_refinement -o 6.bin.refinement -t 60 -A 5.init.binning/metabat2_bins -B 5.init.binning/maxbin2_bins -C 5.init.binning/concoct_bins -c 50 -x 10  &

分箱结果可视化

nohup metawrap blobology -a 3.assembly/final_assembly.fasta -t 60 -o 7.blobology --bins 6.bin.refinement/metawrap_50_10_bins 2.clean.reads/*fastq  & 

分箱丰度

nohup metawrap quant_bins -b 6.bin.refinement/metawrap_50_10_bins -o 8.quant.bins -a 3.assembly/final_assembly.fasta 2.clean.reads/*fastq & 

再次分箱

nohup metawrap reassemble_bins -o 9.bin.reassembly -1 2.clean.reads/all/all_reads_1.fastq -2 2.clean.reads/all/all_reads_2.fastq -t 60 -m 100 -c 50 -x 10 -b 6.bin.refinement/metawrap_50_10_bins &

分箱物种丰度

metawrap classify_bins -b  9.bin.reassembly/reassembled_bins -o 10.bins.classification -t 60

分箱功能注释

mkdir 9.bin.reassembly/renamed.bins
mkdir 9.bin.reassembly/renamed.bins.id
for i in 9.bin.reassembly/reassembled_bins/*; do python3 ~/scripts/rename.fasta.id.py ${i} 9.bin.reassembly/renamed.bins/${i##*/} 9.bin.reassembly/renamed.bins.id/${i##*/}.id.txt; done

 nohup metawrap annotate_bins -o 11.bins.anno -t 60 -b 9.bin.reassembly/renamed.bins &   

eggNOG-mapper注释

mkdir 12.eggNOG-mapper
for i in 11.bins.anno/bin_translated_genes/*; do python3 ~/software/eggNOgmapper/emapper.py -m diamond --cpu 60 -i  11.bins.anno/bin_translated_genes/${i##*/} --output 12.eggNOG-mapper/${i##*/}.eggnogmapper.txt; done

所有脚本

for F in 0.data/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc --skip-bmtagger -1 $F -2 $R -t 70 -o 1.read.qc/$SAMPLE & done

mkdir -p 2.clean.reads/all

for i in 1.read.qc/*; do
b=${i#*/}
mv ${i}/final_pure_reads_1.fastq 2.clean.reads/${b}_1.fastq
mv ${i}/final_pure_reads_2.fastq 2.clean.reads/${b}_2.fastq
done

cat 2.clean.reads/*_1.fastq > 2.clean.reads/all/all_reads_1.fastq
cat 2.clean.reads/*_2.fastq > 2.clean.reads/all/all_reads_2.fastq

metawrap assembly -1 2.clean.reads/all/all_reads_1.fastq -2 2.clean.reads/all/all_reads_2.fastq -m 100 -t 60 --megahit -o 3.assembly

metawrap kraken2 -o 4.kraken2 -s 1000000 2.clean.reads/*.fastq 3.assembly/final_assembly.fasta &      

metawrap binning -o 5.init.binning -t 60 -a 3.assembly/final_assembly.fasta --metabat2 --maxbin2 --concoct 2.clean.reads/*.fastq 

metawrap bin_refinement -o 6.bin.refinement -t 60 -A 5.init.binning/metabat2_bins -B 5.init.binning/maxbin2_bins -C 5.init.binning/concoct_bins -c 50 -x 10

metawrap blobology -a 3.assembly/final_assembly.fasta -t 60 -o 7.blobology --bins 6.bin.refinement/metawrap_50_10_bins 2.clean.reads/*fastq

metawrap quant_bins -b 6.bin.refinement/metawrap_50_10_bins -o 8.quant.bins -a 3.assembly/final_assembly.fasta 2.clean.reads/*fastq

metawrap reassemble_bins -o 9.bin.reassembly -1 2.clean.reads/all/all_reads_1.fastq -2 2.clean.reads/all/all_reads_2.fastq -t 60 -m 100 -c 50 -x 10 -b 6.bin.refinement/metawrap_50_10_bins

metawrap classify_bins -b  9.bin.reassembly/reassembled_bins -o 10.bins.classification -t 60 

mkdir 9.bin.reassembly/renamed.bins
mkdir 9.bin.reassembly/renamed.bins.id
for i in 9.bin.reassembly/reassembled_bins/*; do python3 ~/scripts/rename.fasta.id.py ${i} 9.bin.reassembly/renamed.bins/${i##*/} 9.bin.reassembly/renamed.bins.id/${i##*/}.id.txt; done

metawrap annotate_bins -o 11.bins.anno -t 60 -b 9.bin.reassembly/renamed.bins

mkdir 12.eggNOG-mapper
for i in 11.bins.anno/bin_translated_genes/*; do python3 ~/software/eggNOgmapper/emapper.py -m diamond --cpu 60 -i  11.bins.anno/bin_translated_genes/${i##*/} --output 12.eggNOG-mapper/${i##*/}.eggnogmapper.txt; done