技术贴 | 宏基因组Binning（二）质控、分箱、质检、可视化

原创微文，欢迎转发转载。

只做宏基因组太单调？为什么不试试宏基因组Binning呢？一次测序，“宏基因组”+“Binning”两种分析，微生太帮您一站式处理宏基因组难题。现在，微生太免费向所有人分享Binning的整套分析流程，包含：生信分析代码和R语言绘图代码。我们一共设计了7个课时，每周一次，课表（进度）如下。

上一篇介绍了Binning的基础内容、微生太Binning分析的流程、结果报告。从这篇开始将分享Binning分析的代码。主干步骤的代码将展示在PPT中，主程序源代码和R PYTHON PERL脚本可通过微文尾部的二维码添加微生太老师咨询获取。

对Binning分析、R语言绘图感兴趣的朋友千万别错过。错过也没关系，每次课程不仅有回放，还有技术贴带您回顾课程内容。

下面我们一起来回顾第二节课的主体内容吧。

一, Binning分析工作环境

1.1 搭建环境

conda create -y -n metawrap-env python=2.7
conda activate metawrap-env
conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda
conda config --add channels ursky
conda install -y -c ursky metawrap-mg

Metawrap是2018年发表在microbiome上的Binning分析工具，内部包含做binning的所有必备软件，可通过conda安装十分方便。

Github地址：https://github.com/bxlab/metaWRAP

1.2 下载数据库

mkdir MY_CHECKM_FOLDER
cd MY_CHECKM_FOLDER
wget https://data.ace.uq.edu.au/public/CheckM_databases/checkm_data_2015_01_16.tar.gz 
tar -xvf *.tar.gz
rm *.gz

Metawrap需要多个数据库，可以根据自己的需要选择性的下载。

1.3 一键化分析

bash /PATH/metagenome_binning.sh
# 展示帮助文档
bash /PATH/metagenome_binning.sh -i rawdata \
-p 58 -n demo -m 400g -M 0.8 \
-g metadata.txt -G Category
# 微生太Binning一键化分析

图2

Binning的shell主程序可在每节课后联系微生太老师获取。

二、分析内容

图3

本节课要展示和分享的是标准分析中红框中的部分。标准分析的结题报告可通过联系微生太老师获取，

三、过程：生物信息分析

3.1 输入数据

图4

以上是程序默认可识别的输入fastq文件名命名规则和分组表的格式。通过分组可进行Bin丰度差异分析。

3.2 质控

简介：

测序得到的原始数据会存在一定比例的低质量数据，为了保证后续信息分析结果的准确可靠，首先要对原始数据进行质控及宿主过滤，得到有效数据。分析中将使用Kneaddata软件彻底清除原始数据中的Illumina接头序列、低质量的序列片段和较短序列。质控前和质控后，会用FastQC来检测质控的合理性和效果。质控处理后的数据可进一步通过Bowtie2软件比对到宿主的基因组，没有比对到的序列被保留下来做后续分析。

质控、名称格式化

# 1. 质控【批处理】
for i in ./$rawdata/*_R1.fastq; do
    r1=$i
    r2=${i%_R1.fastq}_R2.fastq
    kneaddata -t $thread --max-memory $memory -v \
-i $r1 \
-i $r2 \
-o Bin_all/kneaddata \
--trimmomatic /home/cheng/pipelines/Genome/softwares/Trimmomatic-0.39/ \
--trimmomatic-options "ILLUMINACLIP:/home/cheng/pipelines/MetaGenome/data/adaptor_Illumina.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:50" \
--fastqc /home/cheng/softwares/fastqc/FastQC/ \
--run-fastqc-start \
--run-fastqc-end \
--remove-intermediate-output
    echo -e "\033[32m$i Done...\033[0m"
done
# 2. 格式化clean data名称
for i in Bin_all/kneaddata/*_R1_kneaddata.trimmed.1.fastq; do
    base=${i%_R1_kneaddata.trimmed.1.fastq}
    name=${base##*/}
    mv $i Bin_all/Clean_data/${name}_1.fastq
done
for i in Bin_all/kneaddata/*_R1_kneaddata.trimmed.2.fastq; do
    base=${i%_R1_kneaddata.trimmed.2.fastq}
    name=${base##*/}
    mv $i Bin_all/Clean_data/${name}_2.fastq
done

# 3. 统计rawdata cleandata信息
python3 /home/cheng/huty/softwares/script_bin/summary_fastq.py rawdata Bin_all/summary_rawdata.txt
python3 /home/cheng/huty/softwares/script_bin/summary_fastq.py Bin_all/Clean_data Bin_all/summary_cleandata.txt 

统计质控前后序列信息

图5

说明：

第一列：样品名

InsertSize：文库插入片段大小

SeqStrategy：采用PE150测序策略

Reads/Base/GC：序列数/碱基数/GC含量

MaxLength/MinLength：最大最小序列长度

3.3 分箱

简介：

使用Megahit组装混合到一起的Clean测序数据，得到contigs。分别用Bowtie2和Samtools进行比对和格式转换。得到conitgs深度数据后，使用Metabat2进行宏基因组分箱。设置Metabat2参数“-m 1500”，即使用1500bp以上的conitgs进行分箱。使用Checkm对每个Bin进行质量检查，获得原始Bin的完成度和污染度信息。

分箱准备、分箱

# 1. 分箱准备：合并、组装、索引、比对、sam2bam、sort bam、计算深度
cat Bin_all/Clean_data/*_1.fastq > Bin_all/Bin_align/R1.fastq
cat Bin_all/Clean_data/*_2.fastq > Bin_all/Bin_align/R2.fastq

time megahit -1 Bin_all/Bin_align/R1.fastq -2 Bin_all/Bin_align/R2.fastq -o Bin_all/Assembly_out --out-prefix final -t $thread -m $megahit_mm

time bowtie2-build --large-index -f Bin_all/Assembly_out/final.contigs.fa Bin_all/Assembly_out/final.contigs --threads $thread --quiet

time bowtie2 --mm -1 Bin_all/Bin_align/R1.fastq -2 Bin_all/Bin_align/R2.fastq -p $thread -x Bin_all/Assembly_out/final.contigs -S Bin_all/Bin_align/final.sam

time samtools view -@ $thread -b -S Bin_all/Bin_align/final.sam -o Bin_all/Bin_align/final.bam

time samtools sort -@ $thread -l 9 -O BAM Bin_all/Bin_align/final.bam -o Bin_all/Bin_align/final.sorted.bam

time jgi_summarize_bam_contig_depths --outputDepth Bin_all/Bin_align/final.depth.txt Bin_all/Bin_align/final.sorted.bam

# 2. 分箱
time metabat2 -m 1500 -t $thread -i Bin_all/Assembly_out/final.contigs.fa -a Bin_all/Bin_align/final.depth.txt -o Bin_all/Bin/bin -v

分箱结果

图6

6个样品合并组装分箱一共得到136个bins，类似细菌基因组组装，binning的结果也是有好有坏，下面通过质检可对bins进行筛选。

3.4 质检、筛选

简介：

使用Checkm软件中的lineage_wf功能对分箱结果进行质量检查，获得每个箱子的完成度、污染度、物种分类等信息。筛选完成度大于70%，同时污染度小于10%的Bin用于后续分析。

质检、筛选

# 1. 质检
time checkm lineage_wf -f Bin_all/Bin_quality/checkm.txt -t $thread -x fa Bin_all/Bin/ Bin_all/Bin_quality/ --tmpdir Bin_all/tmp/
# 2. 格式化
grep 'bin’ Bin_all/Bin_quality/checkm.txt | sed 's/^  //' | awk '{print $1,$13,$14}' | sed 's/\ /\t/g'| sed 's/\./\t/' | sort -n -k 2 | sed 's/\t/./' | sort -k 2 -n -r | sed '1i\BinID\tCompleteness\tContamination' > Bin_all/Bin_quality/checkm_raw.txt
# 3. 筛选：污染<10，完成度>70
grep 'bin' Bin_all/Bin_quality/checkm.txt | sed 's/^  //' | awk '{print $1,$13,$14}' | sed 's/\ /\t/g'| sed 's/\./\t/' | sort -n -k 2 | sed 's/\t/./' | awk '{if($2>=70 && $3<=10) print $0}' | sort -k 2 -n -r | sed '1i\BinID\tCompleteness\tContamination' > Bin_all/Bin_quality/checkm_pick.txt
# 4. 提取高质量Bin ID
awk '{print $1}' Bin_all/Bin_quality/checkm_pick.txt | sed '1d' > Bin_all/Bin_quality/checkm_pick_id.txt
# 5. 提取高质量Bin
for i in `cat Bin_all/Bin_quality/checkm_pick_id.txt`; do
    mv Bin_all/Bin/$i.fa Bin_all/Bin_pick/
done

筛选结果

图7

我们从136bins中筛选得到72个bins的污染<10%，完成度>70%。下游的所有分析都是基于这些高质量的Bins。

统计高质量Bin序列信息

# 6. 统计Size和GC含量
stat_bin_gc.py Bin_all/Bin_pick/ Bin_all/Bin_summary/bin.gc.txt
sed '1d' Bin_all/Bin_summary/bin.gc.txt | sed '1 iBinID\tSize\tGC_percent' > Bin_all/Bin_summary/bin.gc_out.txt
# 7. 统计N50和N90
for i in Bin_all/Bin_pick/*.fa; do
    base=${i##*/}
    name=${base%.fa}
    echo $name >> Bin_all/Bin_summary/N50.N90.txt
    perl /home/cheng/huty/softwares/script_bin/N50.N90.pl $i >> Bin_all/Bin_summary/N50.N90.txt
    echo -e "\033[32m$i Done...\033[0m"
done
# 8. 合并checkm n50 n90 size gc%
R CMD BATCH --args /home/cheng/huty/softwares/script_bin/N50.N90.GC.size.R


stat_bin_gc.py N50.N90.pl N50.N90.GC.size.R脚本可通过联系微生太老师获取。

高质量Bin序列信息统计结果

图8

说明：

第一列：bin编号。

第二列：bin完成度（70%~100%）。

第三列：bin污染度（0%~10%）。

第四列：Contig N50是指将一个Bin的Contig按照从大到小的顺序依次相加，当相加的长度达到Contig总长度的一半时，最后一个加上的Contig长度。

第五列：Contig N90是指将一个Bin的Contig按照从大到小的顺序依次相加，当相加的长度达到Contig总长度的90%时，最后一个加上的Contig长度。

第六列：一个Bin的所有Contig长度之和，即该基因组草图的总长度。

第七列：GC含量是指一个Bin中GC碱基占总碱基的比例。

3.5 可视化

简介：

利用可视化的方法直观展示Bin的由来和所含信息。宏基因组分箱的原理是根据序列（contig/scaffold）四核苷酸频率和序列丰度变化模式将序列分成一个个Bin。每个Bin的每个contig都会在此处理（metabat2分箱）的过程中有一个depth深度数据。另外我们利用自己的python脚本计算每个Bin的contig GC含量。有了contigGC含量和depth数据即可进行Bin可视化，绘制每个Bin中每个contig的散点图。MetaWrap中的blobology模块通过bin、contig、样品fastq三组序列数据进行megablast比对等工作也能够进行Bin可视化，但是这种方法非常耗费计算机运行时间。相较之下我们的方法效率非常高，省去大量的多余工作。

统计合并contig GC%和Depth

# 1. 计算每个Bin中每条contig的GC含量
for i in Bin_all/Bin_pick/bin.*.fa; do
    file=${i##*/}
    fold=${file%.fa}
    mkdir Bin_all/Bin_plot/$fold
    python3 /home/cheng/huty/softwares/script_bin/fasta_gc_percent.py $i Bin_all/Bin_plot/$fold/${fold}.gc.txt
    echo -e "\033[32m\t $i gc percent Done...\033[0m"
done
# 2. 获取contig深度数据
cat Bin_all/Bin_align/final.depth.txt | awk -F"\t" 'BEGIN{OFS="\t"}{print $1, $3}' > Bin_all/Bin_plot/all_contig_avgdepth.txt
# 3. 合并contig gc depth
for i in Bin_all/Bin_plot/bin.*; do
    fold=${i##*/}
    Rscript /home/cheng/huty/softwares/script_bin/merge_gc_depth.R $i/${fold}.gc.txt ../all_contig_avgdepth.txt
    echo -e "\033[32mmerge $i gc depth Done...\033[0m"
done
# 4. 合并all bin gc depth data
for i in Bin_all/Bin_plot/bin*; do fold=${i##*/}; \
cat $i/${fold}.gc.depth.txt | sed '1d' >> Bin_all/Bin_plot/all_bin_gc_depth.txt; done


fasta_gc_percent.py merge_gc_depth.R可通过联系微生太老师获取。

统计合并结果

图9

统计结果可视化：R ggplot scatter

# 5. gc depth 可视化
Rscript /home/cheng/huty/softwares/script_bin/bin_scatter.R Bin_all/Bin_plot/all_bin_gc_depth.txt all_bin_gc_depth.pdf
convert -density 400 -quality 200 Bin_all/Bin_plot/all_bin_gc_depth.pdf Bin_all/Bin_plot/all_bin_gc_depth.png

bin_scatter.R脚本：

#!/usr/bin/env Rscript
library(ggplot2)
args = commandArgs(T)
route_file = unlist(strsplit(args[1], "/"))
route = paste(route_file[1:(length(route_file)-1)], collapse="/")
setwd(route)
file_name = route_file[length(route_file)]
data=read.table(file_name, header=T, sep="\t")
result = ggplot(data, aes(x=GC_percent, y=Depth, color=Bin)) +
geom_point(size = 0.5) +
labs(x="GC Percent", y="Average Depth") +
theme(legend.title=element_text(face=“bold”), panel.grid=element_blank(), panel.background=element_rect(color='black’, fill='transparent’))  # 去背景，加black框
ggsave(result, filename=args[2])

说明：

bin_scatter.R脚本内容，思路如下：

1 加载需要的包

2 写参数传递方法

3 写画图参数

4 保存

可视化结果

图10

说明：

横坐标是contig的GC含量；纵坐标是contig depth；一个点代表一个contig，相同颜色的contig来自同一个Bin。

---

---

---