生物学知识

从CONCOCT入手理解宏基因组binning

2019-03-06  本文已影响101人  UnderStorm

文章会同步发表于本人 github 和 github homepage

目录

    1. 宏基因组binning简介
    1. binning原理
    • 2.1. 可用于binning的特征
    • 2.2. 从哪些序列下手进行binning？
    1. 详述CONCOCT的binning原理
    1. 拆解CONCOCT的流程
    • 4.1. Assembling Metagenomic Reads
    • 4.2. Cutting up contigs
    • 4.3. Map, Remove Duplicate and Quant Coverage

1. 宏基因组binning简介

Metagenome 组装完成后，我们得到的是成千上万的 contigs，我们需要知道哪些 contigs 来自哪一个基因组，或者都有哪些微生物的基因组。所以需要将 contigs 按照物种水平进行分组归类，称为 "bining"

Supervised binning methods: use databases of already sequenced genomes to label contigs into taxonomic classes

Unsupervised (clustering) methods: look for natural groups in the data

Both supervised and unsupervised methods have two main elements: a metric to define the similarity between a given contig and
a bin, and an algorithm to convert those similarities into assignments

一个很容易想到的策略就是，将组装得到的片段与已知物种的参考基因组进行比对，根据同源性进行归类。然而目前大多数的微生物的基因组还没有测序出来，因此限制了这种方法的可行性。

目前主流的 bining 策略利用的是 contigs 的序列组成特点。

2. binning原理

2.1. 可用于binning的特征

2.2. 从哪些序列下手进行binning？

从原始的clean reads，还是从组装成的contig，还是从预测到的gene，都可以。根据基于聚类的序列类型的不同，暂且分为reads binning， contig binning和 genes binning

比较这三种binning的优劣：

  • contig binning

    由于核酸组成和物种丰度变化模式在越长的序列中越显著和稳定，基于contig binning效果可能更好

  • reads binning

    基于reads binning的优势是可以聚类出宏基因组中丰度非常低的物种

    考虑到在宏基因组组装中reads利用率很低，单样品5Gb测序量情况下，环境样品组装reads利用率一般只有10%左右，肠道样品或极端环境样品组装reads利用率一般能达到30%，这样很多物种，尤其是低丰度的物种可能没有被组装出来，没有体现在gene 或者contig 中，因此基于reads binning 才有可能得到低丰度的物种

    如 Brian Cleary 等 (DOI:10.1038/nbt.3329.Detection) 利用基于 reads binning 的 latent strain analysis 可以聚类出丰度低至0.00001%的菌株。此方法虽然得到更全面的 bins，但低丰度 bins 信息依旧不完整。

  • genes binning

    应用非常广泛

    原因可能是（1）基于genes丰度变化模式进行binning可操作性比较强，宏基因组分析中肯定都会计算gene丰度，一般不会计算contig丰度，gene丰度数据可以信手拈来；（2）基于genes binning有很多可参考的文献，过程也并不复杂，可复制性强；（3）对计算机资源消耗比较低

总体来说应用最广泛的就是基于genes binning 和 contig binning

Naseer Sangwan 等 (DOI: 10.1186/s40168-016-0154-5) 总结了 contig binning 的算法和软件（如下表）

基于Genes abundance binning的一般流程

在宏基因组做完组装和基因预测之后，把所有样品中预测到的基因混合在一起，去冗余得到unique genes集合，对这个unique genes集合进行binning，主要是根据gene在各个样品中的丰度变化模式，计算gene之间的相关性，利用这种相关性进行聚类

该图中的聚类过程类似于K-means聚类：随机选择几个seed genes作为诱饵，计算其他基因丰度分布模式与seed genes的相关性，按照固定的相关性值PCC>0.9，将它们归属于不同seed genes所代表的类，然后在聚好的类内重新选择seed genes，进行迭代，最终聚类得到一个个基因集合，较大的集合（超过700个基因）称为 metagenomic species (MGS)，较小的集合称为 co-abundance gene group (CAG)

基于 binning 结果进行单菌组装：

Sequence reads from individual samples that map to the MGS genes and their contigs are then extracted and used to assembly a draft genome sequence for an MGS

3. 详述CONCOCT的binning原理

结合序列组成特征 (sequence composition) 和跨样本覆盖度特征 (coverage across multiple samples) 进行binning

在进行binning之前需要将所有样本的reads进行混拼 (coassembly) 得到contigs

4. 拆解CONCOCT的流程

4.1. Assembling Metagenomic Reads

# 将多个样本的测序数据fastq文件，按照双端分别进行合并
$ cat $CONCOCT_TEST/reads/Sample*_R1.fa > All_R1.fa
$ cat $CONCOCT_TEST/reads/Sample*_R2.fa > All_R2.fa

# 拼接
$ velveth velveth_k71 71 -fasta -shortPaired -separate All_R1.fa All_R2.fa
$ velvetg velveth_k71 -ins_length 400 -exp_cov auto -cov_cutoff auto

velveth:

takes in a number of sequence files, produces a hashtable, then outputs two files in an output directory (creating it if necessary), Sequences and Roadmaps, which are necessary to velvetg.

语法：

./velveth output_directory hash_length [[-file_format][-read_type] filename]

velvetg:

Velvetg is the core of Velvet where the de Bruijn graph is built then manipulated

4.2. Cutting up contigs

将大片段的contigs (>=20kb)，切成一个个10kb的小片段，当切到尾部只剩不到20kb时，停止切割，以防切得过碎

python $CONCOCT/scripts/cut_up_fasta.py -c 10000 -o 0 -m contigs/velvet_71.fa > contigs/velvet_71_c10K.fa

4.3. Map, Remove Duplicate and Quant Coverage

  1. 使用 Bowtie2 执行 mapping 操作

  2. 用 MarkDuplicates（Picard中的一个工具） 去除 PCR duplicates

  3. 用 BEDTools genomeCoverageBed 基于 mapping 得到的 bam 文件计算每个contigs的coverage

（1） Map, Remove Duplicate

其中1、2步操作可以由CONCOCT中提供的脚本map-bowtie2-markduplicates.sh完成

先要自行建好这些contigs的bowtie2索引

# index for contigs
$ bowtie2-build contigs/velvet_71_c10K.fa contigs/velvet_71_c10K.fa

用map-bowtie2-markduplicates.sh脚本完成 mapping -> remove duplicate

for f in $CONCOCT_TEST/reads/*_R1.fa; do
    mkdir -p map/$(basename $f);
    cd map/$(basename $f);
    bash $CONCOCT/scripts/map-bowtie2-markduplicates.sh -ct 1 -p '-f' $f $(echo $f | sed s/R1/R2/) pair $CONCOCT_EXAMPLE/contigs/velvet_71_c10K.fa asm bowtie2;
    cd ../..;
done
  • -c option to compute coverage histogram with genomeCoverageBed
  • -t option is number of threads
  • -p option is the extra parameters given to bowtie2. In this case -f
  • -k 保留中间文件

随后的5个参数：

  • pair1, the fasta/fastq file with the #1 mates
  • pair2, the fasta/fastq file with the #2 mates
  • pair_name, a name for the pair used to prefix output files
  • assembly, a fasta file of the assembly to map the pairs to
  • assembly_name, a name for the assembly, used to postfix outputfiles
  • outputfolder, the output files will end up in this folder

如果要自己逐步执行第1、2两步，则可以通过以下方式实现：

# Index reference, Burrows-Wheeler Transform
$ bowtie2-build SampleA.fasta SampleA.fasta

# Align Paired end, sort and index
bowtie2 \
    -p 32 \
    -x SampleA.fasta \
    -1 $Data/SampleA.1.fastq \
    -2 $Data/SampleA.2.fastq | \
    samtools sort -@ 18 -O BAM -o SampleA.sort.bam
samtools index SampleA.sort.bam

# Mark duplicates and index
java -Xms32g -Xmx32g -XX:ParallelGCThreads=15 -XX:MaxPermSize=2g -XX:+CMSClassUnloadingEnabled \
    -jar picard.jar MarkDuplicates \
    I=./SampleA.sort.bam \
    O=./SampleA.sort.md.bam \
    M=./SampleA.smd.metrics \
    VALIDATION_STRINGENCY=LENIENT \
    MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000 \
    REMOVE_DUPLICATES=TRUE # 该参数默认为false，即在输出中不过滤duplicate，但是会对这些记录的flag进行修改标记
samtools index ./SampleA.sort.md.bam

（2）Quant Coverage

第3步，计算每个contigs的coverage，用gen_input_table.py脚本

# usage: gen_input_table.py [-h] [--samplenames SAMPLENAMES] [--isbedfiles] fastafile bamfiles [bamfiles ...]
# --samplenames 写有样品名的文件，每个文件名一行
# --isbedfiles  如果在上一步map时运行了genomeCoverageBed，则可以加上此参数后直接用 *smds.coverage文件。如果没运行genomeCoverageBed，则不加此参数，依旧使用bam文件。

$ python $CONCOCT/scripts/gen_input_table.py --isbedfiles \
    --samplenames <(for s in Sample*; do echo $s | cut -d'_' -f1; done) \
    ../contigs/velvet_71_c10K.fa */bowtie2/asm_pair-smds.coverage \
    > concoct_inputtable.tsv

注：

这个脚本可以接受两种类型的输入

  • （1）对bamfiles执行genomeCoverageBed (bedtools genomecov得到的*smds.coverage文件，此时要使用--isbedfiles参数，这样脚本只执行下面提到的第2步操作——计算每条contig的平均depth（又称为这条contig的abundance）；
  • （2）原始的bamfiles，则脚本要执行下面提到的两步操作；

也可以自己写命令逐步实现，这样有利于加深对工具的理解

  1. 计算每条contig的depth分布（histograms）

    $ bedtools genomecov -ibam ./SampleA.smds.bam > ./SampleA.smds.coverage

    bedtools genomecov默认计算histograms，如输出为chr1 0 980 1000，则说明在contig chr1上depth=0的碱基数为980bp，该contig长度为1000bp

    例如：

    $ cat A.bed
chr1  10  20
chr1  20  30
chr2  0   500

$ cat my.genome
chr1  1000
chr2  500

$ bedtools genomecov -i A.bed -g my.genome
chr1   0  980  1000  0.98
chr1   1  20   1000  0.02
chr2   1  500  500   1
genome 0  980  1500  0.653333
genome 1  520  1500  0.346667

    输出格式为：

    • chromosome
    • depth of coverage from features in input file
    • number of bases on chromosome (or genome) with depth equal to column 2
    • size of chromosome (or entire genome) in base pairs
    • size of chromosome (or entire genome) in base pairs

  2. 计算每条contig的平均depth

    有两种计算方法：

    或

    第二种计算方法本质上就是加权平均

    awk 'BEGIN {pc=""} 
{
    c=$1;
    if (c == pc) {
        cov=cov+$2*$5;
    } else {
      print pc,cov;
      cov=$2*$5;
    pc=c}
} END {print pc,cov}' SampleA.smds.coverage | tail -n +2 > SampleA.smds.coverage.percontig

（3）Generate linkage table

接着要构建 linkage per sample between contigs，目前不是很理解它这一步的目的

尝试作简单的理解：

# usage: bam_to_linkage.py [-h] [--samplenames SAMPLENAMES] [--regionlength REGIONLENGTH] [--fullsearch] [-m MAX_N_CORES] [--readlength READLENGTH] [--mincontiglength MINCONTIGLENGTH] fastafile bamfiles [bamfiles ...]
# --samplenames 写有样品名的文件，每个文件名一行
# --regionlength contig序列中用于linkage的两端长度 [默认 500]
# --fullsearch 在全部contig中搜索用于linkage
# -m 最大线程数，每个ban文件对应一个线程
# --readlength untrimmed reads长度 [默认 100]
# --mincontiglength 识别的最小contig长度 [默认 0]

cd $CONCOCT_EXAMPLE/map
python bam_to_linkage.py -m 8 --regionlength 500 --fullsearch --samplenames sample.txt $DATA/SampleA.fasta ./SampleA.smds.bam > SampleA_concoct_linkage.tsv
mv SampleA_concoct_linkage.tsv ../concoct-input

# 输出文件格式
# 共2+6*i列 (i样品数)，依次为contig1、contig2、nr_links_inward_n、nr_links_outward_n、nr_links_inline_n、nr_links_inward_or_outward_n、read_count_contig1_n、read_count_contig2_n
# where n represents sample name. 
# Links只输出一次，如 contig1contig2 输出，则 contig2contig1 不输出

# contig1: Contig linking with contig2
# contig2: Contig linking with contig1
# nr_links_inward: Number of pairs confirming an inward orientation of the contigs -><-
# nr_links_outward: Number of pairs confirming an outward orientation of the contigs <--> 
# nr_links_inline: Number of pairs confirming an outward orientation of the contigs ->->
# nr_links_inward_or_outward: Number of pairs confirming an inward or outward orientation of the contigs. This can be the case if the contig is very short and the search region on both tips of a contig overlaps or the --fullsearch parameter is used and one of the reads in the pair is outside
# read_count_contig1/2: Number of reads on contig1 or contig2. With --fullsearch read count over the entire contig is used, otherwise only the number of reads in the tips are counted.

参考资料：

(1) Quince C, Walker A W, Simpson J T, et al. Shotgun metagenomics, from sampling to analysis[J]. Nature Biotechnology, 2017, 35(9):833.

(2) Nielsen H B, Almeida M, Juncker A S, et al. Identification and assembly of genomes and genetic elements in complex metagenomic samples without using reference genomes[J]. Nature Biotechnology, 2014, 32(8):822-828.

(3) Sangwan N, Xia F, Gilbert J A. Recovering complete and draft population genomes from metagenome datasets[J]. Microbiome, 2016, 4(1):8.

(4) Abubucker, S. et al. Metabolic reconstruction for metagenomic data and its application to the human microbiome. PLoS Comput. Biol. 8, e1002358(2012).

(5) Beaulaurier J, Zhu S, Deikus G, et al. Metagenomic binning and association of plasmids with bacterial host genomes using DNA methylation.[J]. Nature Biotechnology, 2017, 36(1).

(6) Alneberg, J. et al. Binning metagenomic contigs by coverage and composition. Nat. Methods 11, 1144–1146 (2014).

(7) CONCOCT’s documentation

(8) Manual for Velvet

(9) BEDtools官网

(10) 【Yue Zheng博客】宏基因组binning-CONCOCT

(11) AI研习社《数据科学中必须熟知的5种聚类算法》

上一篇 下一篇

猜你喜欢

热点阅读