当前位置： 首页 > news >正文

宏基因组分析流程（Metagenomic workflow）202405｜持续更新

news 2026/2/11 1:32:25

Logs

  • 增加R包pctax内的一些帮助上游分析的小脚本（2024.03.03）
  • 增加Mmseqs2用于去冗余，基因聚类的速度非常快，且随序列量线性增长（2024.03.12）
  • 更新全文细节（2024.05.29）

注意：

  • 篇幅有限，本文一般不介绍具体的软件安装配置，数据库下载等，也不提供示例数据，主要帮助快速了解整个宏基因组的上游分析流程。所以软件安装请参考软件自身介绍，运行时遇到具问题也最好先Google或者去github issue里面问问。

  • 宏基因组分析技术还是在不停的更新迭代，很多更新更好用的软件在不断推出，所以我在这也会持续更新（如果我还一直做这个方向的话）。这里介绍的也只是比较基本的分析流程，但是掌握了之后，自己进一步去做个性化分析也会顺手很多。

  • 完成上游分析后我们可以得到各种物种或者功能的profile，后续一般用python，R进行数据分析和可视化，可参考我的其他博文。

  • 绝大多数这里介绍的软件都是仅支持linux平台的，我们做测序文件上游分析也肯定是在服务器上做，个人PC一般很难满足需求，所以在做这些分析前必须先学习linux基础知识如文件系统，shell脚本编写，软件安装等。安装软件建议使用conda或mamba（新建环境和管理），有很多参考方法。

  • 我使用的集群是slurm作业管理系统，下面的示例脚本也是适用与slurm的，尽量先学习一下slurm的使用再尝试提交作业。如果不用slurm的话可以只参考#############中间的软件部分。

Introduction

宏基因组（Metagenome）是指对一个生态系统中的所有微生物进行DNA分析的过程，可以帮助研究人员了解微生物的多样性、功能和互作关系。

宏基因组的应用非常广泛，包括：

  • 生物多样性研究：通过对宏基因组进行分析，可以了解不同生态系统中微生物的多样性和分布情况。
  • 生态学研究：宏基因组可以帮助研究人员了解微生物在生态系统中的功能、互作关系和生态位等。
  • 生物技术：宏基因组可以用于筛选具有特定功能的微生物，例如，寻找能够降解有害物质的微生物。

宏基因组的分析一般包括以下步骤：

  1. DNA提取与建库。
  2. 高通量测序：使用高通量测序技术对扩增后的DNA进行测序，得到原始序列数据。
  3. 数据清洗和组装：对原始数据进行质量控制、去除低质量序列和冗余序列，将序列拼接成较长的连续序列（contigs）。
  4. 基因注释：将contigs中的基因进行注释，得到基因功能信息。
  5. 数据分析：了解微生物多样性、群落结构、功能特征等信息（更多是指获取了物种丰度表或功能丰度表之后的进一步分析）。
  6. MAGs binning， 进化动态等进一步分析

下图是我常用的一套上游基本流程，当然在面对不同项目时应该有不同的侧重点和适用的分析方法，可以在此基础上添加或修改。

最早的时候，这方面的分析我都是参考刘永鑫老师的EasyMetagenome，现在这套流程也发文章了 (1)，值得参考，对上手16S测序数据或宏基因组数据都很有帮助。

最近（2024.3.3）我整理了一下流程放在我的R包pctax里了，方便整个宏基因组的上下游流程结合：

install.packages("pctax")
#使用micro_sbatch可以快速获得每一步的slurm脚本：
pctax::micro_sbatch(work_dir = "~/work/",step = "fastp")

preprocess

一般把所有样本的测序双端文件放在一个data文件夹下，然后把所有样本名写入到一个samplelist文件中，方便后续各种软件的批量操作。

质控：fastp

测序文件质控是指在下游分析前对测序数据进行质量控制和清理，以确保数据的准确性和可靠性。

fastp是一个非常快速、多功能的测序数据质控工具，支持自动化的质控处理，包括去除低质量序列、剪切接头序列和生成质控报告。

#!/bin/bash
#SBATCH --job-name=fastp
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=2G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################
echo "SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
START=$SLURM_ARRAY_TASK_IDNUMLINES=40 #how many sample in one arraySTOP=$((SLURM_ARRAY_TASK_ID*NUMLINES))
START="$(($STOP - $(($NUMLINES - 1))))"echo "START=$START"
echo "STOP=$STOP"####################
for (( N = $START; N <= $STOP; N++ ))
dosample=$(head -n "$N" $samplelist | tail -n 1)echo $samplemkdir -p data/fastp~/miniconda3/envs/waste/bin/fastp -w 8 -i data/rawdata/${sample}_f1.fastq.gz -o data/fastp/${sample}_1.fq.gz \-I data/rawdata/${sample}_r2.fastq.gz -O data/fastp/${sample}_2.fq.gz -j data/fastp/${i}.json
done##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s

fastp的json文件中统计了测序数据的基本指标，我们可以将其整理为表格：
把所有的.json文件移到一个文件夹里report/下，使用我写的R包pctax读取json文件并整理成表格：

pctax::pre_fastp("report/")

另外，FastQC，Cutadapt，Trimmomatic等也是常用的测序质控工具。

去宿主：bowtie2

去宿主的过程其实就是将序列比对到宿主基因组上，然后没有比对到的序列整合成新文件就是去宿主后的了。宿主基因组需要自己先下载好并用 bowtie2-build 建立索引，以人类为例：

#!/bin/bash
#SBATCH --job-name=rm_human
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=2G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################
echo "SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
START=$SLURM_ARRAY_TASK_IDNUMLINES=40 #how many sample in one arraySTOP=$((SLURM_ARRAY_TASK_ID*NUMLINES))
START="$(($STOP - $(($NUMLINES - 1))))"echo "START=$START"
echo "STOP=$STOP"####################
for (( N = $START; N <= $STOP; N++ ))
dosample=$(head -n "$N" $samplelist | tail -n 1)echo $samplemkdir -p data/rm_human~/miniconda3/envs/waste/bin/bowtie2 -p 8 -x ~/db/humangenome/hg38 -1 data/fastp/${sample}_1.fq.gz \-2 data/fastp/${sample}_2.fq.gz -S data/rm_human/${sample}.sam \--un-conc data/rm_human/${sample}.fq --very-sensitiverm data/rm_human/${sample}.sam
done##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s

另外也可以使用bwa，kneaddata等软件去宿主。

fastq信息统计

可以用一些小工具如FastqCount，seqkit等来统计reads数，碱基数，GC，Q20，Q30等指标：

~/biosoft/FastqCount-master/FastqCount_v0.5 xx.fastq.gzTotal Reads     Total Bases     N Bases Q20     Q30     GC
11568822 (11.57 M)      1702829127 (1.70 G)     0.00%   98.00%  94.00%  54.00%

还可以把cleandata再跑一下fastp😂，但是只看处理前的统计指标。因为fastp确实非常快，纯粹用来统计也行。

reads-based

reads-based宏基因组分析是指直接利用测序读长（reads）进行宏基因组数据分析的方法，而不是先进行基因组组装。该方法通过对短读长进行比对和注释，从中提取功能信息和物种组成，常用于快速评估环境样本中的微生物群落结构和功能潜力。

  • 优点：
    • 快速：无需组装，处理速度较快。
    • 全面：能够检测到低丰度的微生物和功能基因，适合大规模样本分析。
  • 局限性：
    • 片段化：由于不进行组装，分析基于短读长，可能会错过长距离的基因联系和结构信息。
    • 依赖数据库：分析结果依赖于参考数据库的全面性和准确性。

物种注释：kraken2

Kraken2是一个用于对高通量测序数据进行分类和标识物种的软件。它使用参考数据库中的基因组序列来进行分类，并使用k-mer方法来实现快速和准确的分类。

使用Kraken2进行基本分类的简单步骤：

  1. 准备参考数据库：Kraken2需要一个参考数据库，以便对测序数据进行分类。可以直接下载官方构建的标准库，也可以从NCBI、Ensembl或其他数据库下载相应的基因组序列，并使用Kraken2内置的工具来构建数据库。

kraken2-build --standard --threads 24 --db ./

--standard标准模式下只下载5种数据库：古菌archaea、细菌bacteria、人类human、载体UniVec_Core、病毒viral。

  1. 安装Kraken2：可以从Kraken2官方网站下载并安装Kraken2软件。

  2. 运行Kraken2：使用Kraken2对测序数据进行分类需要使用以下命令：

kraken2 --db <path_to_database> <input_file> --output <output_file>

这里，<path_to_database>是参考数据库的路径，<input_file>是需要进行分类的输入文件，<output_file>是输出文件的名称。Kraken2将输出一个分类报告文件和一个序列文件。

需要注意的是kraken运行至少要提供数据库大小的内存大小（运行内存），因为它会把整个数据库载入内存后进行序列的注释，所以如果发现无法载入数据库的报错，可以尝试调大内存资源。

#!/bin/bash
#SBATCH --job-name=kraken2
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=40G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################
echo "SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
START=$SLURM_ARRAY_TASK_IDNUMLINES=40 #how many sample in one arraySTOP=$((SLURM_ARRAY_TASK_ID*NUMLINES))
START="$(($STOP - $(($NUMLINES - 1))))"####################
for (( N = $START; N <= $STOP; N++ ))
dosample=$(head -n "$N" $samplelist | tail -n 1)echo $samplemkdir -p result/kraken2~/miniconda3/envs/waste/bin/kraken2 --threads 8 \--db ~/db/kraken2/stand_krakenDB \--confidence 0.05 \--output result/kraken2/${sample}.output \--report result/kraken2/${sample}.kreport \--paired \data/rm_human/${sample}.1.fq \data/rm_human/${sample}.2.fq
done##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s

另外，kraken数据库是可以自己构建的，所以适用于各种项目的物种注释，我做的比较多的是环境样本的宏基因组，就可能需要更全面的物种数据库（甚至除了各种微生物，还要动植物数据等），我们实验室的WX师姐收集构建了一个超大的物种库。

kraken软件运行时载入数据库是一个十分耗时的步骤，而每条序列的鉴定时间差不多，所以我们可以将很多样本的fastq文件合并成一个大文件后输入kraken注释，之后再按照序列的数量拆分结果文件，这样多个样本也只需要载入一次数据库，节省时间，以下是用自定义的kraken2M.py脚本实现。

#!/bin/bash
#SBATCH --job-name=kraken2M
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/kraken/%x_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/kraken/%x_%a.err
#SBATCH --time=14-00:00:00
#SBATCH --partition=mem
#SBATCH --cpus-per-task=32
#SBATCH --mem-per-cpu=100G
echo start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################mkdir -p result/krakenpython /share/home/jianglab/shared/krakenDB/K2ols/kraken2M.py -t 16 \-i data/rm_human \-c 0.05 \-s .1.fq,.2.fq \-o result/kraken \-d /share/home/jianglab/shared/krakenDB/mydb2 \-k ~/miniconda3/envs/waste/bin/kraken2 \-kt /share/home/jianglab/shared/krakenDB/K2ols/KrakenToolsmkdir -p result/kraken/kreportmv result/kraken/*_kreport.txt result/kraken/kreport/python ~/db/script/format_kreports.py -i result/kraken/kreport \-kt /share/home/jianglab/shared/krakenDB/K2ols/KrakenTools --save-name2taxid##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
Kraken输出格式

标准输出格式output文件(五列表):

  • C/U代表分类classified或非分类unclassifed
  • 序列ID
  • 物种注释
  • 比序列注释的区域，如98|94代表左端98bp，右端94bp比对至数据库
  • LCA比对结果，如”562:13 561:4”代表13 k-mer比对至物种#562，4 k-mer比对至#561物种

报告输出格式report文件(包括6列，方便整理下游分析):

  1. 百分比
  2. count
  3. count最优
  4. (U)nclassified, ®oot, (D)omain, (K)ingdom, §hylum, ©lass, (O)rder, (F)amily, (G)enus, or (S)pecies. “G2”代表位于属一种间
  5. NCBI物种ID
  6. 科学物种名

常用的物种丰度表格式除了kraken report，还有mpa，spf，krona等格式，关于kraken结果的整理以及格式转换方式，有一些现成的脚本或者自己写。

KrakenTools (jhu.edu) 就是一套很好用的kraken工具包，其中常用的有：

  1. extract_kraken_reads.py
    此程序提取读取在任何用户指定的分类id处分类的内容。用户必须指定Kraken输出文件、序列文件和至少一个分类法ID。下面指定了其他选项。截至2021年4月19日，此脚本与KrakenUniq/Kraken2Uniq报告兼容。

  2. combine_kreports.py
    这个脚本将多个Kraken报告合并到一个合并的报告文件中。

  3. kreport2krona.py
    这个程序需要一个Kraken报告文件，并打印出一个krona兼容的文本文件，换成krona文件好画图。

转换后的krona使用 ktImportText MYSAMPLE.krona -o MYSAMPLE.krona.html得到网页可视化结果。

  1. kreport2mpa.py
    这个程序需要一个Kraken报告文件，并打印出一个mpa (MetaPhlAn)风格的文本文件。

  2. combine_mpa.py
    这个程序合并来自kreport2mpa.py的多个输出。

物种+功能：HUMAnN

HUMAnN（The HMP Unified Metabolic Analysis Network）是一款用于分析人类微生物组的功能和代谢能力的工具。

它通过将宏基因组序列与参考基因组数据库比对，利用MetaCyc代谢通路数据库和UniRef蛋白质序列数据库，分析微生物组在功能和代谢通路水平上的组成和活性。HUMAnN2还提供了多样性分析、关联分析和可视化工具，可用于深入研究人类微生物组对宿主健康的影响和治疗策略的制定等方面。

HUMAnN是由美国国家人类微生物组计划（HMP）开发的，目前最新版本为HUMAnN3，于2020年发布。与HUMAnN2相比，HUMAnN3改进了基因家族注释的方法，提高了注释精度和速度，并提供了新的功能和工具，如功能韧度分析、代谢指纹识别和多样性分析等。

但是HUMAnN的数据库基本都是与人相关的微生物，比较适合做各种人体微生物组（肠道，肺部，口腔，皮肤等等），对于环境样本可能unclassified比较多。

HUMAnN要求双端序列合并的文件作为输入，for循环根据实验设计样本名批量双端序列合并。

  • 物种组成调用MetaPhlAn3, bowtie2比对至核酸序列库，解决有哪些微生物存在的问题；
  • 功能组成为humann3调用diamond比对至蛋白库，解决这些微生物参与哪些功能通路的问题；
cd data/rm_human
for i in `cat ~/work/asthma/samplelist`
do
echo $i
cat ${i}_f1.fastq ${i}_r2.fastq >${i}_paired.fastq
done

#!/bin/bash
#SBATCH --job-name=humann
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --nodes=1
#SBATCH --tasks-per-node=1
#SBATCH --cpus-per-task=1
#SBATCH --mem-per-cpu=2G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################
echo "SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
START=$SLURM_ARRAY_TASK_IDNUMLINES=40 #how many sample in one arraySTOP=$((SLURM_ARRAY_TASK_ID*NUMLINES))
START="$(($STOP - $(($NUMLINES - 1))))"#set the min and max
if [ $START -lt 1 ]
thenSTART=1
fi
if [ $STOP -gt 40 ]
thenSTOP=40
fiecho "START=$START"
echo "STOP=$STOP"####################
for (( N = $START; N <= $STOP; N++ ))
dosample=$(head -n "$N" $samplelist | tail -n 1)echo $samplemkdir -p data/pairedcat data/rm_human/${sample}.1.fq data/rm_human/${sample}.2.fq > data/paired/${sample}_paired.fqmkdir -p result/humann~/miniconda3/envs/waste/bin/humann --input data/paired/${sample}_paired.fq \--output result/humann/ --threads 24ln result/humann/${sample}_paired_humann_temp/${sample}_paired_metaphlan_bugs_list.tsv result/humann/rm -rf result/humann/${sample}_paired_humann_temp
done##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s

mkdir -p result/metaphlan2
## 合并、修正样本名
merge_metaphlan_tables2.py \result/humann/*_metaphlan_bugs_list.tsv | \sed 's/_metaphlan_bugs_list//g' \> result/metaphlan2/taxonomy.tsv

contigs-based

contigs-based宏基因组分析是指先将测序读长（reads）组装成较长的序列片段（contigs），然后对这些contigs进行进一步的分析。这种方法可以提供更长的基因组片段，有助于更准确地进行基因注释和微生物群落结构分析。

  • 优点：
    • 更高的解析度：通过组装生成较长的序列片段，可以更准确地进行基因和基因组注释。
    • 结构信息：能够获得长距离的基因联系和基因组结构信息，有助于更全面的功能和进化分析。
  • 局限性：
    • 计算资源需求高：组装过程需要较高的计算资源和时间。
    • 复杂性：组装和binning步骤增加了分析的复杂性，需要更多的技术经验和工具支持。

组装（assembly）

组装（assembly）是指将测序产生的短读长（reads）拼接成更长的连续序列（contigs）的过程。这个过程在基因组学和宏基因组学研究中至关重要，因为它能够将片段化的信息整合成更完整的基因组序列，便于后续的功能和分类分析。

注意assembly有单样本拼装和混拼等策略，可以自行根据计算资源和研究目的选择。

  • 单样本拼装（Single-sample Assembly）

定义：单样本拼装是指对来自单一环境或条件下的样本进行组装。这种策略适用于相对简单的样本，或者希望单独分析每个样本的基因组组成。

优点：

  1. 独立分析：每个样本单独组装，有助于深入分析样本特定的基因组特征。
  2. 数据简单：适用于数据复杂度较低的样本，组装结果较为清晰。

缺点：

  1. 资源需求高：对每个样本进行独立组装，计算资源和时间需求较高。
  2. 有限的覆盖度：对于低丰度微生物，单样本组装可能无法提供足够的覆盖度和组装质量。

适用场景：

  • 环境样本较为简单的研究。

  • 需要对每个样本进行独立比较和分析的研究。

  • 混拼策略（Co-assembly）

定义：混拼是指将多个样本的数据合并在一起进行组装。这种策略适用于复杂的宏基因组样本，通过整合多个样本的数据，可以提高组装的覆盖度和质量。

优点：

  1. 提高覆盖度：合并多个样本的数据，增加了序列的覆盖度，有助于更完整的组装。
  2. 处理复杂样本：适用于复杂的环境样本，能够捕捉到更多的低丰度微生物和基因组信息。

缺点：

  1. 复杂性增加：混拼的数据复杂度高，组装和后续分析更为复杂。
  2. 样本间差异模糊化：样本间的独特特征可能在混拼过程中被模糊化，影响对个体样本的具体分析。

适用场景：

  • 环境样本复杂，包含多种微生物群落。
  • 需要提高低丰度物种的检测能力。
  • 目标是获取整体微生物群落的综合信息。
megahit

Megahit是一个用于对高通量测序数据进行de novo组装的软件。
它使用了一种基于短读比对和图形构建的算法来组装基因组，可以高效地处理大规模的数据集。
以下是Megahit的一些优点和适用情况：

  1. 速度快：Megahit的算法非常高效，可以处理大规模的数据集，通常比其他de novo组装工具更快。
  2. 高质量的组装：Megahit在组装结果的连通性和准确性方面表现优异，尤其在处理高GC含量基因组时效果显著。
  3. 适用于不同类型的测序数据：Megahit支持多种不同类型的测序数据，包括 Illumina HiSeq/MiSeq、IonTorrent和PacBio等平台。
  4. 易于使用：Megahit具有简单的命令行语法，方便用户进行组装操作，且具有中断点，避免失败后全部重跑。
#!/bin/bash
#SBATCH --job-name=megahit
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=10G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################
echo "SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
START=$SLURM_ARRAY_TASK_IDNUMLINES=40 #how many sample in one arraySTOP=$((SLURM_ARRAY_TASK_ID*NUMLINES))
START="$(($STOP - $(($NUMLINES - 1))))"echo "START=$START"
echo "STOP=$STOP"####################
for (( N = $START; N <= $STOP; N++ ))
dosample=$(head -n "$N" $samplelist | tail -n 1)echo $sample#single samplemkdir -p result/megahitmkdir -p result/megahit/contigs~/miniconda3/envs/waste/bin/megahit -t 8 -1 data/rm_human/${sample}.1.fq \-2 data/rm_human/${sample}.2.fq \-o result/megahit/${sample} --out-prefix ${sample}sed -i "/>/s/>/>${sample}_/" result/megahit/${sample}/${sample}.contigs.famv result/megahit/${sample}/${sample}.contigs.fa result/megahit/contigs/
done##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s

混拼的话：

#!/bin/bash
#SBATCH --job-name=megahit2
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=mem
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=100G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`#####################multi-sample混拼#需要大内存，建议3倍fq.gz内存mkdir -p data/com_readmkdir -p result/megahit2for i in `cat samplelist`doecho ${i}cat data/rm_human/${i}.1.fq >> data/com_read/R1.fqcat data/rm_human/${i}.2.fq >> data/com_read/R2.fqdone~/miniconda3/envs/waste/bin/megahit -t 8 -1 data/com_read/R1.fq \-2 data/com_read/R2.fq \-o result/megahit2/ --out-prefix multi_sample##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
metaSPAdes

SPAdes 是一款多功能工具包，专为测序数据的组装和分析而设计。 SPAdes 主要是为 Illumina 测序数据而开发的，但也可用于 IonTorrent。大多数 SPAdes 管道支持混合模式，即允许使用长读取（PacBio 和 Oxford Nanopore）作为补充数据。

而metaSPAdes是目前宏基因组领域组装指标较好的软件，尤其在株水平组装优势明显。相比Megahit表现出更好的基因长度和读长比较率，但是更好的组装质量也伴随着更长时间和内存消耗，同时也有错误组装上升的风险。

#!/bin/bash
#SBATCH --job-name=asthma_megahit
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/megahit/log/%x_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/megahit/log/%x_%a.err
#SBATCH --array=1-33
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=32echo start: `date +'%Y-%m-%d %T'`
start=`date +%s`
echo "SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
sample=$(head -n $SLURM_ARRAY_TASK_ID ~/work/asthma/data/namelist | tail -1)
#sample=$(head -n 1 namelist | tail -1)
echo handling: $sample	
####################
for (( N = $START; N <= $STOP; N++ ))
dosample=$(head -n "$N" $samplelist | tail -n 1)echo $sample#single samplemkdir -p result/metaspadesmkdir -p result/metaspades/contigs~/miniconda3/envs/metawrap/bin/metaspades.py –t 8 -k 21,33,55,77,99,127 --careful \-1 data/rm_human/${sample}.1.fq \-2 data/rm_human/${sample}.2.fq \-o result/metaspades/${sample}sed -i "/>/s/>/>${sample}_/" result/metaspades/${sample}/scaffolds.fastamv result/metaspades/${sample}/scaffolds.fasta result/metaspades/contigs/
done
####################
echo end: `date +'%Y-%m-%d %T'`
end=`date +%s`
echo TIME:`expr $end - $start`s

组装评估：QUAST

QUAST是一个质量评估工具。 QUAST可以使用参考基因组以及不使用参考基因组来评估装配。 QUAST生成详细的报告，表格和图解，以显示装配的不同方面。

基因预测：Prodigal

基因预测是指在基因组序列（如contigs）中识别出编码区序列（CDS, Coding DNA Sequences）的过程。这个过程对于理解基因组功能和注释基因组中的基因是至关重要的。常用的软件有Prodigal，GeneMark（GeneMarkS用于细菌，GeneMark-ES用于真核生物），Glimmer等。

这里以Prodigal为例：
输入文件：拼装好的序列文件 megahit/final.contigs.fa
输出文件：prodigal预测的基因序列 prodigal/gene.fa

prodigal不支持多线程运行，所以我们可以自行分割序列文件调用多个prodigal程序分别跑实现伪多线程。

#!/bin/bash
#SBATCH --job-name=prodigal
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=1
#SBATCH --mem-per-cpu=2G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################
echo "SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
START=$SLURM_ARRAY_TASK_IDNUMLINES=40 #how many sample in one arraySTOP=$((SLURM_ARRAY_TASK_ID*NUMLINES))
START="$(($STOP - $(($NUMLINES - 1))))"echo "START=$START"
echo "STOP=$STOP"####################
for (( N = $START; N <= $STOP; N++ ))
dosample=$(head -n "$N" $samplelist | tail -n 1)echo $samplemkdir -p result/prodigalmkdir -p result/prodigal/gene_famkdir -p result/prodigal/gene_gff~/miniconda3/envs/waste/bin/prodigal -i result/megahit/contigs/${sample}.contigs.fa \-d result/prodigal/gene_fa/${sample}.gene.fa \-o result/prodigal/gene_gff/${sample}.gene.gff \-p meta -f gffmkdir -p result/prodigal/fullgene_fa## 提取完整基因(完整片段获得的基因全为完整，如成环的细菌基因组)grep 'partial=00' result/prodigal/gene_fa/${sample}.gene.fa | cut -f1 -d ' '| sed 's/>//' > result/prodigal/gene_fa/${sample}.fullidseqkit grep -f result/prodigal/gene_fa/${sample}.fullid result/prodigal/gene_fa/${sample}.gene.fa > result/prodigal/fullgene_fa/${sample}.gene.fadone##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s

去冗余

基因集去冗余是指在预测到的基因集中去除重复的基因或高度相似的基因，以获得一个非冗余的基因集合。这一步骤在宏基因组分析中非常重要，因为宏基因组样本中通常包含来自多种微生物的重复或相似序列，去冗余可以简化数据集并提高后续分析的效率和准确性。常用软件有CD-HIT，MMseqs2，UCLUST（集成在QIIME和USEARCH中的工具）等。

上面产生了n个样本的基因预测结果文件，gene.fa文件要先想办法整合为一个文件再去去冗余。

CD-HIT

之前大家常用的是CD-HIT来去冗余：

#!/bin/bash
#SBATCH --job-name=cluster
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --nodes=1
#SBATCH --tasks-per-node=1
#SBATCH --cpus-per-task=1
#SBATCH --mem-per-cpu=2G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################echo "start merge"cat result/prodigal/gene_fa/*.gene.fa > result/prodigal/all.gene.facat result/prodigal/fullgene_fa/*.gene.fa > result/prodigal/all.fullgene.faecho "done"~/miniconda3/envs/waste/bin/cd-hit-est -i result/prodigal/all.gene.fa \-o result/NR/nucleotide.fa \-aS 0.9 -c 0.9 -G 0 -g 0 -T 0 -M 0mv result/NR/nucleotide_rep_seq.fasta result/NR/nucleotide.farm result/NR/nucleotide_all_seqs.fasta~/miniconda3/envs/waste/bin/seqkit translate result/NR/nucleotide.fa > result/NR/protein.fased 's/\*//g' result/NR/protein.fa > result/NR/protein_no_end.farm result/NR/protein.fa##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
Mmseqs2

最近发现了Mmseqs2这个神器，这个聚类要比CD-HIT快非常多，其他应用参见Mmseqs2的基础使用：

#!/bin/bash
#SBATCH --job-name=cluster
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --nodes=1
#SBATCH --tasks-per-node=1
#SBATCH --cpus-per-task=1
#SBATCH --mem-per-cpu=2G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################echo "start merge"cat result/prodigal/gene_fa/*.gene.fa > result/prodigal/all.gene.facat result/prodigal/fullgene_fa/*.gene.fa > result/prodigal/all.fullgene.faecho "done"mkdir -p result/NRmmseqs easy-linclust result/prodigal/all.gene.fa result/NR/nucleotide mmseqs_tmp \--min-seq-id 0.9 -c 0.9 --cov-mode 1  --threads 8mv result/NR/nucleotide_rep_seq.fasta result/NR/nucleotide.farm result/NR/nucleotide_all_seqs.fasta~/miniconda3/envs/waste/bin/seqkit translate result/NR/nucleotide.fa > result/NR/protein.fased 's/\*//g' result/NR/protein.fa > result/NR/protein_no_end.farm result/NR/protein.fa##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s

fasta信息统计

使用seqkit统计刚刚得到的一些关键fasta文件的基础信息。

test_file=`head -n1 $samplelist`
if [ -f result/megahit/contigs/${test_file}.contigs.fa ]; then
~/miniconda3/envs/waste/bin/seqkit stats result/megahit/contigs/*.contigs.fa >result/megahit/contig_stats
fi
if [ -f result/prodigal/gene_fa/${test_file}.gene.fa ]; then
~/miniconda3/envs/waste/bin/seqkit stats result/prodigal/gene_fa/*.gene.fa >result/prodigal/gene_fa_stats
fi
if [ -f result/prodigal/fullgene_fa/${test_file}.gene.fa ]; then
~/miniconda3/envs/waste/bin/seqkit stats result/prodigal/fullgene_fa/*.gene.fa >result/prodigal/fullgene_fa_stats
fi
if [ -f result/NR/nucleotide.fa ]; then
~/miniconda3/envs/waste/bin/seqkit stats result/NR/nucleotide.fa >result/NR/nucleotide_stat
fi

基因定量：salmon

  1. 建立索引
#!/bin/bash
#SBATCH --job-name=salmon-index
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=4
#SBATCH --mem-per-cpu=10G
echo start: `date +"%Y-%m-%d %T"`
start=`date +%s`###################### 建索引, -t序列, -i 索引# 大点内存mkdir -p result/salmon~/miniconda3/envs/waste/share/salmon/bin/salmon index \-t result/NR/nucleotide.fa \-p 4 \-i result/salmon/index##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
  1. 对每个样本定量
#!/bin/bash
#SBATCH --job-name=salmon-quant
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=4
#SBATCH --mem-per-cpu=2G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################
echo "SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
START=$SLURM_ARRAY_TASK_IDNUMLINES=40 #how many sample in one arraySTOP=$((SLURM_ARRAY_TASK_ID*NUMLINES))
START="$(($STOP - $(($NUMLINES - 1))))"echo "START=$START"
echo "STOP=$STOP"####################
## 输入文件：去冗余后的基因和蛋白序列：NR/nucleotide.fa
## 输出文件：Salmon定量后的结果：salmon/gene.count; salmon/gene.TPM
## 定量，l文库类型自动选择，p线程，--meta宏基因组模式
for (( N = $START; N <= $STOP; N++ ))
dosample=$(head -n "$N" $samplelist | tail -n 1)echo $samplemkdir -p result/salmon/quant~/miniconda3/envs/waste/share/salmon/bin/salmon quant \--validateMappings \-i result/salmon/index -l A -p 4 --meta \-1 data/rm_human/${sample}.1.fq \-2 data/rm_human/${sample}.2.fq \-o result/salmon/quant/${sample}.quant
done##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
  1. 合并各样本结果
#!/bin/bash
#SBATCH --job-name=salmon-merge
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=1
#SBATCH --mem-per-cpu=2G
samplelist=/share/home/jianglab/pengchen/work/asthma/samplelistecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################ls result/salmon/quant/*.quant/quant.sf |awk -F'/' '{print $(NF-1)}' |sed 's/.quant//' >tmp_finisheddiff_rows -f samplelist -s tmp_finished >tmp_diff# 计算结果的行数line_count=$( cat tmp_diff| wc -l)# 检查行数是否大于1，如果是则终止脚本if [ "$line_count" -gt 1 ]; thenecho 'sf文件和samplelist数量不一致'exit 1fi##mapping ratescat result/salmon/quant/*.quant/logs/salmon_quant.log |grep 'Mapping rate = '|awk '{print $NF}'> tmppaste samplelist tmp > result/salmon/mapping_rates## combine~/miniconda3/envs/waste/share/salmon/bin/salmon quantmerge \--quants result/salmon/quant/*.quant \-o result/salmon/gene.TPM~/miniconda3/envs/waste/share/salmon/bin/salmon quantmerge \--quants result/salmon/quant/*.quant \--column NumReads -o result/salmon/gene.countsed -i '1 s/.quant//g' result/salmon/gene.*##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s

功能基因注释

功能基因注释是对非冗余基因集进行分类和功能描述的过程，以便理解这些基因在生物体内可能的角色和作用。

上一步已经有了所有的基因和每个样本所有基因的read count定量结果，
我们只需要对上一步的基因序列（或蛋白质序列）进行不同数据库的注释，
然后合并注释结果得到的就是功能丰度表。

很多数据库自带软件都是用diamond比对，如果没有专用软件的数据库我们也可以自己用diamond比对。

diamond选择–outfmt 6的输出结果和blastp格式一样。

eggNOG (COG/KEGG/CAZy)

EggNOG数据库收集了COG（Clusters of Orthologous Groups of proteins，直系同源蛋白簇）,构成每个COG的蛋白都是被假定为来自于一个祖先蛋白，因此是orthologs或者是paralogs。通过把所有完整基因组的编码蛋白一个一个的互相比较确定的。在考虑来自一个给定基因组的蛋白时，这种比较将给出每个其他基因组的一个最相似的蛋白（因此需要用完整的基因组来定义COG），这些基因的每一个都轮番地被考虑。如果在这些蛋白（或子集）之间一个相互的最佳匹配关系被发现，那么那些相互的最佳匹配将形成一个COG。这样，一个COG中的成员将与这个COG中的其他成员比起被比较的基因组中的其他蛋白更相像。

EggNOG里面已经包含了GO，KEGG，CAZy等。

## 下载常用数据库，注意设置下载位置
mkdir -p ${db}/eggnog5 && cd ${db}/eggnog5
## -y默认同意，-f强制下载，eggnog.db.gz 7.9G+4.9G
download_eggnog_data.py -y -f --data_dir ./## 下载方式2(可选)：链接直接下载
wget -c http://eggnog5.embl.de/download/emapperdb-5.0.0/eggnog.db.gz ## 7.9G
wget -c http://eggnog5.embl.de/download/emapperdb-5.0.0/eggnog_proteins.dmnd.gz ## 4.9G
gunzip *.gz

#!/bin/bash
#SBATCH --job-name=eggnog
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=10G
echo start: `date +"%Y-%m-%d %T"`
start=`date +%s`
###################### 大点内存, 数据库有50G左右source ~/miniconda3/etc/profile.d/conda.shconda activate func## diamond比对基因至eggNOG 5.0数据库, 1~9h，默认diamond 1e-3mkdir -p result/eggnogemapper.py --no_annot --no_file_comments --override \--data_dir ~/db/eggnog5 \-i result/NR/protein_no_end.fa \--cpu 8 -m diamond \-o result/eggnog/protein## 比对结果功能注释, 1hemapper.py \--annotate_hits_table result/eggnog/protein.emapper.seed_orthologs \--data_dir ~/db/eggnog5 \--cpu 8 --no_file_comments --override \-o result/eggnog/output## 添表头, 1列为ID，9列KO，16列CAZy，21列COG，22列描述sed '1 i Name\tortholog\tevalue\tscore\ttaxonomic\tprotein\tGO\tEC\tKO\tPathway\tModule\tReaction\trclass\tBRITE\tTC\tCAZy\tBiGG\ttax_scope\tOG\tbestOG\tCOG\tdescription' \result/eggnog/output.emapper.annotations \> result/eggnog/eggnog_anno_output##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
dbCAN2（碳水化合物）
## dbCAN2 http://bcb.unl.edu/dbCAN2
## 创建数据库存放目录并进入
mkdir -p ${db}/dbCAN2 && cd ${db}/dbCAN2
## 下载序列和描述
wget -c http://bcb.unl.edu/dbCAN2/download/CAZyDB.07312020.fa
wget -c http://bcb.unl.edu/dbCAN2/download/Databases/CAZyDB.07302020.fam-activities.txt
## 备用数据库下载地址并解压 
#wget -c http://210.75.224.110/db/dbcan2/CAZyDB.07312020.fa.gz
#gunzip CAZyDB.07312020.fa.gz
## diamond建索引
time diamond makedb \--in CAZyDB.07312020.fa \--db CAZyDB.07312020

#!/bin/bash
#SBATCH --job-name=cazy
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=2G
echo start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################mkdir -p result/dbcan2diamond blastp   \--db ~/db/dbcan2/CAZyDB.07312020  \--query result/NR/protein_no_end.fa \--threads 8 -e 1e-5 --outfmt 6 \--max-target-seqs 1 --quiet \--out result/dbcan2/gene_diamond.f6##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
CARD（抗生素抗性基因ARGs）

RGI软件Github主页: https://github.com/arpcard/rgi，可以参考之前写的Antibiotics resistance gene。

#!/bin/bash
#SBATCH --job-name=rgi
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=2Gecho start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################source ~/miniconda3/etc/profile.d/conda.shconda activate rgimkdir -p result/card~/miniconda3/envs/rgi/bin/rgi main \--input_sequence result/NR/protein_no_end.fa \--output_file result/card/protein \--input_type protein --num_threads 8 \--clean --alignment_tool DIAMOND # --low_quality #partial genes##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
VFDB（毒力因子）

官网地址：http://www.mgc.ac.cn/VFs/

在官网下载数据库时，带有setA 的库为VFDB数据库核心库(set A)，而setB为全库(setB), 其中setA仅包含经实验验证过的毒力基因，而setB则在setA的基础上增加了预测的毒力基因，选择好数据库后，直接用blast/diamond即可完成注释。

#!/bin/bash
#SBATCH --job-name=vfdb
#SBATCH --output=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error=/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array=1
#SBATCH --partition=cpu
#SBATCH --cpus-per-task=8
#SBATCH --mem-per-cpu=2G
echo start: `date +"%Y-%m-%d %T"`
start=`date +%s`####################mkdir -p result/vfdbdiamond blastp   \--db ~/db/VFDB/VFDB_setB_pro \--query result/NR/protein_no_end.fa \--threads 8 -e 1e-5 --outfmt 6 \--max-target-seqs 1 --quiet \--out result/vfdb/gene_diamond.f6##############
echo end: `date +"%Y-%m-%d %T"`
end=`date +%s`
echo TIME:`expr $end - $start`s
其他数据库
  • NCBI NR（Non-redundant）:包含广泛的已知蛋白序列，是常用的参考数据库。
  • Pfam：包含蛋白质家族和功能域的信息，通过Hidden Markov Models (HMMs) 描述蛋白质功能域。
  • 元素循环数据库：专注于环境中的元素（如碳、氮、硫、磷等）的生物地球化学循环相关基因和途径注释。如NCycDB，SCycDB，PCycDB，MCycDB等

这些数据库我们可以自行下载蛋白数据库，然后用blast建库比对。

功能注释合并

写一个python脚本，将表1（基因-功能的对应表）与表2（基因丰度表）合并，即不同基因可能注释到相同功能，把它们的丰度加在一起得到新表3（功能丰度表）。

可以参考https://github.com/YongxinLiu/EasyMicrobiome/blob/master/script/summarizeAbundance.py。

mkdir -p result/summ_table
if [ -f result/eggnog/eggnog_anno_output ]; then
# 汇总，9列KO，16列CAZy按逗号分隔，21列COG按字母分隔，原始值累加
~/db/script/summarizeAbundance.py \\
-i result/salmon/gene.count \\
-m result/eggnog/eggnog_anno_output \\
-c '9,16,21' -s ',+,+*' -n raw \\
-o result/summ_table/eggnog
sed -i 's/^ko://' summ_table/eggnog.KO.raw.txt
fi

分箱（binning）

宏基因组binning是指将不同的序列集合（如metagenome序列集合）根据它们的物种归类到不同的bins中，以便进一步研究它们的组成和功能。这个过程可以将类似的序列组合在一起，形成代表不同物种或基因组的bins，以便进行后续分析，如物种注释、基因组组装等。

以下是常用的宏基因组binning方法：

  1. 基于聚类的方法：该方法使用序列聚类将相似序列分到同一个bin中。一般来说，聚类算法可分为两类：无监督聚类（如k-means、DBSCAN等）和有监督聚类（如CAMI、MyCC等）。

  2. 基于组装的方法：该方法使用de novo组装来将相似序列组装成连续的序列，再根据这些序列的基因组信息来将其分类到不同的bins中。这种方法的优点是可以更好地处理重复序列，缺点是需要大量的计算资源和时间。

  3. 基于分类器的方法：该方法使用机器学习分类器来将序列分配到不同的bins中。这种方法的优点是可以自动学习特征并在处理大规模数据时效率高，缺点是需要先建立一个分类器并进行训练。

在进行宏基因组binning时，通常需要使用多个方法进行比较，以选择最适合数据集的方法。可以使用一些流行的工具来进行binning，如MetaBAT、MaxBin、CONCOCT和MEGAN等。这些工具通常包含各种binning方法，可以根据数据集和分析目的选择适合的方法。

篇幅限制，具体的方法放在另一篇里面讲解吧，查看宏基因组分箱流程。

Reference

  1. Y.-X. Liu, Y. Qin, T. Chen, M. Lu, et al., A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell. 12, 315–330 (2021).

关注公众号，获取最新推送

关注公众号 ‘bio llbug’,获取最新推送。

相关文章：

宏基因组分析流程（Metagenomic workflow）202405｜持续更新

Logs 增加R包pctax内的一些帮助上游分析的小脚本&#xff08;2024.03.03&#xff09;增加Mmseqs2用于去冗余&#xff0c;基因聚类的速度非常快&#xff0c;且随序列量线性增长&#xff08;2024.03.12&#xff09;更新全文细节&#xff08;2024.05.29&#xff09; 注意&#x…...

编程日记 2024/5/30 12:14:21

一千题，No.0037（组个最小数）

给定数字 0-9 各若干个。你可以以任意顺序排列这些数字&#xff0c;但必须全部使用。目标是使得最后得到的数尽可能小&#xff08;注意 0 不能做首位&#xff09;。例如&#xff1a;给定两个 0&#xff0c;两个 1&#xff0c;三个 5&#xff0c;一个 8&#xff0c;我们得到的最…...

编程日记 2024/5/30 12:12:20

PV PVC

默写 1 如何将pod创建在指定的Node节点上 node亲和、pod亲和、pod反亲和: 调度策略 匹配标签 操作符 nodeAffinity 主机 In,NotIn,Exists,DoesNotExist&#xff0c;Gt&#xff0c;Lt podAffinity …...

编程日记 2024/5/30 12:11:18

深入理解Nginx配置文件：全面指南

Nginx 是一个高性能的 HTTP 服务器和反向代理服务器&#xff0c;也是一个电子邮件&#xff08;IMAP/POP3&#xff09;代理服务器。由于其高效性和灵活性&#xff0c;Nginx 被广泛应用于各种 web 服务中。本文将详细介绍 Nginx 配置文件的结构和主要配置项&#xff0c;帮助你深入…...

编程日记 2024/5/30 12:09:15

【传知代码】自监督高效图像去噪（论文复现）

前言&#xff1a;在数字化时代&#xff0c;图像已成为我们生活、工作和学习的重要组成部分。然而&#xff0c;随着图像获取方式的多样化&#xff0c;图像质量问题也逐渐凸显出来。噪声&#xff0c;作为影响图像质量的关键因素之一&#xff0c;不仅会降低图像的视觉效果&#xf…...

编程日记 2024/5/30 12:08:14

linnux上安装php zip(ZipArchive)、libzip扩展

安装顺序&#xff1a; 安装zip&#xff08;ZipArchive&#xff09;&#xff0c;需要先安装libzip扩展 安装libzip&#xff0c;需要先安装cmake 按照cmake、libzip、zip的先后顺序安装 下面的命令都是Linux命令 1、安装cmake 确认是否已安装 cmake --version cmake官网 未安装…...

编程日记 2024/5/30 12:07:13

油封制品中各种橡胶材料的差异

在机械系统中&#xff0c;油封起着关键的作用&#xff0c;其主要功能是防止润滑剂泄漏和污染物进入。油封的性能很大程度上取决于所用的橡胶材料。不同的橡胶化合物各有其独特的特性、优点和应用场景。本文将详细探讨油封制品中各种橡胶材料的差异&#xff0c;重点分析其特性、…...

编程日记 2024/5/30 12:06:12

梳理清楚的echarts地图下钻和标点信息组件

效果图 说明 默认数据没有就是全国地图&#xff0c; $bus.off("onresize")是地图容器变化刷新地图适配的&#xff0c;可以你们自己写 getEchartsFontSize是适配字体大小的&#xff0c;getEchartsFontSize(0.12) 12 mapScatter是base64图片就是图上那个标点的底图 Ge…...

编程日记 2024/5/30 12:05:10

【busybox记录】【shell指令】readlink

目录 内容来源&#xff1a; 【GUN】【readlink】指令介绍 【busybox】【readlink】指令介绍 【linux】【readlink】指令介绍 使用示例&#xff1a; 打印符号链接或规范文件名的值 - 默认输出 打印符号链接或规范文件名的值 - 打印规范文件的全路径 打印符号链接或规范文…...

编程日记 2024/5/30 12:04:09

C++之vector

1、标准库的vector类型 2、vector对象的初始化 3、vector常用成员函数 #include <vector> #include <algorithm> #include <iostream> using namespace std;typedef vector<int> INTVEC;// 普通方法 //void showVec(const INTVEC& vec) // 这边如…...

编程日记 2024/5/30 12:03:07

【简单介绍下idm有那些优势】

&#x1f3a5;博主&#xff1a;程序员不想YY啊 &#x1f4ab;CSDN优质创作者&#xff0c;CSDN实力新星&#xff0c;CSDN博客专家 &#x1f917;点赞&#x1f388;收藏⭐再看&#x1f4ab;养成习惯 ✨希望本文对您有所裨益&#xff0c;如有不足之处&#xff0c;欢迎在评论区提出…...

编程日记 2024/5/30 11:57:00

MyBatis系统学习 - 使用Mybatis完成查询单条，多条数据，模糊查询，动态设置表名，获取自增主键

上篇博客我们围绕Mybatis链接数据库进行了相关概述&#xff0c;并对Mybatis的配置文件进行详细的描述&#xff0c;本篇博客也是建立在上篇博客之上进行的&#xff0c;在上面博客搭建的框架基础上&#xff0c;我们对MyBatis实现简单的增删改查操作进行重点概述&#xff0c;在MyB…...

编程日记 2024/5/30 11:55:59

Generative Action Description Prompts for Skeleton-based Action Recognition

标题&#xff1a;基于骨架的动作识别的生成动作描述提示 源文链接&#xff1a;https://openaccess.thecvf.com/content/ICCV2023/papers/Xiang_Generative_Action_Description_Prompts_for_Skeleton-based_Action_Recognition_ICCV_2023_paper.pdfhttps://openaccess.thecvf.c…...

编程日记 2024/5/30 11:53:56

动手学深度学习（Pytorch版）代码实践 -深度学习基础-02线性回归基础版

02线性回归基础版 主要内容 数据生成&#xff1a;使用线性模型 ( y X*w b ) 加上噪声生成人造数据集。数据读取&#xff1a;通过小批量读取数据集来实现批量梯度下降&#xff0c;打乱数据顺序并逐批返回特征和标签。模型参数初始化&#xff1a;随机初始化权重和偏置&#x…...

编程日记 2024/5/30 11:48:50

信息学奥赛初赛天天练-15-阅读程序-深入解析二进制原码、反码、补码，位运算技巧，以及lowbit的神奇应用

更多资源请关注纽扣编程微信公众号 1 2021 CSP-J 阅读程序1 阅读程序&#xff08;程序输入不超过数组或字符串定义的范围&#xff1b;判断题正确填 √&#xff0c;错误填&#xff1b;除特 殊说明外&#xff0c;判断题 1.5 分&#xff0c;选择题 3 分&#xff09; 源码 #in…...

编程日记 2024/5/30 11:47:49

期权具体怎么交易详细的操作流程？

期权就是股票&#xff0c;唯一区别标的物上证指数&#xff0c;会看大盘吧&#xff0c;交易两个方向认购做多&#xff0c;认沽做空&#xff0c;双向t0交易&#xff0c;期权具体交易流程可以理解选择方向多和空&#xff0c;选开仓的合约&#xff0c;买入开仓和平仓没了&#xff0…...

编程日记 2024/5/30 11:46:46

系统架构设计师【第3章】: 信息系统基础知识 (核心总结)

文章目录 3.1　信息系统概述3.1.1　信息系统的定义3.1.2　信息系统的发展3.1.3　信息系统的分类3.1.4　信息系统的生命周期3.1.5　信息系统建设原则3.1.6　信息系统开发方法 3.2　业务处理系统&#xff08;TPS&#xff09;3.2.1　业务处理系统的概念3.2.2　业务处理系统的功能 …...

编程日记 2024/5/30 11:43:43

Linux 驱动设备匹配过程

一、Linux 驱动-总线-设备模型 1、驱动分层 Linux内核需要兼容多个平台&#xff0c;不同平台的寄存器设计不同导致操作方法不同&#xff0c;故内核提出分层思想&#xff0c;抽象出与硬件无关的软件层作为核心层来管理下层驱动&#xff0c;各厂商根据自己的硬件编写驱动…...

编程日记 2024/5/30 11:41:40

游戏子弹类python设计与实现详解

新书上架~&#x1f447;全国包邮奥~ python实用小工具开发教程http://pythontoolsteach.com/3 欢迎关注我&#x1f446;&#xff0c;收藏下次不迷路┗|&#xff40;O′|┛ 嗷~~ 目录 一、引言 二、子弹类设计思路 1. 属性定义 2. 方法设计 三、子弹类实现详解 1. 定义子弹…...

编程日记 2024/5/30 11:40:39

Python基础学习笔记（六）——列表

目录 一、一维列表的介绍和创建二、序列的基本操作1. 索引的查询与返回2. 切片3. 序列加 三、元素的增删改1. 添加元素2. 删除元素3. 更改元素 四、排序五、列表生成式 一、一维列表的介绍和创建 列表&#xff08;list&#xff09;&#xff0c;也称数组&#xff0c;是一种有序、…...

编程日记 2024/5/30 11:38:37

HTML 语义化

目录 HTML 语义化HTML5 新特性HTML 语义化的好处语义化标签的使用场景最佳实践 HTML 语义化 HTML5 新特性 标准答案&#xff1a; 语义化标签&#xff1a; <header>&#xff1a;页头<nav>&#xff1a;导航<main>&#xff1a;主要内容<article>&#x…...

编程新知 2026/1/15 21:24:41

Java 语言特性(面试系列1)

一、面向对象编程 1. 封装&#xff08;Encapsulation&#xff09; 定义&#xff1a;将数据&#xff08;属性&#xff09;和操作数据的方法绑定在一起&#xff0c;通过访问控制符&#xff08;private、protected、public&#xff09;隐藏内部实现细节。示例&#xff1a; public …...

编程新知 2025/10/6 4:03:40

PHP和Node.js哪个更爽?

先说结论&#xff0c;rust完胜。 php&#xff1a;laravel&#xff0c;swoole&#xff0c;webman&#xff0c;最开始在苏宁的时候写了几年php&#xff0c;当时觉得php真的是世界上最好的语言&#xff0c;因为当初活在舒适圈里&#xff0c;不愿意跳出来&#xff0c;就好比当初活在…...

编程新知 2025/9/14 7:59:52

iPhone密码忘记了办？iPhoneUnlocker，iPhone解锁工具Aiseesoft iPhone Unlocker 高级注册版​分享

平时用 iPhone 的时候&#xff0c;难免会碰到解锁的麻烦事。比如密码忘了、人脸识别 / 指纹识别突然不灵&#xff0c;或者买了二手 iPhone 却被原来的 iCloud 账号锁住&#xff0c;这时候就需要靠谱的解锁工具来帮忙了。Aiseesoft iPhone Unlocker 就是专门解决这些问题的软件&…...

编程新知 2026/1/29 10:22:28

鸿蒙中用HarmonyOS SDK应用服务 HarmonyOS5开发一个生活电费的缴纳和查询小程序

一、项目初始化与配置 1. 创建项目 ohpm init harmony/utility-payment-app 2. 配置权限 // module.json5 {"requestPermissions": [{"name": "ohos.permission.INTERNET"},{"name": "ohos.permission.GET_NETWORK_INFO"…...

编程新知 2025/11/26 13:15:57

Mac下Android Studio扫描根目录卡死问题记录

环境信息 操作系统: macOS 15.5 (Apple M2芯片)Android Studio版本: Meerkat Feature Drop | 2024.3.2 Patch 1 (Build #AI-243.26053.27.2432.13536105, 2025年5月22日构建) 问题现象 在项目开发过程中&#xff0c;提示一个依赖外部头文件的cpp源文件需要同步&#xff0c;点…...

编程新知 2026/2/6 1:00:00

以光量子为例，详解量子获取方式

光量子技术获取量子比特可在室温下进行。该方式有望通过与名为硅光子学&#xff08;silicon photonics&#xff09;的光波导&#xff08;optical waveguide&#xff09;芯片制造技术和光纤等光通信技术相结合来实现量子计算机。量子力学中&#xff0c;光既是波又是粒子。光子本…...

编程新知 2026/2/7 11:51:32

A2A JS SDK 完整教程：快速入门指南

目录 什么是 A2A JS SDK?A2A JS 安装与设置A2A JS 核心概念创建你的第一个 A2A JS 代理A2A JS 服务端开发A2A JS 客户端使用A2A JS 高级特性A2A JS 最佳实践A2A JS 故障排除 什么是 A2A JS SDK? A2A JS SDK 是一个专为 JavaScript/TypeScript 开发者设计的强大库&#xff…...

编程新知 2025/12/12 8:42:24

安全突围：重塑内生安全体系：齐向东在2025年BCS大会的演讲

文章目录 前言第一部分&#xff1a;体系力量是突围之钥第一重困境是体系思想落地不畅。第二重困境是大小体系融合瓶颈。第三重困境是“小体系”运营梗阻。 第二部分&#xff1a;体系矛盾是突围之障一是数据孤岛的障碍。二是投入不足的障碍。三是新旧兼容难的障碍。 第三部分&am…...

编程新知 2026/1/25 2:04:42

华为OD机考-机房布局

import java.util.*;public class DemoTest5 {public static void main(String[] args) {Scanner in new Scanner(System.in);// 注意 hasNext 和 hasNextLine 的区别while (in.hasNextLine()) { // 注意 while 处理多个 caseSystem.out.println(solve(in.nextLine()));}}priv…...

编程新知 2026/1/31 6:01:08