宏基因组单样品vamb分箱，gtdb物种注释与建树

个人记录，其实很多文件夹的设置并不是那么合理。

先放一些分箱相关的来来来，一起来pick宏基因组binning分析工具 (baidu.com)

存放已有序列的文件夹

• ./qc 质控后的双端测序文件.fq.gz
• ./ 当前文件夹，放.fa后缀的contig文件
• 我的文件名编码是S1-1，S1-2这样的。

文件夹

• ./concatenate 整理好的序列长度大于2000的
• ./map 放.mmi的索引文件
• ./bam 放.bam的比对文件
• ./sort 排序后的.bam文件
• ./bins 把所有的bin的.fa文件放在一个文件夹
• ./checkmout 放checkm的输出文件，这个文件夹在代码运行前要是空的
• ./bin_cds 放prodigal output

软件位置（已安装和自己安的）

• 软件minimap2 v2.24 
• 软件samtools v1.6  samtools(1) manual page (htslib.org)
• 软件vamb v4.1.1 自己下载
  • 从GitHub - RasmussenLab/vamb: Variational autoencoder for metagenomic binning下载code（vamb-master.zip)→unzip vamb-master.zip → cd vamb-master → pip install -e .
  • avamb/workflow_avamb at avamb_new · RasmussenLab/avamb (github.com)vamb的更多信息可以参考一下这个
• 软件checkM v1.0.12
• 软件prodigal(如果需要比对功能基因）
• 软件diamond(如果需要比对功能基因）
• 软件gtbdtk

GitHub - Ecogenomics/GTDBTk: GTDB-Tk: a toolkit for assigning objective taxonomic classifications to bacterial and archaeal genomes.https://github.com/Ecogenomics/GTDBTk

如果要混样binning，cat合并contigs前参考一下这个:处理contig合并后有序列名相同的问题-CSDN博客

1.取序列长度大于2000的序列（根据推荐）

for file in 2023tire_contig/*.fa; do
filename=$(basename "$file" .fa)
python soft/vamb/vamb-master/src/concatenate.py ./concatenate/$filename.fna.gz $file
done

28个样4分钟（25个线程）

2.创建索引文件mmi

for file in *.fa; do
filename=$(basename "$file" .fa)
python ~/soft/vamb/vamb-master/src/concatenate.py ./concatenate/$filename.fna.gz $file
done

-I 参数用于设置索引文件的最大内存占用量

3. 比对

for file in map/*.mmi; do 
filename=$(basename "$file" .mmi) minimap2 -aI 64g -t 25 -N 5 -ax sr $file --split-prefix mmsplit qc/${filename}_1.fq.gz qc/${filename}_2.fq.gz | samtools view -F 3584 -b --threads 25 > bam/$filename.bam 
done

一个样品9分钟左右，28个样品4个小时（25个线程）。.bam的文件很大，我的文件10G+

4. sort——把得到的bam给sort一下

for file in bam/*.bam; do
filename=$(basename "$file" .bam)
samtools sort --threads 25 -o sort/${filename}.sorted.bam $file
done

一个样品4分钟左右，28个样品2小时（25个线程）。

5.vamb

for file in concatenate/*.fna.gz; do
filename=$(basename "$file" .fna.gz)
vamb -o C --outdir ${filename}-OUT --fasta $file --bamfiles sort/${filename}.sorted.bam --minfasta 200000
done

一个样品6分钟左右，一共3个小时（25个线程）。

6.整理生成的所有bin文件，将前缀命名为样品名，并把它们全部移到binning文件夹

下面这个mv最好改成cp就是了

mkdir -p bins # 创建名为bins的文件夹，如果创过了就不用了 
for dir in S*-OUT/bins/S1; do 
sample_id=$(echo "$dir" | cut -d- -f1-2) # 提取样本ID 
for file in "$dir"/*; do 
mv "$file" "bins/${sample_id}_$(basename "$file")" # 重命名文件并保存到bins文件夹 
done

瞬间完成

解析cut -d- -f1-2（如果使用vamb时把结果文件夹直接名为样品名，就不用这么麻烦了）

-d选项指定了字段的分隔符，这里设置为-；-f选项指定要提取的字段编号或范围，这里设置为1-2，表示提取第一个和第二个字段。

例如，对于字符串S1-1-OUT/bins/S1，使用cut -d- -f1-2命令将返回S1-1，即将-作为分隔符，提取字符串中的前两个字段。

如果不放到文件夹，直接放当前目录用下面这个替换第五行（我用的是下面这个）

mv "$file" "${sample_id}_$(basename "$file")"

7. checkM检测

CheckM ：用于评估从分离株、单细胞或宏基因组中获得的基因组质量的工具。 - 简书 (jianshu.com)

source activate metawrap-env
checkm lineage_wf ./bins ./checkmout -t 25 --pplacer_threads 8 --tab_table -f bin_checkm.tsv

checkm lineage_wf默认的bins后缀是.fna，如果你的文件后缀是.fa，需要增加一个参数（-x fa）

28个样品耗时25分钟（25个线程）

8.把completeness大于50 and contamination小于 10的挑进文件夹bins_90_10

手动从表格里筛选completeness大于50 and contamination小于10的（可以用excel），另存为表格bin_checkm_50_10.tsv。

mkdir -p bins_90_10 # 创建名为bins_90_10的文件夹 
while read prefix _; do 
find ./bins -name "${prefix}*" -exec cp {} bins_50_10/ \; # 复制符合前缀的文件到bins_90_10文件夹 
done < bin_checkm_50_10.tsv

可以把挑出来的这些再checkm一次（但没必要）

checkm lineage_wf ./bins_50_10 ./checkmout -t 25 --pplacer_threads 8 --tab_table -f bin_checkm_50_10.tsv

十分钟左右（25个线程）

9.prodigal 基因预测

其实checkm里面就直接做了这个了，但是文件还得移来移去重命名，我就直接自己prodigal了

for file in bins/*.fna; do

filename=$(basename "$file" .fna) prodigal -i $file -f gff -p meta -o bin_cds/$filename.gff3 -a bin_cds/$filename.faa

done

28分钟

prodigal虽然没有在这一命令里输入线程值，但它会根据提交任务设置的CPU使用数自动分配并加速。

10.diamond比对功能基因

用blastp比对自己的库

11.物种注释

source activate gtdbtk

gtdbtk classify_wf --pplacer_cpus 50 --cpus 50 --genome_dir bins_50_10 --full_tree -x fna --out_dir gtdbtk 
tail -n+2 gtdbtk/gtdbtk.bac120.summary.tsv|cut -f 1-2|sed 's/;/\t/g'|sed '1 s/^/ID\tDomain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies\n/' > gtdbtk/tax.txt

50线程，CPU，1h（感觉提高线程对增速不明显）

进入./gtdbtk/align

gunzip gtdbtk.bac120.user_msa.fasta.gz

这一步的可能报错如下（对于我的427个样品），要注意提交服务器节点的运行内存性能。

WARNING: pplacer requires ~320 GB of RAM to fully load the bacterial tree into memory.

12.对细菌们建树（同样是用gtdbtk，利用gtdbtk classify_wf生成的多序列比对文件）

gtdbtk infer --out_dir infer --cpus 50 --msa_file gtdbtk/align/gtdbtk.bac120.user_msa.fasta

14. 画进化树图

使用生成的gtdbtk.unrooted.tree文件作图，另用一个excel导入进化树图的注释信息。

circle tree (chiplot.online)