易宏基因组分析流程--记录

evekankan
2,207 阅读19分钟 
参考：易宏基因组分析：https://github.com/YongxinLiu/EasyMetagenome/blob/master/1Pipeline.sh

一、数据预处理

1、准备工作

# 创建工作目录
mkdir easy_metagenome && cd easy_metagenome
mkdir -p seq temp result
# 下载元数据，到result目录下
wget http://www.imeta.science/github/EasyMetagenome/result/metadata.txt
1.1 查看工作目录
# 查看当前准备的工作目录结构
tree easy_metagenome

image.png

1.2 检查元数据文件格式
# 检查文件格式，^I为制表符，$为Linux换行，^M$为Windows回车，^M为Mac换行符
cat -A result/metadata.txt

image.png

1.3 元数据细节优化
# 转换Windows回车为Linux换行
sed -i 's/\r//' result/metadata.txt
# 去除数据中的一个多余空格
sed -i 's/Male  /Male/' result/metadata.txt
cat -A result/metadata.txt

image.png

1.4 下载序列数据
cd seq/
awk '{system("wget -c http://www.imeta.science/github/EasyMetagenome/seq/"$1"_1.fq.gz")}' <(tail -n+2 ../result/metadata.txt)
awk '{system("wget -c http://www.imeta.science/github/EasyMetagenome/seq/"$1"_2.fq.gz")}' <(tail -n+2 ../result/metadata.txt)
cd ..
1.5、查看当前目录结构，检查数据下载结果
tree ./easy_metagenome

image.png

1.6、查看序列数据信息##### 1.6、查看序列数据信息
ls -lsh seq/*.fq.gz

image.png

1.7、查看序列文件格式
# 使用zless和zcat可以查看压缩文件内容
zless seq/C1_1.fq.gz | head -n4
zcat seq/C1_2.fq.gz | head -n4

2、FastQC质量评估报告

质量控制软件kneaddata本身就是一个宏基因组分析流程：
FastQC进行质量控制；
Trimmomatic进行过滤；
BowTie2进行去宿主；
2.1、启动质量控制环境kneaddata
# 启动conda虚拟环境--kneaddata
conda activate kneaddata
2.2、检查虚拟环境
fastqc --version    # FastQC v0.12.1
multiqc --version   # multiqc, version 1.16
kneaddata --version     # kneaddata v0.12.0
2.3、fastqc质量控制
# time统计运行时间，fastqc质量控制，-t 指定线程数
time fastqc seq/*.gz -t 2

image.png

2.4、查看当前目录结构
tree easy_metagenome

image.png

2.5、使用multiqc整理质控报告
# 整理fastqc报告，输出multiqc_report.html至result/qc目录
multiqc -d seq/ -o result/qc -f 
# 查看当前工作目录结构
tree ../easy_metagenome
# 查看mutiqc_report.html看查看质量控制报告信息

image.png

2.6、multiqc说明
MultiQC（Multiple Quality Control）是一个用于生物信息学和数据分析的开源工具；
旨在简化和改进数据质量控制报告的生成；
它可以自动整合和分析多种数据分析工具和流程的输出，生成易于理解的、交互式的报告；
帮助研究人员更好地理解和评估其实验数据的质量。

3、质量控制

3.1、创建fastp质量控制输出路径
mkdir -p temp/qc
3.2、fastp版本检查
fastp --version    # fastp 0.23.4
3.3、单样本质量控制
i=C1
fastp -i seq/${i}_1.fq.gz  -I seq/${i}_2.fq.gz \
  -j temp/qc/${i}_fastp.json -h temp/qc/${i}_fastp.html \
  -o temp/qc/${i}_1.fastq -O temp/qc/${i}_2.fastq
报错，缺少数据
上面在下载数据的步骤中少下载了数据，两行下载数据的代码，只运行了第一个，这里补充上，重新下载了。
# 下载完成后查看文件项目组成
tree ../easy_metagenomic

image.png

3.4、重新运行单样本质量控制
# 运行结果：
tree ../easy_metagenome

image.png

3.5 多样本并行
 tail -n+2 result/metadata.txt | cut -f1 | rush -j 2 \
  "fastp -i seq/{1}_1.fq.gz -I seq/{1}_2.fq.gz \
    -j temp/qc/{1}_fastp.json -h temp/qc/{1}_fastp.html \
    -o temp/qc/{1}_1.fastq  -O temp/qc/{1}_2.fastq"

# 运行结果查看
tree ../easy_metagenome

image.png

4、Kneaddata质量控制和去宿主

4.1、Kneaddata
Kneaddata本身就是一个流程，依赖trimmomatic质量控制和去接头，bowtie2比对宿主，然后筛选非宿主序列用于下游分析。
4.2、核心版本检查
Kneaddata --version    # kneaddata v0.12.0
trimmomatic -version    # version 2.2.5
bowtie2 --version    # 0.39
4.3、单样品质量控制
i=C1
kneaddata -i1 seq/${i}_1.fq.gz -i2 seq/${i}_2.fq.gz \
  -o temp/qc -v -t 2 --remove-intermediate-output \
  --trimmomatic ${soft}/envs/kneaddata/share/trimmomatic/ \
  --trimmomatic-options "ILLUMINACLIP:${soft}/envs/kneaddata/share/trimmomatic/adapters/TruSeq2-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50" \
  --reorder --bowtie2-options "--very-sensitive --dovetail" \
  -db ${db}/kneaddata/human_genome/hg37dec_v0.1
报错
报错信息：
Error message returned from Trimmomatic :
Error: Invalid or corrupt jarfile /home/users-data/liukai/anaconda3/envs/py20/share/trimmomatic/trimmomatic

问题原因：
由于biobakery工作流程和conda安装之间存在冲突，~/anaconda3/envs/kneaddata/bin/trimmomatic只是java文件的纯文本包装；
但是biobabery_workflows尝试将此纯文本调用为jar文件；
解决办法：
conda activate kneaddata
cd ~/anaconda3/envs/kneaddata/bin
ln -s ../share/trimmomatic/trimmomatic.jar . sed -i 's/trimmomatic_jar="trimmomatic\*"/trimmomatic_jar="trimmomatic.jar"/' ../lib/python3.9/site-p ackages/kneaddata/config.py 
python ../lib/python3.9/site-packages/kneaddata/config.py
4.4、查看运行结果
# 查看单样品控制运行结果
tree ../easy_metagenome

image.png

# 查看质量控制生成文件大小
ls -shtr temp/qc/${i}_1_kneaddata_paired_?.fastq

image.png

4.5、多样品并行质量控制
# 使用rush进行并行管理
tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \
  "kneaddata -i1 seq/{1}_1.fq.gz -i2 seq/{1}_2.fq.gz \
  -o temp/qc -v -t 3 --remove-intermediate-output \
  --trimmomatic ${soft}/envs/kneaddata/share/trimmomatic/ \
  --trimmomatic-options 'ILLUMINACLIP:~/miniconda3/envs/kneaddata/share/trimmomatic/adapters/TruSeq2-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50' \
  --reorder --bowtie2-options '--very-sensitive --dovetail' \
  -db ${db}/kneaddata/human_genome/hg37dec_v0.1"
代码解释
tail -n+2 result/metadata.txt ：这部分命令从名为“metadata.txt”的文件中提取数据，并丢弃第一行
cut -f1 ：cut命令用于从数据中选择第一个字段
Kneaddata -i1 seq/... -i2 seq/...  指定输入数据
    -o temp/qc   指定输出目录
    -v    用于启用详细的输出，以便查看处理进度和日志信息
    -t 3    指定线程数
    --remove-intermediate-output    这个选项将删除中间产生的临时文件，以节省存储空间
    --trimmomatic ${soft}/envs ...     指定trimmomatic质量控制工具的路径
    --trimmomatic-options     这个选项制定了trimmomatic的参数，包括：
        适配器裁剪ILLUMINACLIP
        滑动窗口裁剪SLIDINGWINDOW
        最小长度MINLEN
    --reorder    将对测序数据进行重新排序
    --bowtie2--options    指定了bowtie2比对工具的参数，包括
        非常敏感的比对模式 --very-sensitive
        交叉比对模式--dovetail
    -db ${db}/kneaddata/human...    指定去宿主使用的数据库
4.6、质量控制结果整理
# 先查看当前temp/qc目录大小
ls -lsh temp/qc

image.png

# 大文件清理，高宿主含量样本可节约>90%空间
/bin/rm -rf temp/qc/*contam* temp/qc/*unmatched*  temp/qc/*.fq
ls -lsh temp/qc/

image.png

# 列出temp/qc中所有文件或目录
ls -ltrsh temp/qc

image.png

# awk的system命令批处理改名，与fastp结果统一
awk '{system("/bin/mv -f temp/qc/"$1"_1_kneaddata_paired_1.fastq temp/qc/"$1"_1.fastq")}' <(tail -n+2 result/metadata.txt)
awk '{system("/bin/mv -f temp/qc/"$1"_1_kneaddata_paired_2.fastq temp/qc/"$1"_2.fastq")}' <(tail -n+2 result/metadata.txt)

image.png

4.7、质量控制结果汇总
# 创建结果汇总目录
mkdir -p result/qc
# 生成所有fastq文件摘要报告
kneaddata_read_count_table --input temp/qc \
    --output temp/kneaddata.txt
# 筛选重点结果列
cut -f 1,2,4,12,13 temp/kneaddata.txt | sed 's/_1_kneaddata//' > result/qc/sum.txt
cat result/qc/sum.txt
# 对齐方式查看表格
csvtk -t pretty result/qc/sum.txt

image.png

4.8、R语言统计质量控制结果
# 使用R语言生成统计摘要信息
Rscript -e "data=read.table('result/qc/sum.txt', header=T, row.names=1, sep='\t'); summary(data)"

image.png

# R转换宽表格为长表格
Rscript -e "library(reshape2); data=read.table('result/qc/sum.txt', header=T,row.names=1, sep='\t'); write.table(melt(data), file='result/qc/sum_long.txt',sep='\t', quote=F, col.names=T, row.names=F)"
cat result/qc/sum_long.txt

image.png

警告信息
# No id variables; using all as measure variables
提示该警告是因为函数未找到表示变量（id variables），于是将所有变量当做测量变量（measure variables）处理。

5、质量控制后质量评估

5.1、整理bowtie2, trimmomatic, fastqc报告
fastqc temp/qc/*.fastq -t 2
multiqc -d temp/qc/ -o result/qc/

image.png

结果查看：result/qc/multiqc_report_1.html

二、基于读长分析

1、Humann3软件介绍

Humann3是一种根据宏基因组或宏转录组数据高效准确的分析微生物代谢途径和其他分子功能丰度的方法；
使用与任何类型的微生物群落，不仅仅是人类微生物；

2、准备Humann3输入文件

Humann3要求双端序列合并的文件作为输入！
# 创建合并文件工作目录
mkdir -p temp/concat
# 双端合并为单个文件
for i in `tail -n+2 result/metadata.txt|cut -f1`;do 
  cat temp/qc/${i}_?.fastq \
  > temp/concat/${i}.fq; done
# 使用tree 查看目录结构
tree ../eaay_metagenome/temp/concat
# 查看样品数量和大小
ls -lsh temp/concat/*.fq

image.png

image.png

3、Humann3计算物种和功能组成

3.1、结构说明
输入文件：
    temp/concat/*.fq    每个样品质量控制后双端合并后的fastq序列
输出文件：
    temp/humann3/ 目录下
        C1_pathabundance.tsv
        C1_pathcoverage.tsv
        C1_genefamilies.tsv
整合后的输出：
    result/metaphlan3/taxonomy.tsv    物种丰度表
    result/metaphlan3/taxonomy.spf    物种丰度表（用于stamp分析）
    result/humann3/pathabundance_relab_unstratified.tsv    通路丰度表
    result/humann3/pathabundance_relab_stratified.tsv     通路物种组成丰度表
    stratified   每个菌对此功能通路组成的贡献
    unstratified    功能组成
3.2、humann3环境
# 激活conda虚拟环境humann3
conda activate humann3
# 设置metaphlan4数据库位置
DEFAULT_DB_FOLDER=${db}/metaphlan4
# 查看版本
humann3 --version    # humann v3.8
# 创建humann3工作目录
mkdir -p temp/humann3
3.3、单样本运行测试
# 单样本1.25M PE150运行测试
i=C1
time humann --input temp/concat/${i}.fq --output temp/humann3 --threads 1
# 比较慢，通过调整threads的值加快处理速度（取多大取决于CPU性能）

# 运行进程
Running metaphlan .........
Creating custom ChocoPhlAn database .......
Running bowtie2-build ........
Running bowtie2 ........
Running diamond ........
Computing gene families ...
Computing pathways abundance and coverage ...
运行结果
# 运行结果

image.png

# 目录结构
tree temp/humann3

image.png

3.4、多样本并行计算
# 多样本并行计算
tail -n+2 result/metadata.txt | cut -f1 | rush -j 2 \
  'humann --input temp/concat/{1}.fq  \
  --output temp/humann3/ --threads 18'

# 运行进程
Running metaphlan ........
Creating custom ChocoPhlAn database ........
Running bowtie2-build ........
Running bowtie2 ........
Running diamond ........
Computing gene families ...
Computing pathways abundance and coverage ...
Running metaphlan ........
Creating custom ChocoPhlAn database ........
Running bowtie2-build ........
Running bowtie2 ........
Running diamond ........
Computing gene families ...
Computing pathways abundance and coverage ...
运行结果
运行结果

image.png

image.png

# 查看工作目录运行结果
tree temp/humann3

image.png

3.5、链接重要文件到Humann3目录
# 链接重要文件至humann2目录
for i in $(tail -n+2 result/metadata.txt | cut -f1); do 
   ln -f temp/humann3/${i}_humann_temp/${i}_metaphlan_bugs_list.tsv temp/humann3/
done    
# 删除临时文件，极占用空间
/bin/rm -rf temp/concat/* temp/humann3/*_humann3_temp

4、物种组成表

4.1、样品结果合并
mkdir -p result/metaphlan4
# 合并、修正样本名、预览
merge_metaphlan_tables.py temp/humann3/*_metaphlan_bugs_list.tsv | \
  sed 's/_metaphlan_bugs_list//g' > result/metaphlan4/taxonomy.tsv
head -n5 result/metaphlan4/taxonomy.tsv
# 查看结果
tree result

image.png

4.2、转换为stamp的spf格式
# metaphlan4较2增加更多unclassifed和重复结果，用sort和uniq去除
metaphlan_to_stamp.pl result/metaphlan4/taxonomy.tsv \
  | sort -r | uniq > result/metaphlan4/taxonomy.spf


# 检查执行结果
tree result

image.png

# 查看taxonomy.spf文件内容
head result/metaphlan4/taxonomy.spf

image.png

5、功能组成分析

5.1、合并通路丰度
# 合并通路丰度(pathabundance)
mkdir -p result/humann3
humann_join_tables --input temp/humann3 \
  --file_name pathabundance \
  --output result/humann3/pathabundance.tsv
# 样本名调整：删除列名多余信息
head result/humann3/pathabundance.tsv
sed -i 's/_Abundance//g' result/humann3/pathabundance.tsv
# 查看pathabundance.tsv文件生成结果
tree result

image.png

# 预览和统计
head -n 3 result/humann3/pathabundance.tsv
csvtk -t stat result/humann3/pathabundance.tsv

image.png

5.2、标准化为相对丰度
humann_renorm_table \
  --input result/humann3/pathabundance.tsv \
  --units relab \
  --output result/humann3/pathabundance_relab.tsv
head -n5 result/humann3/pathabundance_relab.tsv

image.png

5.3、分层结果
# 功能和对应物种表(stratified)和功能组成表(unstratified)
humann_split_stratified_table \
  --input result/humann3/pathabundance_relab.tsv \
  --output result/humann3/
# 查看处理后文件生成结果
tree result/humann3

image.png

5.4、差异比较与柱状图
两样本无法组间比较，在pcl层面替换为HMP数据进行统计和可视化
流程概况：
    输入数据：通路丰度表格 result/humann3/pathabundance.tsv
    输入数据：实验设计信息 result/metadata.txt
    中间数据：包含分组信息的通路丰度表格文件 result/humann3/pathabundance.pcl
    输出结果：result/humann3/associate.txt
在通路丰度中添加分组
# 提取样品列表
head -n1 result/humann3/pathabundance.tsv | sed 's/# Pathway/SampleID/' | tr '\t' '\n' > temp/header
# 合成样本、分组+数据
cat <(head -n1 result/humann3/pathabundance.tsv) temp/group <(tail -n+2 result/humann3/pathabundance.tsv) \
  > result/humann3/pathabundance.pcl
head -n5 result/humann3/pathabundance.pcl
tail -n5 result/humann3/pathabundance.pcl
# 查看文件生成结果
tree result

image.png

组间比较，样本量少无差异，无法统计，使用hmp_pathabund.pcl统计和绘图
# 下载数据：
wget -c http://www.imeta.science/github/EasyMetagenome/result/humann2/hmp_pathabund.pcl
# 存放位置为result/humann3
# 查看文件位置：
tree result

image.png

# 统计表格行、列数量
pcl=result/humann3/hmp_pathabund.pcl
csvtk -t stat ${pcl}
head -n3 ${pcl} | cut -f 1-5
# 使用humann_associate进行关联分析
humann_associate --input ${pcl} \
    --focal-metadatum Group --focal-type categorical \
    --last-metadatum Group --fdr 0.05 \
    --output result/humann3/associate.txt
wc -l result/humann3/associate.txt
head -n5 result/humann3/associate.txt
barplot展示通路的物种组成
# barplot展示通路的物种组成，如：腺苷核苷酸合成
# --sort sum metadata 按丰度和分组排序
# 指定差异通路，如 P163-PWY / PWY-3781 / PWY66-409 / PWY1F-823
path=P163-PWY
humann_barplot --sort sum metadata \
    --input ${pcl} --focal-feature ${path} \
    --focal-metadata Group --last-metadata Group \
    --output result/humann3/barplot_${path}.pdf

# 检验文件生成
tree result

image.png

# 查看文件内容

image.png

6、转换为KEGG注释

6.1、重组基因家族文件为uniref90_ko文件
# 为什么重组
将代谢功能数据与代谢通路的功能注释相关联，从而提供更具生物学意义的数据分析。
# 什么是uniref90_ko
uniref90是一个用于减少蛋白质数据冗余的高分辨率系统；
uniref90_ko是uniref90的一个变种，具有相同的高分辨率，但是增加了与代谢通路功能的关联信息。

# 基因家族转换为uniref90_ko
for i in `tail -n+2 result/metadata.txt|cut -f1`; do
  humann_regroup_table \
    -i temp/humann3/${i}_genefamilies.tsv \
    -g uniref90_ko \
    -o temp/humann3/${i}_ko.tsv
done

# 检查文件生成结果
tree temp

image.png

6.2、合并uniref_ko文件
humann_join_tables \
  --input temp/humann3/ \
  --file_name ko \
  --output result/humann3/ko.tsv

# 检查文件生成结果
tree result

image.png

6.3、修改标题名称
sed -i '1s/_Abundance-RPKs//g' result/humann3/ko.tsv
tail result/humann3/ko.tsv
6.4、分层
# 对ko合并文件进行分层
humann_split_stratified_table \
  --input result/humann3/ko.tsv \
  --output result/humann3/ 

# 检查文件生成结果

image.png

6.5、summarizeAbundance生成COG/KO/CAZy丰度汇总表
# 使用summarizeAbundance生成KO丰度汇总表
summarizeAbundance.py \
  -i result/humann3/ko_unstratified.tsv \
  -m ${db}/EasyMicrobiome/kegg/KO1-4.txt \
  -c 2,3,4 -s ',+,+,' -n raw \
  -o result/humann3/KEGG

# 显示文件生成结果
tree result

image.png

7、GraPhlAn图

7.1、使用export2graphlan可视化物种丰度信息
export2graphlan.py --skip_rows 1,2 -i result/metaphlan4/taxonomy.tsv \
  --tree temp/merged_abundance.tree.txt \
  --annotation temp/merged_abundance.annot.txt \
  --most_abundant 1000 --abundance_threshold 20 --least_biomarkers 10 \
  --annotations 3,4 --external_annotations 7
7.2、使用graphlan_annotate.py生成注释信息
graphlan_annotate.py --annot temp/merged_abundance.annot.txt \
   temp/merged_abundance.tree.txt  temp/merged_abundance.xml
7.3、使用graphlan.py生成图像pdf
graphlan.py temp/merged_abundance.xml result/metaphlan4/graphlan.pdf \
   --external_legends 
# 查看文件生成结果
tree result

image.png

7.4、查看文件内容
查看graphlan.pdf

image.png

查看graphlan_annot.pdf

image.png

查看graphlan_legend.pdf

image.png

8、LEfSe差异分析物种

8.1、流程概况
输入文件：物种丰度表result/metaphlan4/taxonomy.tsv
输入文件：样品分组信息 result/metadata.txt
中间文件：整合后用于LefSe分析的文件 result/metaphlan4/lefse.txt，这个文件可以提供给www\.ehbio.com/ImageGP 用于在线LefSE分析
LefSe结果输出：result/metaphlan4/目录下lefse开头和feature开头的文件

前面演示数据仅有2个样本，无法进行差异比较。下面使用result12目录中由12个样本生成的结果表进行演示
8.2、准备输入数据
# 创建result/metaphlan2目录
mkdir -p result/metaphlan2
# 使用EasyMetagenome中的数据
cp ~/database/EasyMetagenome/result12/metaphlan2/taxonomy.tsv ~/workspace/easy_metagenome/result/metaphlan2/
# 查看当前目录结构（检查输入数据）
tree result

image.png

# 预览数据
head -n3 result/metaphlan2/taxonomy.tsv

image.png

# 提取样本行，替换为每个样本一行，修改ID为SampleID
head -n1 result/metaphlan2/taxonomy.tsv|tr '\t' '\n'|sed '1 s/ID/SampleID/' >temp/sampleid
head -n3 temp/sampleid

image.png

# 提取SampleID对应的分组Group(假设为metadata.txt中第二列$2)，替换换行\n为制表符\t，再把行末制表符\t替换回换行
awk 'BEGIN{OFS=FS="\t"}NR==FNR{a[$1]=$2}NR>FNR{print a[$1]}' result/metadata.txt temp/sampleid|tr '\n' '\t'|sed 's/\t$/\n/' >groupid
cat groupid

image.png

# 合并分组和数据(替换表头)
cat groupid <(tail -n+2 result/metaphlan2/taxonomy.tsv) > result/metaphlan2/lefse.txt
head -n3 result/metaphlan2/lefse.txt

image.png

8.3 使用ImageGP一键绘制
网址：https://www.bic.ac.cn/ImageGP/
输入：lefse.txt文本内容
绘制结果：

image.png

8.4、LEfSe分析
# 激活lefse环境
conda activate lefse
# 格式转换为lefse内部格式
lefse-format_input.py result/metaphlan2/lefse.txt \
  temp/input.in -c 1 -o 1000000
# 运行lefse(样本无重复、分组将报错)
run_lefse.py temp/input.in temp/input.res
8.5、绘制物种树注释差异
lefse-plot_cladogram.py temp/input.res \
   result/metaphlan2/lefse_cladogram.pdf --format pdf
# 查看绘制文件生成结果
tree result

image.png

# 查看绘制结果：lefse_cladogram.pdf

image.png

8.5、绘制所有差异features柱状图
lefse-plot_res.py temp/input.res \
  result/metaphlan2/lefse_res.pdf --format pdf
# 查看生成结果

image.png image.png

8.6、绘制单个features柱状图
# 查看显著差异features，按丰度排序
grep -v '-' temp/input.res | sort -k3,3n 
# 手动选择指定feature绘图，如Firmicutes
lefse-plot_features.py -f one --format pdf \
  --feature_name "k__Bacteria.p__Firmicutes" \
  temp/input.in temp/input.res \
  result/metaphlan2/lefse_Firmicutes.pdf

image.png image.png

8.7、批量绘制所有差异features柱状图
lefse-plot_features.py -f diff \
  --archive none --format pdf \
  temp/input.in temp/input.res \
  result/metaphlan2/lefse_

查看生成文件结果：
查看生生文件内容（省略）

image.png

9、Kraken2+Bracken物种注释和丰度估计

9.1、简介与环境
Kraken2可以快速完成读长(read)层面的物种注释和定量；
还可以进行重叠群(contig)、基因(gene)、宏基因组组装基因组(MAG/bin)层面的序列物种注释。

# 启动kraken2环境
conda activate kraken2
# 查看核心版本：
kraken2 --version  # Kraken version 2.1.2
9.2、Kraken2物种注释流程概况
输入：temp/qc/{1}\_?.fastq 质控后的数据，{1}代表样本名
参考数据库：-db ~、database/kraken2/pluspfp/，pluspfp代表最全库
输出结果：每个样本单独输出，temp/kraken2/{1}\_report和temp/kraken2/{1}\_output
整合输出结果：result/kraken2/taxonomy\_count.txt 物种丰度表
9.3、Kraken2物种注释
# 创建中间文件存储目录
mkdir -p temp/kraken2
单样本注释
i=C1
time kraken2 --db ~/database/kraken2/pluspf/ --paired temp/qc/${i}_?.fastq \
  --threads 18 --use-names --report-zero-counts \
  --report temp/kraken2/${i}.report \
  --output temp/kraken2/${i}.output

# 查看文件生成结果
tree temp

image.png

多样本并行生成report
tail -n+2 result/metadata.txt | cut -f1 | rush -j 2 \
  "kraken2 --db ~/database/kraken2/pluspf/ --paired temp/qc/{1}_?.fastq \
  --threads 18 --use-names --report-zero-counts \
  --report temp/kraken2/{1}.report \
  --output temp/kraken2/{1}.output"

image.png

9.4、使用krakentools转换report为mpa格式
for i in `tail -n+2 result/metadata.txt | cut -f1`;do
  kreport2mpa.py -r temp/kraken2/${i}.report \
  --display-header -o temp/kraken2/${i}.mpa; done

image.png

9.5、合并样本为表格
# 创建输入结果存储目录
mkdir -p result/kraken2
# 输出结果行数相同，但不一定顺序一致，要重新排序
tail -n+2 result/metadata.txt|cut -f1|rush -j 1 \
  'tail -n+2 temp/kraken2/{1}.mpa | LC_ALL=C sort | cut -f 2 | sed "1 s/^/{1}\n/" > temp/kraken2/{1}_count '
# 提取第一样本品行名为表行名
header=`tail -n 1 result/metadata.txt | cut -f 1`
echo $header
tail -n+2 temp/kraken2/${header}.mpa | LC_ALL=C sort | cut -f 1 | \
  sed "1 s/^/Taxonomy\n/" > temp/kraken2/0header_count
head -n3 temp/kraken2/0header_count
# paste合并样本为表格
ls temp/kraken2/*count
paste temp/kraken2/*count > result/kraken2/tax_count.mpa
# 检查表格及统计
csvtk -t stat result/kraken2/tax_count.mpa
head -n 5 result/kraken2/tax_count.mpa

image.png

10、Bracken丰度估计

10.1、参数说明
-d为数据库，-i为输入kraken2报告文件
r是读长，此处为100，通常为150，o输出重新估计的值
l为分类级，可选域D、门P、纲C、目O、科F、属G、种S级别丰度估计
t是阈值，默认为0，越大越可靠，但可用数据越少
10.2、运行Bracken进行风度分析
# 输入文件为kraken2生成report
# 设置分类级别
tax=P   # P 门
# 创建中间文件存储目录
mkdir -p temp/bracken
i=C1

for i in `tail -n+2 result/metadata.txt | cut -f1`;do
  bracken -d ~/database/kraken2/pluspf/ \
    -i temp/kraken2/${i}.report \
    -r 100 -l ${tax} -t 0 \
    -o temp/bracken/${i}; done
    
# 查看文件生成结果
tree temp

image.png

10.3、bracken结果合并成表
# 输出结果行数相同，但不一定顺序一致，要按表头排序
# 仅提取第6列reads count，并添加样本名
tail -n+2 result/metadata.txt | cut -f1 | rush -j 1 \
  'tail -n+2 temp/bracken/{1} | LC_ALL=C sort | cut -f6 | sed "1    s/^/{1}\n/" > temp/bracken/{1}.count'
# 提取第一样本品行名为表行名
h=`tail -n1 result/metadata.txt|cut -f1`
tail -n+2 temp/bracken/${h} | sort | cut -f1 | \
  sed "1 s/^/Taxonomy\n/" > temp/bracken/0header.count
# 检查文件数，为n+1
ls temp/bracken/*count | wc
# paste合并样本为表格，并删除非零行
paste temp/bracken/*count > result/kraken2/bracken.${tax}.txt
# 统计行列，默认去除表头
csvtk -t stat result/kraken2/bracken.${tax}.txt
# 按频率过滤，-r可标准化，-e过滤(microbiome_helper)
Rscript ~/database/EasyMicrobiome/script/filter_feature_table.R \
  -i result/kraken2/bracken.${tax}.txt \
  -p 0.01 \
  -o result/kraken2/bracken.${tax}.0.01
csvtk -t stat result/kraken2/bracken.${tax}.0.01
10.4、个性化结果筛选
# 门水平去除脊索动物(人)
grep 'Chordata' result/kraken2/bracken.P.0.01
grep -v 'Chordata' result/kraken2/bracken.P.0.01 >  result/kraken2/bracken.P.0.01-H


# 上面没有进行种水平的处理，这里没有bracken.S.0.01文件

# 按物种名手动去除宿主污染，以人为例(需按种水平计算相关结果)
# 种水平去除人类P:Chordata,S:Homo sapiens
grep 'Homo sapiens' result/kraken2/bracken.S.0.01
grep -v 'Homo sapiens' result/kraken2/bracken.S.0.01 > result/kraken2/bracken.S.0.01-H

# 分析后清理每条序列的注释大文件
/bin/rm -rf temp/kraken2/*.output

三、组装分析流程

1、MegaHit组装

1.1、megahit组装环境
# 启动工作环境
conda activate megahit

megahit -v   # MEGAHIT v1.2.9
metaspades.py -v   # SPAdes genome assembler v3.15.4 [metaSPAdes mode]
quast.py -v   # QUAST v5.0.2
prodigal -v   #Prodigal V2.6.3: February, 2016
1.2、megahit组装
megahit -t 3 \
  -1 `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/qc\//;s/$/_1.fastq/'|tr '\n' ','|sed 's/,$//'` \
  -2 `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/qc\//;s/$/_2.fastq/'|tr '\n' ','|sed 's/,$//'` \
  -o temp/megahit 

# 统计大小
seqkit stat temp/megahit/final.contigs.fa

image.png

# 预览重叠群最前6行，前60列字符
head -n6 temp/megahit/final.contigs.fa | cut -c1-60

image.png

# 备份重要结果
mkdir -p result/megahit/
ln -f temp/megahit/final.contigs.fa result/megahit/
# 删除临时文件
rm -rf temp/megahit/intermediate_contigs

2、metaSPAdes拼装

2.1、metaSPAdes特点
metaSPAdes提供更准确的拼装结果，同时要占用更多的内存资源，花费更长的处理时间；
metaSPAdes支持二、三代混合组装。
2.2、metaSPAdes组装
# 精细但使用内存和时间更多
/usr/bin/time -v -o metaspades.py.log metaspades.py -t 3 -m 100 \
  `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/qc\//;s/$/_1.fastq/'|sed 's/^/-1 /'| tr '\n' ' '` \
  `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/qc\//;s/$/_2.fastq/'|sed 's/^/-2 /'| tr '\n' ' '` \
  -o temp/metaspades
2.3、组装结果
# 查看软件时间User time和内存Maximum resident set size
cat metaspades.py.log

image.png

# 组装的contigs体积更大
ls -sh temp/metaspades/contigs.fasta
seqkit stat temp/metaspades/contigs.fasta

image.png

# 备份重要结果
mkdir -p result/metaspades/
ln -f temp/metaspades/contigs.fasta result/metaspades/
# 删除临时文件
/bin/rm -rf temp/metaspades

3、QUAST评估

3.1、QUAST评估
QUAST评估可以生成report文本tsv/txt、网页html、PDF等格式报告
3.2、对megahit组装结果进行评估
quast.py result/megahit/final.contigs.fa \
  -o result/megahit/quast -t 2

# 查看评估结果生成文件
tree result

image.png

3.3、megahit和metaSPAdes进行比较
quast.py --label "megahit,metapasdes" \
    result/megahit/final.contigs.fa \
    result/metaspades/contigs.fasta \
    -o result/quast

# 查看文件生成结果
tree result

image.png

4、metaquast评估

4.1、MetaQUAST和QUAST
QUAST计算和报告了各种组装质量统计指标，如N50、L50、GC含量等，以帮助用户评估组装的质量。
MetaQUAST与QUAST类似，但它可以处理混合物中的多个微生物基因组的组装；
MetaQUAST还可以生成用于评估混合物中各种微生物的组装质量的统计信息。
MetaQUAST评估，更全面，但需下载相关数据库，速度会慢很多。
4.2、MetaQUAST评估
# 使用MetaQUAST对megahit组装结果进行评估
time metaquast.py result/megahit/final.contigs.fa \
  -o result/megahit/metaquast

5、基因预测、去冗余和定量

5.1、MetaProdigal基因预测概况
# 输入文件：组装的序列 result/megahit/final.contigs.fa
# 输出文件：prodigal预测的基因序列 temp/prodigal/gene.fa
5.2、MetaProdigal基因预测
# 创建中间文件生成目录
mkdir -p temp/prodigal
# 执行prodigal基因预测
prodigal -i result/megahit/final.contigs.fa \
    -d temp/prodigal/gene.fa \
    -o temp/prodigal/gene.gff \
    -p meta -f gff > temp/prodigal/gene.log 2>&
# 检查生成日志文件，查看执行是否发生错误
tail temp/prodigal/gene.log
# 统计基因数量
seqkit stat temp/prodigal/gene.fa 
# 统计完整基因数量，数据量大可只用完整基因部分
grep -c 'partial=00' temp/prodigal/gene.fa 
# 提取完整基因(完整片段获得的基因全为完整，如成环的细菌基因组)
grep 'partial=00' temp/prodigal/gene.fa | cut -f1 -d ' '| sed 's/>//' > temp/prodigal/full_length.id

seqkit grep -f temp/prodigal/full_length.id temp/prodigal/gene.fa > temp/prodigal/full_length.fa

seqkit stat temp/prodigal/full_length.fa

image.png

6、cd-hit基因聚类/去冗余

6.1、cd-hit处理概况
# 输入文件：prodigal预测的基因序列 temp/prodigal/gene.fa
# 输出文件：去冗余后的基因和蛋白序列：result/NR/nucleotide.fa, result/NR/protein.fa
6.2、cd-hit基因预测
# 创建处理结果文件生成目录
mkdir -p result/NR
# aS覆盖度，c相似度，G局部比对，g最优解，T多线程，M内存0不限制
# 2万基因2m，2千万需要2000h，多线程可加速
cd-hit-est -i temp/prodigal/gene.fa \
    -o result/NR/nucleotide.fa \
    -aS 0.9 -c 0.95 -G 0 -g 0 -T 0 -M 0
# 统计非冗余基因数量，单次拼接结果数量下降不大，多批拼接冗余度高
grep -c '>' result/NR/nucleotide.fa
# 翻译核酸为对应蛋白序列, --trim去除结尾的*
seqkit translate --trim result/NR/nucleotide.fa \
    > result/NR/protein.fa

# 两批数据去冗余使用cd-hit-est-2d加速，见附录

7、salmon基因定量

7.1、Salmon处理概况
# 输入文件：去冗余后的基因序列：result/NR/nucleotide.fa
# 输出文件：Salmon定量：result/salmon/gene.count, gene.TPM
7.2、Salmon定量
# 创建Salmon定量处理中间文件生成目录
mkdir -p temp/salmon
# 建索引, -t序列, -i 索引，10s
salmon index -t result/NR/nucleotide.fa \
  -p 3 -i temp/salmon/index

# 定量，l文库类型自动选择，p线程，--meta宏基因组模式, 2个任务并行2个样
tail -n+2 result/metadata.txt | cut -f1 | rush -j 2 \
  "salmon quant -i temp/salmon/index -l A -p 3 --meta \
    -1 temp/qc/{1}_1.fastq -2 temp/qc/{1}_2.fastq \
    -o temp/salmon/{1}.quant"
7.3、结果整合
# 合并，创建整合结果存储目录
mkdir -p result/salmon
# 使用salmon quantmerge合并基因定量结果
salmon quantmerge --quants temp/salmon/*.quant \
    -o result/salmon/gene.TPM
# 在合并基因定量结果是指令合并列
salmon quantmerge --quants temp/salmon/*.quant \
    --column NumReads -o result/salmon/gene.count
# 修改文件的标题行或文件名
sed -i '1 s/.quant//g' result/salmon/gene.*
# 预览结果表格
head -n3 result/salmon/gene.*

8、功能基因注释

8.1、功能基因注释概况
# 输入数据：上一步预测的蛋白序列 result/NR/protein.fa
# 中间结果：temp/eggnog/protein.emapper.seed_orthologs
#    temp/eggnog/output.emapper.annotations
#    temp/eggnog/output

# COG定量表：result/eggnog/cogtab.count
#    result/eggnog/cogtab.count.spf (用于STAMP)

# KO定量表：result/eggnog/kotab.count
#    result/eggnog/kotab.count.spf  (用于STAMP)

# CAZy碳水化合物注释和定量：result/dbcan2/cazytab.count
#    result/dbcan2/cazytab.count.spf (用于STAMP)

# 抗生素抗性：result/resfam/resfam.count
#    result/resfam/resfam.count.spf (用于STAMP)

# 这部分可以拓展到其它数据库
8.2、eggNOG基因注释
# 激活环境与版本检查
# 使用的是emapper功能注释和分析蛋白质序列工具
conda activate eggnog
emapper.py --version   # emapper-2.1.6

# 创建工作目录
mkdir -p temp/eggnog
mkdir -p result/eggnog

# 运行emapper
time emapper.py --data_dir ${db}/eggnog \
  -i result/NR/protein.fa --cpu 3 -m diamond --override \
  -o temp/eggnog/output
8.3、运行结果处理
# 查看注释生成中间文件结果
tree temp/eggnog

image.png

# 格式化结果并显示表头
grep -v '^##' temp/eggnog/output.emapper.annotations | sed '1 s/^#//' \
  > temp/eggnog/output
csvtk -t headers -v temp/eggnog/output
8.4、丰度汇总表
# 使用summarizeAbundance生成COG/KO/CAZy丰度汇总表
# 7列COG_category按字母分隔，12列KEGG_ko和19列CAZy
summarizeAbundance.py \
  -i result/salmon/gene.TPM \
  -m temp/eggnog/output \
  -c '7,12,19' -s '*+,+,' -n raw \
  -o result/eggnog/eggnog
  
sed -i 's#^ko:##' result/eggnog/eggnog.KEGG_ko.raw.txt
sed -i '/^-/d' result/eggnog/eggnog*

# 添加注释生成STAMP的spf格式
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
  ${db}/EasyMicrobiome/kegg/KO_description.txt \
  result/eggnog/eggnog.KEGG_ko.raw.txt | \
  sed 's/^\t/Unannotated\t/' \
  > result/eggnog/eggnog.KEGG_ko.TPM.spf
head -n 5 result/eggnog/eggnog.KEGG_ko.TPM.spf

# KO to level 1/2/3
summarizeAbundance.py \
  -i result/eggnog/eggnog.KEGG_ko.raw.txt \
  -m ${db}/EasyMicrobiome/kegg/KO1-4.txt \
  -c 2,3,4 -s ',+,+,' -n raw \
  -o result/eggnog/KEGG

image.png

# CAZy
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
   ${db}/EasyMicrobiome/dbcan2/CAZy_description.txt result/eggnog/eggnog.CAZy.raw.txt | \
  sed 's/^\t/Unannotated\t/' > result/eggnog/eggnog.CAZy.TPM.spf

# COG
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2"\t"$3} NR>FNR{print a[$1],$0}' \
  ${db}/EasyMicrobiome/eggnog/COG.anno result/eggnog/eggnog.COG_category.raw.txt > \
  result/eggnog/eggnog.COG_category.TPM.spf

image.png

9、CAZy碳水化合物酶

9.1、注释碳水化合物酶
# 创建中间文件存储目录
mkdir -p temp/dbcan2
# 执行diamond blastp比对
diamond blastp \
  --db ${db}/dbcan2/CAZyDB.08062022 \
  --query result/NR/protein.fa \
  --threads 18 -e 1e-3 --outfmt 6 --max-target-seqs 1 --quiet \
  --out temp/dbcan2/gene_diamond.f6

# 统计文本行数
wc -l temp/dbcan2/gene_diamond.f6
9.2、结果处理
# 创建输出结果存储目录
mkdir -p result/dbcan2

# 提取基因与dbcan分类对应表
format_dbcan2list.pl \
  -i temp/dbcan2/gene_diamond.f6 \
  -o temp/dbcan2/gene.list

# 按对应表累计丰度，依赖
summarizeAbundance.py \
  -i result/salmon/gene.TPM \
  -m temp/dbcan2/gene.list \
  -c 2 -s ',' -n raw \
  -o result/dbcan2/TPM

# 添加注释生成STAMP的spf格式，结合metadata.txt进行差异比较
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
   ${db}/EasyMicrobiome/dbcan2/CAZy_description.txt result/dbcan2/TPM.CAZy.raw.txt | \
  sed 's/^\t/Unannotated\t/' \
  > result/dbcan2/TPM.CAZy.raw.spf

head result/dbcan2/TPM.CAZy.raw.spf
# 检查未注释数量，有则需要检查原因
grep 'Unannotated' result/dbcan2/TPM.CAZy.raw.spf|wc -l

10、CARD耐药基因

10.1、CARD说明
CARD在线分析平台：https://card.mcmaster.ca/ 
本地软件使用教程: https://github.com/arpcard/rgi
参考文献：http://doi.org/10.1093/nar/gkz935
10.2、环境准备
# 启动RGI环境
conda activate rgi
# 查看核心版本
rgi main -v   # 5.2.0
10.3、简化蛋白ID
# 创建结果生成文件存储目录
mkdir result/card
# 简化蛋白ID
cut -f 1 -d ' ' result/NR/protein.fa > temp/protein.fa

grep '>' result/NR/protein.fa | head -n 3
grep '>' temp/protein.fa | head -n 3
10.4、蛋白层面注释ARG
rgi main -i temp/protein.fa -t protein \
  -n 9 -a DIAMOND --include_loose --clean \
  -o result/card/protein
head -n3 result/card/protein.txt
10.5、基因层面注释ARG
cut -f 1 -d ' ' result/NR/nucleotide.fa > temp/nucleotide.fa
grep '>' temp/nucleotide.fa | head -n3
rgi main -i temp/nucleotide.fa -t contig \
  -n 9 -a DIAMOND --include_loose --clean \
  -o result/card/nucleotide
head -n3 result/card/nucleotide.txt
10.6、重叠群层面注释ARG
cut -f 1 -d ' ' result/megahit/final.contigs.fa > temp/contigs.fa
grep '>' temp/contigs.fa | head -n3
rgi main -i temp/contigs.fa -t contig \
  -n 9 -a DIAMOND --include_loose --clean \
  -o result/card/contigs
head result/card/contigs.txt
10.7、结果说明
- protein.json，在线可视化
- protein.txt，注释基因列表

image.png

11、基因物种注释

11.1、使用Kraken2进行物种注释
conda activate kraken2
kraken2 --db ${db}/kraken2/pluspf \
  result/NR/nucleotide.fa \
  --threads 18 \
  --report temp/NRgene.report \
  --output temp/NRgene.output

# 对输出结果进行整理
grep '^C' temp/NRgene.output | cut -f 2,3 | sed '1 i Name\ttaxid' \
  > temp/NRgene.taxid

awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$0} NR>FNR{print $1,a[$2]}' \
  /db/EasyMicrobiome/kraken2/taxonomy.txt \
  temp/NRgene.taxid \
  > result/NR/nucleotide.tax

四、分箱挖掘单菌基因组

1、MetaWRAP混合样本分箱

1.1、MetaWrap介绍
主页：https://github.com/bxlab/metaWRAP
教程: https://github.com/bxlab/metaWRAP/blob/master/Usage_tutorial.md

挖掘单菌基因组，需要研究对象复杂度越低、测序深度越大，结果质量越好；
要求单样本6GB+，复杂样本如土壤推荐数据量30GB+，至少3个样本。

演示数据12个样仅140MB，无法获得单菌基因组，这里使用官方测序数据演示；
软件和数据库布置需1-3天，演示数据分析过程超10h，30G样也需1-30天，由服务器性能决定。
1.2、工作环境准备
# 创建工作环境目录、进入工作目录
mkdir -p meta/binning && cd meta/binning

# 初始化项目，创建详细工作目录结构
mkdir -p temp/hr seq result

# 启动metawrap环境
conda activate metawrap
1.3、概要说明
# 1.要求
这里基于质控clean数据和拼接好的重叠群contigs，基于上游结果继续分析；
如果上游测试数据过小，分箱无结果；
软件推荐至少7G数据。

# 2.输入输出文件
输入：
    - 质控后序列，文件名格式为*_1.fastq和*_2.fastq；
    - 组装的重叠群文件：megahit/final_configs.fa
输出：
    - Binning结果：temp/binning
    - 提纯后的Bin统计结果：temp/bin_refinement/metawrap_50_10_bins.stats
    - Bin定量结果文件和图：binning/temp/bin_quant/bin_abundance_table.tab 和 bin_abundance_heatmap.png
    - Bin物种注释：binning/temp/bin_classify/bin_taxonomy.tab
    - Prokka基因预测：binning/temp/bin_annotate/prokka_out/bin.*.ffn 核酸序列
    - Bin可视化图表：binning/temp/bloblogy/final.contigs.binned.blobplot (数据表) 和 blobplot_figures (可视化图)
1.3、准备输入文件
# 1.准备质控后的数据
cd temp/hr
for i in `seq 7 9`;do
    wget -c ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR01134${i}/ERR01134${i}_1.fastq.gz
    wget -c ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR01134${i}/ERR01134${i}_2.fastq.gz
done
# 解压数据文件、-k表示保留原文件
gunzip -k *.gz
# 批量修改扩展名fq为fastq
# rename .fq .fastq *.fq

image.png

# 2.megahit拼接结果
# 创建拼接工作目录
cd ~/workspace/easy_metagenome/meta/binning
mkdir -p temp/megahit
cd temp/megahit
# 下载megahit组装数据
wget -c http://www.imeta.science/db/metawrap/final.contigs.fa.gz
# 解压数据文件
gunzip -k *.gz

image.png

2、分箱Binning

2.1、使用MetaWrap进行分箱操作
# 进入工作目录并加载工作环境
cd ~/workspace/easy_metagenome/meta/binning
conda activate metawrap
# 检查核心版本
metawrap -v
# 检查当前工作目录结构
tree ../binning

image.png

# 分箱
nohup metawrap binning -o temp/binning \
  -t 36 -a temp/megahit/final.contigs.fa \
  --metabat2 --maxbin2 --concoct \
  temp/hr/*.fastq &

# 代码解读
# nohup command & 将命令或进程置于后台运行，关闭终端并不会中断进程
# -o  指定输出目录
# -t  指定线程数
# -a  指定用于分箱的输入序列文件
# --metabat2 --maxbin2 --concoct用于分箱的工具

# 查看进程运行状态
jobs
# 查询结果如下图所示，Running表示进程仍在运行中

image.png

# 获取更详细的进程运行信息
htop
# 查询结果如下图所示：
    所有CPU资源跑满；
    已占用CPU时间为7h

image.png

# 运行过程记录会保存在 nohup.out文件中
tail nohup.out

# 检查运行成功后、可以删除日志文件nohup.out
rm -f nohup.out

# 查看运行结果-文件生成结果
tree ../binning

image.png

2.2、分箱提纯
# 创建分箱提纯工作目录
mkdir -p temp/bin_refinement
# 分箱提纯
nohup metawrap bin_refinement \
  -o temp/bin_refinement \
  -A temp/binning/metabat2_bins/ \
  -B temp/binning/maxbin2_bins/ \
  -C temp/binning/concoct_bins/ \
  -c 50 -x 10 -t 36 &

# 代码解释
# -o  指定输出目录
# -A -B -C  指定分箱过程分箱工具的输出目录，分别对应于MetaBAT2、MaxBin2、CONCOCT
# -c 50  表示至少 50% 的相似性时才会被合并
# -x 10  表示最多进行10次重新拆分操作，目的是执行多次重新拆分以提高分箱的准确性和优化分箱结果

# 统计高质量Bin的数量，10个
wc -l temp/bin_refinement/metawrap_50_10_bins.stats
# 分析比较图见 temp/bin_refinement/figures/

image.png

intermediate_binning_results

intermediate_binning_results.png

binning_results

binning_results.png

2.3、分箱结果整理
# 将所有分箱结果保存到同一目录
mkdir -p temp/drep_in
ln -s ~/workspace/easy_metagenome/meta/binning/temp/bin_refinement/metawrap_50_10_bins/bin.* temp/drep_in/
ls -l temp/drep_in/
# 改名Ubuntu
rename s/bin./Mx_All_/ temp/drep_in/bin.*
ls -l temp/drep_in/

image.png

image.png

3、dRep去冗余种/株基因组集

3.1、工作环境准备
# 进入虚拟环境drep和工作目录
conda activate drep
cd ~/workspace/easy_metagenome/meta/binning
3.2、按种水平去冗余
mkdir -p temp/drep95
# rm -rf temp/drep95/data/checkM
dRep dereplicate temp/drep95/ \
  -g temp/drep_in/*.fa  \
  -sa 0.95 -nc 0.30 -comp 50 -con 10 -p 5

# drep去冗余处理代码解释
# drep dereplicate  去除冗余命令
# temp/drep95  指定结果文件和输出的目标目录
# -g  指定要进行去冗余的输入文件
# -sa 0.95  设置去冗余的相似性阈值，相似性为95%的被认为是重复基因组
# -nc 0.30  设置去冗余的基因组长度阈值
# -comp 50  设置去冗余的基因组的完整性阈值
# -con 10   设置去冗余的基因组的污染程度的阈值
# -p 10   设置线程数，即使用处理器数量
处理结束

image.png

主要结果-生成文件概况
非冗余基因组集：temp/drep95/dereplicated_genomes/\*.fa
聚类信息表：temp/drep95/data_tables/Cdb.csv
聚类和质量图：temp/drep95/figures/\*clustering\*

image.png

image.png

image.png

3.3、按株水平去冗余
mkdir -p temp/drep99
dRep dereplicate temp/drep99/ \
  -g temp/drep_in/*.fa \
  -sa 0.99 -nc 0.30 -comp 50 -con 10 -p 18
ls -l temp/drep99/dereplicated_genomes/ | grep '.fa' | wc -l
处理结束

image.png

4、CoverM基因组定量

# 基因组定量是宏基因组学分析的一个重要步骤，涉及到确定从宏基因组样本中检测到的
微生物基因组的数量和相对丰度。

# 基因组定量操作的目的是测量或估计给定样本中不同微生物基因组的相对丰度或存在程度。

# CoverM是一个用于进行基因组覆盖度计算的工具，可以帮助确定样本中不同微生物的
相对覆盖情况，从而提供关于这些微生物的相对丰度信息。
4.1、环境准备
# 启动环境
conda activate coverm
# 创建工作目录
mkdir -p temp/coverm
4.2、单样本测试
i=ERR011347
time coverm genome --coupled temp/hr/${i}_1.fastq temp/hr/${i}_2.fastq \
  --genome-fasta-directory temp/drep95/dereplicated_genomes/ -x fa \
  -o temp/coverm/${i}.txt
cat temp/coverm/${i}.txt

# 代码解析
# time  用于测量命令执行所需的时间
# coverm genome  运行coverm命令，执行基因组覆盖度计算
# --coupled  用于指定用于计算覆盖度的样本文件
# --genome-fasta-directory  这个选项指定了参考基因组的目录，用于计算覆盖度
# -x fa   这个选项指定参考基因组文件的扩展名为“.fa”，告诉coverm工具参考基因组文件的类型
# -o  指定输出文件的位置和名称
4.3、并行计算
tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \
  "coverm genome --coupled temp/hr/{}_1.fastq temp/hr/{}_2.fastq \
  --genome-fasta-directory temp/drep95/dereplicated_genomes/ -x fa \
  -o temp/coverm/{}.txt"
4.4、结果合并
mkdir -p result/coverm
conda activate humann3
sed -i 's/_1.fastq Relative Abundance (%)//' temp/coverm/*.txt
humann_join_tables --input temp/coverm \
  --file_name txt \
  --output result/coverm/abundance.tsv

5、GTDB-tk物种注释和进化树

5.1、工作环境准备
# 激活conda虚拟环境
conda activate gtdbtk
# 设置数据库路径
export GTDBTK_DATA_PATH="~/database/gtdb"
# 版本检查
gtdbtk -v
5.2、代表性细菌基因组物种注释
mkdir -p temp/gtdb_all
# 10000个基因组，32p，100min
gtdbtk classify_wf \
    --genome_dir temp/drep_in/ \
    --out_dir temp/gtdb_all \
    --extension fa --skip_ani_screen \
    --prefix tax \
    --cpus 36
5.3、代表种注释
mkdir -p temp/gtdb_classify
# 10个基因组，24p，100min 152G内存; 6p, 22基因组，1h
gtdbtk classify_wf \
    --genome_dir temp/drep95/dereplicated_genomes \
    --out_dir temp/gtdb_classify \
    --extension fa --skip_ani_screen \
    --prefix tax \
    --cpus 6
# less -S按行查看，按q退出
# less -S temp/gtdb_classify/tax.bac120.summary.tsv
head temp/gtdb_classify/tax.bac120.summary.tsv
head temp/gtdb_classify/tax.ar53.summary.tsv
5.4、多序列对齐结果建树
mkdir -p temp/gtdb_infer
gtdbtk infer --msa_file temp/gtdb_classify/align/tax.bac120.user_msa.fasta.gz \
    --out_dir temp/gtdb_infer --prefix tax --cpus 36
    
树文件`tax.unrooted.tree`可使用iTOL在线美化，也可使用GraphLan本地美化

制作树注释文件：以gtdb-tk物种注释(tax.bac120.summary.tsv)和drep基因组评估(Widb.csv)信息为注释信息
mkdir -p result/itol
# 制作分类学表
tail -n+2 temp/gtdb_classify/tax.bac120.summary.tsv|cut -f 1-2|sed 's/;/\t/g'|sed '1 s/^/ID\tDomain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies\n/' \
  > result/itol/tax.txt
head result/itol/tax.txt
# 基因组评估信息
sed 's/,/\t/g;s/.fa//' temp/drep95/data_tables/Widb.csv|cut -f 1-7,11|sed '1 s/genome/ID/' \
  > result/itol/genome.txt
head result/itol/genome.txt
# 整合注释文件
awk 'BEGIN{OFS=FS="\t"} NR==FNR{a[$1]=$0} NR>FNR{print $0,a[$1]}' result/itol/genome.txt result/itol/tax.txt|cut -f 1-8,10- > result/itol/annotation.txt
head result/itol/annotation.txt
# 添加各样本相对丰度(各组替换均值)
awk 'BEGIN{OFS=FS="\t"} NR==FNR{a[$1]=$0} NR>FNR{print $0,a[$1]}' <(sed '1 s/Genome/ID/' result/coverm/abundance.tsv) result/itol/annotation.txt|cut -f 1-15,17- > result/itol/annotation2.txt
head result/itol/annotation2.txt
1 Gtdbtk报错,重新安装
安装gtdbtk
conda create -y -n gtdbtk
conda install -c bioconda getdbtk
提示信息
GTDB-Tk v2.3.2 requires ~78G of external data which needs to be downloaded
and extracted. This can be done automatically, or manually.

Automatic:

    1. Run the command "download-db.sh" to automatically download and extract to:
        /home/users-data/liukai/anaconda3/envs/gtdbtk/share/gtdbtk-2.3.2/db/

Manual:

    1. Manually download the latest reference data:
        wget https://data.gtdb.ecogenomic.org/releases/release214/214.0/auxillary_files/gtdbtk_r214_data.tar.gz

    2. Extract the archive to a target directory:
        tar -xvzf gtdbtk_r214_data.tar.gz -C "/path/to/target/db" --strip 1 > /dev/null
        rm gtdbtk_r214_data.tar.gz

    3. Set the GTDBTK_DATA_PATH environment variable by running:
        conda env config vars set GTDBTK_DATA_PATH="/path/to/target/db"
解决办法
GTDBtk安装和使用：https://juejin.cn/post/7295613501342269459
avatar
文章
阅读
粉丝