仪器设备
- 计算服务器(操作系统:Linux主流发行版本,如CentOS 7+/Ubuntu 16.04+;CPU:8核+;内存:32G+;硬盘:> 30 GB,且大于原始数据大小3倍),网络访问畅通
- 个人电脑(Windows用户需安装XShell、Git bash或Putty等终端类软件,Mac使用系统内置终端)即可远程访问计算服务器
软件和数据库
- 远程文件传输工具FileZilla客户端3.49.1+:https://filezilla-project.org/
- (可选) Windows远程访问服务器终端工具Xshell 6.0.0197p+:https://www.netsarang.com/zh/free-for-home-school/
- 软件管理器Miniconda2 Linux 64-bit (Python 2.7): https://conda.io/miniconda.html
- 测序数据质量评估FastQC v0.11.9:https://www.bioinformatics.babraham.ac.uk/projects/download.html
- 质量评估报告汇总MultiQC version 1.6 (Ewels等,2016):https://multiqc.info/
- 宏基因组质量控制和去宿主分析流程KneadData v0.7.4: http://huttenhower.sph.harvard.edu/kneaddata
- (可选)并行任务队列管理Parallel 20200522 (Tange,2020):https://www.gnu.org/software/parallel/
- 常用宿主基因组下载Ensembl Genome:http://ensemblgenomes.org/,如人类基因组(International Human Genome Sequencing,2001),拟南芥基因组(The Arabidopsis Genome,2000)。
- 流程参考代码和示例结果详见:https://github.com/YongxinLiu/MicrobiomeProtocol/blob/master/e1.KneadData/,包括本文的参考代码文档 (*.sh)、质量评估报告 (*.html)和结果统计表 (*.txt)。
软件安装和数据库部署
Windows/Mac用户安装FileZilla客户端,用于上传测序数据至服务器或数据中心,也可下载分析结果本地查看。Windows用户安装Xshell用于远程访问服务器并开展分析,Mac用户可使用系统自带Terminal中的ssh命令远程访问服务器。
在Linux系统的计算服务器端,以Miniconda2软件和Python2虚拟环境安装所需软件,在将来随着软件的更新可能需要新建Python3虚拟环境才能安装新版本;然后下载人类基因组索引,同时以拟南芥为例介绍下载基因组并建立索引的步骤。
注:代码行添加灰色底纹背景,其中需要根据系统环境修改的部分标为蓝色。
- 安装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 - 配置Conda环境,添加Bioconda生物频道以方便安装生物学相关的分析软件。
conda config --add channels bioconda
conda config --add channels conda-forge - Conda新建Python 2.7环境,命名为qc2 (quality control python2),然后进入。
conda create -n qc2 python=2.7
conda activate qc2
注:新建虚拟环境,然后在新建的环境下安装工作流程,可以防止新装的软件或者其依赖软件与系统默认环境中的版本相互冲突。另外,将整个分析流程的软件存放在虚拟环境并放置在指定目录下,不用时可以轻松移除,不会对系统产生任何影响。 - 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。 - 宿主基因组数据库下载。
为了方便指定接下来的文件路径,我们首先使用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。 - 准备输入数据
通常测序公司会返回原始(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 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
使用wget下载单个样本,-c为支持断点续传,-O指定保存位置并可重命名,每个双端样本需要下载两个文件。
结合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. 宏基因组测序数据分析流程演示视频和讲解
实验步骤