学习资料：

• EasyMetagenome
• De novo genome assembly
• Assembling a metagenome and recovering "genomes" with Anvi'o

使用的软件：

• megahit
• metaSPAdes
• QUAST
• Prodigal

megahit组装

 megahit -1 r3_1.fa  -2 r3_2.fa  -o out2

metaSPAdes精细拼接

    # 精细但使用内存和时间更多，15~65m
    memusg -t 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
    # 23M，contigs体积更大
    seqkit stat temp/metaspades/contigs.fasta

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

注：metaSPAdes支持二、三代混合组装，见附录，此外还有OPERA-MS组装二、三代方案

QUAST评估

There are some general things we can look at, like N50 or largest contig, or fraction of reads that successfully recruit to our assembly, and more.

Quast （Quality Assessment Tool for Genome Assemblies）是用于评估单个细菌或真核生物基因组组装质量的工具。它比较组装结果与参考基因组之间的差异，并生成详细的统计报告。Quast 提供了各种评估指标，如基因组覆盖率、N50 值、误配率、插入/删除错误等，以帮助评估组装的准确性和完整性。详细文档见QUAST 5.2.0 manual。

MetaQUAST 是基于 Quast 开发的工具，专门用于评估宏基因组（metagenome）组装的质量。宏基因组是从环境样品中提取的混合细菌、古菌和真核微生物的基因组。MetaQUAST 能够处理混合的、多物种的组装结果，并提供对宏基因组组装的评估。它支持多个参考基因组、分离的和混合的组装结果，并提供详细的统计信息和可视化报告。

使用 Quast 对给定的四个组装结果进行评估，并将结果保存在指定的输出目录中。评估结果将包括与参考基因组的比对情况、N50 值、误配率、插入/删除错误等统计信息。

quast -o quast_B_cep_out \
	-R reference_genome/BCep_ref.fna \
	-G reference_genome/BCep_ref.gff \
	-l "spades_default, spades_kmers_careful, megahit_default, megahit_min_count_3" \
	-t 4 -m 1000 spades_default_assembly/contigs.fasta \
	spades_kmers_set_careful_assembly/contigs.fasta \
	megahit_default_assembly/final.contigs.fa \
	megahit_min_count_3_assembly/final.contigs.fa

• -o quast_B_cep_out：指定输出目录为 "quast_B_cep_out"，Quast 将在该目录中生成评估结果。
• -R reference_genome/BCep_ref.fna：指定参考基因组的 FASTA 文件路径。
• -G reference_genome/BCep_ref.gff：指定参考基因组的 GFF 文件路径。
• -l "spades_default, spades_kmers_careful, megahit_default, megahit_min_count_3"：指定待评估的组装结果的标签列表，以逗号分隔。
• -t 4：指定线程数为 4，即使用 4 个线程并行运行。
• -m 1000：指定最大内存使用量为 1000 MB。
• spades_default_assembly/contigs.fasta：指定第一个组装结果的 FASTA 文件路径。
• spades_kmers_set_careful_assembly/contigs.fasta：指定第二个组装结果的 FASTA 文件路径。
• megahit_default_assembly/final.contigs.fa：指定第三个组装结果的 FASTA 文件路径。
• megahit_min_count_3_assembly/final.contigs.fa：指定第四个组装结果的 FASTA 文件路径。

第一件要注意的事情是，所有的组装结果都重构了超过98%的参考基因组，我认为这相当出色，但没有一个达到100%。这可能有几个原因，但最有可能的原因是，仅通过短Illumina读段无法组装超过配对读片段长度的重复区域。

我们还可以看到，我们的组装结果与参考基因组中注释的7,705个基因中的约7,550个基因对齐。

通过观察每100 kbp的错误配对数量，我们可以看到每100 kbp的SNV（单核苷酸变异）约为2-4个，这无疑是测序误差、组装误差和实际生物变异的混合体。

在表格中的与参考基因组无关的部分，我们可以看到哪些组装结果以最少的连续序列覆盖了最多的碱基对。这并不意味着一切，但是较少的连续序列覆盖了大约相同数量的碱基对，提供了更多有关连锁性的信息，这可能非常重要。

    quast.py result/megahit/final.contigs.fa -o result/megahit/quast -t 2
    # 生成report文本tsv/txt、网页html、PDF等格式报告
    
    # (可选) megahit和metaspades比较
    quast.py --label "megahit,metapasdes" \
        result/megahit/final.contigs.fa \
        result/metaspades/contigs.fasta \
        -o result/quast
    
    # (可选)metaquast评估，更全面，但需下载相关数据库，受网速影响可能时间很长(我很少成功)
    # metaquast based on silva, and top 50 species genome to access
    time metaquast.py result/megahit/final.contigs.fa -o result/megahit/metaquast

可视化组装结果

下面脚本用适用于anvi'o的方式组织我们的contig，并生成一些基本统计信息，然后使用程序Prodigal识别开放阅读框架（open-reading frames）。

anvi-gen-contigs-database \
	-f spades_kmers_set_careful_assembly/contigs.fasta \
	-o contigs.db \
	-n B_cepacia_assembly

现在我们有了保存序列和一些基本信息的contigs.db文件，我们可以开始添加更多内容。这是灵活性发挥作用的一个方面，但现在我们将继续进行一般anvi'o工作流的一些部分，包括：

• 使用HMMER程序和概要隐藏马尔可夫模型（HMM）来搜索细菌单拷贝基因和核糖体RNA
• 使用NCBI COGs（Clusters of Orthologous Groups）对由Prodigal预测的开放阅读框进行功能注释，可以使用BLAST或DIAMOND进行注释。
• 使用名为Centrifuge的工具对已识别的开放阅读框进行分类学分类。

  # HMM searching for single-copy genes and rRNAs
    # note, if using a newer version of anvi'o, this may need to be 
    # changed to "Bacteria_71" instead of "Campbell_et_al"
    # can see those available by running `anvi-run-hmms -h`
anvi-run-hmms -I Campbell_et_al -c contigs.db -T 4
anvi-run-hmms -I Ribosomal_RNAs -c contigs.db -T 4

  # functional annotation with DIAMOND against NCBI's COGs
anvi-setup-ncbi-cogs -T 4 # only needed the first time
anvi-run-ncbi-cogs -c contigs.db --num-threads 4

  # exporting Prodigal-identified open-reading frames from anvi'o
anvi-get-sequences-for-gene-calls -c contigs.db -o gene_calls.fa

  # setting up and running centrifuge for taxonomy
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/centrifuge/data/p_compressed+h+v.tar.gz
tar -xzvf p_compressed+h+v.tar.gz && rm -rf p_compressed+h+v.tar.gz

centrifuge -f -x p_compressed+h+v gene_calls.fa -S centrifuge_hits.tsv -p 4

  # importing the taxonomy results into our anvi'o contigs database
anvi-import-taxonomy-for-genes -c contigs.db -i centrifuge_report.tsv centrifuge_hits.tsv -p centrifuge

我们现在要添加的最后一件事情是将我们的读取映射信息映射到组装结果上。以下是使用bowtie2生成映射信息并为anvi'o准备的步骤：

  # building bowtie index from our selected assembly fasta file
bowtie2-build spades_kmers_set_careful_assembly/contigs.fasta \
              spades_kmers_set_careful_assembly.btindex
  # mapping our reads
bowtie2 -q -x spades_kmers_set_careful_assembly.btindex \
        -1 BCep_R1_QCd_err_cor.fastq.gz \
        -2 BCep_R2_QCd_err_cor.fastq.gz \
        -p 4 -S spades_kmers_set_careful_assembly.sam
  # converting to a bam file
samtools view -bS spades_kmers_set_careful_assembly.sam > B_cep_assembly.bam
  # sorting and indexing our bam file (can be done with samtools also)
anvi-init-bam B_cep_assembly.bam -o B_cep.bam

然后，我们可以使用anvi-profile程序将这些映射信息与anvi'o集成，该程序生成anvi'o称之为“profile database”的另一种数据库。在这一步骤中涉及到很多内容，我们可以进行各种操作，并在这里进行讨论，但就我们目前的目的而言，主要是为每个contig提供覆盖度信息。

  # integrating the read-mapping info into anvi'o
anvi-profile -i B_cep.bam -c contigs.db -M 1000 -T 4 --cluster-contigs -o B_cep_profiled/

我们刚刚生成并将关于我们的组装的大量信息放入了我们的contigs和profile database中。作为anvi'o，现在我们可以很容易地在需要时访问这些信息的具体细节。

例如，我们刚刚运行的命令之一是anvi-run-hmms，它搜索可能感兴趣的核糖体RNA和单拷贝基因。我们可以使用anvi-get-sequences-for-hmm-hits命令来获取与这些HMMs匹配的序列（添加-h标志可以告诉我们如何使用该程序）。假设我们想要查看被识别为核糖体RNA的内容：

anvi-get-sequences-for-hmm-hits -c contigs.db --hmm-sources Ribosomal_RNAs -o rRNAs.fa

这将所有核糖体RNA的匹配结果写入了一个名为rRNAs.fa的新文件中，如果我们查看该文件，我们会发现找到了2个匹配，一个是16S和一个是23S。为了进一步确认，我们可以快速使用BLAST对它们进行比对，而在这种情况下，我们会感到高兴，因为它们与 B. cepacia ATCC 25416完全一致

类似地，我们可以提取出我们的细菌profile HMM搜索到的所有单拷贝基因，并将其转换为氨基酸序列：

anvi-get-sequences-for-hmm-hits -c contigs.db --hmm-sources Campbell_et_al --get-aa-sequences -o bacterial_SCGs.faa --no-wrap

什么是co-assembly？

"Co-assembly"是指执行一种组装过程，其中输入文件是来自多个样本的reads。这与对每个样本单独进行组装的方式形成对比，后者的输入仅为该个体样本的reads。合并组装具有三个主要优点：1）更高的reads深度（这可以实现更稳健的组装，捕获系统中更多的多样性，但并非总是如此）；2）通过提供一个参考组装方便进行跨样本比较；3）由于差异覆盖度的巨大作用，它可以显著提高从宏基因组中恢复基因组的能。

什么时候使用co-assemble？

虽然合并组装有其优点，但并非在所有情况下都是理想的选择。在生物信息学中，没有一项黄金规则可以确定何时最好将多个样本进行合并组装，而何时最好对每个样本进行独立组装。这可能取决于许多因素（例如样本之间的变异性、群落的多样性、您使用的组装工具、您的整体问题等等），并且可能存在某些数据集，合并组装会得到比独立组装更差的组装结果。因此，在处理每个数据集时，您应该仔细考虑，并且如果您不确定并且有时间/资源，那么尝试并评估结果总是一个好的做法。