为进一步提高《微生物组实验手册》稿件质量，本项目新增大众评审环节。文章在通过同行评审后，采用公众号推送方式分享全文，任何人均可在线提交修改意见。公众号格式显示略有问题，建议点击文末阅读原文下载PDF审稿。在线文档(https://kdocs.cn/l/cL8RRqHIL)大众评审页面登记姓名、单位和行号索引的修改建议。修改意见的征集截止时间为推文发布后的72小时，文章将会结合有建设性的修改意见进一步修改后获得DOI在线发表，同时将致谢审稿人。感谢广大同行提出宝贵意见。

随机宏基因组测序数据质量控制和去宿主的分析流程和常见问题

Analysis pipeline and frequently asked questions of quality control and host removal in shotgun metagenomic sequencing

刘永鑫1, 2, 3, #, *，刘芳1, 2, 3, #，陈同4，白洋1, 2, 3, 5, *

1中国科学院遗传与发育生物学研究所，植物基因组学国家重点实验室，北京；2中国科学院大学，生物互作卓越创新中心，北京；3中国科学院遗传与发育生物学研究所，中国科学院–英国约翰英纳斯中心植物和微生物科学联合研究中心，北京；4中国中医科学院，中药资源中心，北京；5中国科学院大学现代农学院，北京

*通讯作者邮箱: [email protected] ; [email protected]

#共同第一作者/同等贡献

摘要: 随机宏基因组测序，也称鸟枪法宏基因组测序，是指对环境样品的总DNA进行高通量测序以获得微生物群落的物种组成及其潜在功能，抑或通过序列拼接和分箱得到其微生物的基因组。宏基因组测序数据预处理包括两方面：一方面，与转录组、基因组测序等分析相似的数据质量控制过程，包括质量评估，去除低质量、引物和接头序列；另一方面，涉及到宿主相关微生物的宏基因组样本易受宿主序列的污染，需要去除宿主序列并评估宿主比例，以获得高质量的微生物组相关数据以方便开展下游分析。本文主要介绍FastQC、MultiQC、KneadData（涵盖并调用Trimmomatic + Bowtie 2）等软件组合分析流程的安装、使用方法和结果解读，实现数据质量评估、质量控制和去宿主污染、质量再评估的分析过程，同时对各步骤常见问题和解决方法进行总结，方便同行更准确、高效地实现宏基因组数据的预处理，为下游分析提供高质量的宏基因组数据。

关键词: 宏基因组测序，质量控制，去宿主，FastQC，KneadData

仪器设备

1. 计算服务器（操作系统：Linux主流发行版本，如CentOS 7+ / Ubuntu 16.04+；CPU：8核+；内存：32G+；硬盘：> 30 GB，且大于原始数据大小3倍），网络访问畅通。

2. 个人电脑（Windows用户需安装XShell或Putty等终端类软件，Mac使用系统内置终端）即可远程访问计算服务器。

软件和数据库

1. 远程文件传输工具FileZilla客户端3.49.1+：https://filezilla-project.org/

2. （可选）Windows远程访问服务器终端工具Xshell 6.0.0197p+：https://www.netsarang.com/zh/free-for-home-school/

3. 软件管理器Miniconda2 Linux 64-bit （Python 2.7）: https://conda.io/miniconda.html

4. 测序数据质量评估FastQC v0.11.9：https://www.bioinformatics.babraham.ac.uk/projects/download.html

5. 质量评估报告汇总MultiQC version 1.6 （Ewels等，2016）：https://multiqc.info/

6. 宏基因组质量控制和去宿主分析流程KneadData v0.7.4: http://huttenhower.sph.harvard.edu/kneaddata

7. （可选）并行任务队列管理Parallel 20200522 （Tange，2020）：https://www.gnu.org/software/parallel/     

8. 常用宿主基因组下载Ensembl Genome：http://ensemblgenomes.org/ ，如人类基因组（International Human Genome Sequencing，2001），拟南芥基因组（The Arabidopsis Genome，2000）。

9. 流程参考代码详见：https://github.com/YongxinLiu/MicrobiomeProtocol/blob/master/e1.KneadData/QualityControl_HostRemoval_Pipelie.sh

软件安装和数据库部署

Windows/Mac用户安装FileZilla客户端，用于上传测序数据至服务器或数据中心，也可下载分析结果本地查看。Windows用户安装Xshell用于远程访问服务器并开展分析，Mac用户可使用系统自带Terminal中的ssh命令远程访问服务器。

在Linux系统的计算服务器端，以Miniconda2软件和Python2虚拟环境安装所需软件，在将来随着软件的更新可能需要新建Python3虚拟环境才能安装新版本；然后下载人类基因组索引，同时以拟南芥为例介绍下载基因组并建立索引的步骤。

注：代码行添加灰色底纹背景，其中需要根据系统环境修改的部分标为蓝色。

1. 安装Miniconda2 Linux 64-bit（Python 2.7），已经安装Conda可跳过此步骤。

wget -c https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh

bash Miniconda2-latest-Linux-x86_64.sh

2. 配置Conda环境，添加Bioconda生物频道以方便安装生物学相关的分析软件。

conda config --add channels bioconda

conda config --add channels conda-forge

3. Conda新建Python 2.7环境，命名为qc2（quality control python2），然后进入。

conda create -n qc2 python=2.7

conda activate qc2

注：新建虚拟环境，然后在新建的环境下安装工作流程，可以防止新装的软件或者其依赖软件与系统默认环境中的版本相互冲突。另外，将整个分析流程的软件存放在虚拟环境并放置在指定目录下，不用时可以轻松移除，不会对系统产生任何影响。

4. Conda安装相关软件，-y默认同意直接安装，不再提示是否确认。

conda install fastqc -y

conda install multiqc -y

conda install kneaddata -y

conda install parallel -y

注：如果软件下载慢或无法下载，详见常见问题1。Conda默认安装Bioconda中的最新版本或所处系统环境支持的最新版本；如果无法安装或安装后使用存在问题，可使用conda remove xxx移除某软件，再指定版本安装，如指定安装KneadData的0.6.1版本：conda install kneaddata=0.6.1。

5. 宿主基因组数据库下载。

为了方便指定接下来的文件路径，我们首先使用mkdir命令为整个分析流程建立一个文件夹，并命名为meta_preprocess（参数-p允许建立多级文件夹、多个文件夹且不报错）。然后使用cd命令进入该文件夹。

mkdir -p meta_preprocess

cd meta_preprocess

为了去除宿主序列，我们需要建立宿主序列的索引以供KneadData通过序列比对找到并去除宿主序列。KneadData提供了多个预先建立的常用的宿主序列索引。下面的命令可供我们查看KneadData软件整理好的可用的数据库索引，包括人类基因组、小鼠基因组、人类转录组和核糖体数据库等。

            kneaddata_database

以人类基因组为例，下载Bowtie 2格式索引，此类索引文件通过包含多个文件，推荐建立文件夹并指定下载位置。

mkdir -p db

kneaddata_database --download human_genome bowtie2 db/

如果默认数据库下载速度慢或无法下载，可使用国内备份链接，详见常见问题2。

KneadData包括的数据库种类有限，用户可自行下载参考基因组并建索引，以拟芥为例的实例详见常见问题3。

6. 准备输入数据

通常测序公司会返回原始（raw）或纯净（clean）数据两类数据：原始数据为下机后按测序文库的索引（Index）拆分获得的样本序列，纯净数据是去除了明显的低质量、测序引物和接头污染序列后的结果。推荐大家使用体积更小、质量更高的纯净序列进行下游分析和提交数据中心。此外，涉及人类研究的数据，需要上传去除人类相关序列后再上传数据中心（即本文的输出结果）。

本文使用的数据来自人类口腔癌症研究的文章（Schmidt等，2014），NCBI的SRA项目号为PRJEB4953。为方便演示流程的使用，我们从中选取4个样本，并且随机抽取了75000对序列作为软件的测序数据，可以从中国科学院基因组研究所的原始数据归档库（Genome Sequence Archive，GSA，https://bigd.big.ac.cn/gsa/ ）（Wang等，2017）中按批次编号CRA002355搜索并下载，也可通过wget并结合for循环通过批次和样本编号批量下载至seq目录（代码如下）。

    mkdir -p seq

使用wget下载单个样本，-c为支持断点续传，-O指定保存位置并可重命名，每个双端样本需要下载两个文件。

wget -c ftp://download.big.ac.cn/gsa/CRA002355/CRR117732/CRR117732_ f1.fq.gz -O seq/C2_1.fq.gz

wget -c ftp://download.big.ac.cn/gsa/CRA002355/CRR117732/CRR117732_ r2.fq.gz -O seq/C2_2.fq.gz

结合for循环再下载3个样本，seq命令产生连续序列，$i替换命令中可变部分，结尾加保证变量名结束而被识别。

    for i in `seq 3 5`;do

 wget -c ftp://download.big.ac.cn/gsa/CRA002355/CRR11773$i/CRR11773$i_f1.fq.gz

        -O seq/C$i_1.fq.gz

 wget -c ftp://download.big.ac.cn/gsa/CRA002355/CRR11773$i/CRR11773$i_r2.fq.gz

       -O seq/C$i_2.fq.gz

done

视频1. 宏基因组测序数据分析流程演示视频和讲解

（https://v.qq.com/x/page/a3128efr2t3.html ）

实验步骤

开始分析前，我们应处于项目所在目录（如meta_preprocess），并启动软件所在的Conda环境。

cd meta_preprocess

conda activate qc2

1. FastQC测序数据质量评估。

fastqc seq/*.fq.gz -t 3

*.fq.gz代表所有以.fq.gz结尾的文件，即所有测序数据；-t 3指定3个线程，即同时对3个文件进行并行分析。

图1. FastQC质量评估报告中的主要结果和注意事项。A. 序列中每个碱基的质量分布（Per base sequence quality）。B. 所有序列的GC含量（Per sequence GC content）分布（红色）与理论值分布（蓝色）曲线。C. 接头含量（Adapter content）。本图数据为样本C3右端序列为列对fastQC的评估结果进行说明，完整评估报告详见seq/C3_2_fastqc.html。

FastQC质量评估包括基本统计（比如对应样本总序列数，序列长度和GC含量等简要总结）、单碱基位点测序质量、GC含量及接头含量等10大类的评估。我们以C3样本右端报告为例，首先查看基本统计中的总序列数（Total Sequences）和GC含量（%GC）等。其次查看每个碱基位点的质量分数的箱线图（图1A），每个箱体中间的红线代表此位置上所有序列的测序质量的中位数，然后黄色箱体代表25%-75%百分位数内的质量分布，而两端黑线顶端对应10%和90%百分位的质量数，另外连接每个箱体的蓝色线代表的是平均值。根据Y轴序列质量，整个图片区域被划分为高（绿色，得分>=28）、中（黄色，