宏基因组实战2. 数据质控fastqc, Trimmomatic, MultiQC, khmer

前情提要

如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章

宏基因组分析理论教程
微生物组入门圣经+宏基因组分析实操课程
1. 背景知识-Shell入门与本地blast实战

测试数据

刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：

下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
方便结果预览和跳过问题步骤；服务器Linux在不同平台和版本下，软件安装和兼容性问题还是很多的，而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助；想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中，有每一步的运行结果，读者可以下载查看分析结果，并可基于此结果进一步分析。不要纠结于某一步无法通过，重点是了解整个流程的分析思路。

最后送上本教程使用到的所有文件同步共享文件夹链接：http://pan.baidu.com/s/1hsIjosk 密码：y0tb 。

数据质控

https://2017-cicese-metagenomics.readthedocs.io/en/latest/quality.html # 有时连接不稳定打不开，等会就会好。或访问它更早版本的链接如下：

https://2017-dibsi-metagenomics.readthedocs.io/en/latest/quality.html

安装软件

安装依赖关系

sudo apt-get -y update && \
sudo apt-get -y install trimmomatic python-pip \samtools zlib1g-dev ncurses-dev python-dev unzip \python3.5-dev python3.5-venv make \libc6-dev g++ zlib1g-de

安装 fastqc

wget -c http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
unzip fastqc_v0.11.5.zip
cd FastQC
chmod +x fastqc
cd

创建Python3.5虚拟环境

cd
python3.5 -m venv ~/py3
. ~/py3/bin/activate
pip install -U pip
pip install -U Cython
pip install -U jupyter jupyter_client ipython pandas matplotlib scipy scikit-learn khmer
pip install -U https://github.com/dib-lab/sourmash/archive/master.zip

运行Jupyter Notebook

# 配置
jupyter notebook --generate-config -y
cat >>~/.jupyter/jupyter_notebook_config.py <<EOF
c = get_config()
c.NotebookApp.ip = '*'
c.NotebookApp.open_browser = False
c.NotebookApp.password = u'sha1:5d813e5d59a7:b4e430cf6dbd1aad04838c6e9cf684f4d76e245c'
c.NotebookApp.port = 8888EOF# 
jupyter notebook &

1. 测序数据准备

我们分析采用 Hu et al., 2016. 文章中数据的子集，下载数据

# 创建数据文件夹
mkdir data
cd data
# 下载测试数据
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_1.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_2.fastq.gz
# 如果无法科学上网而下载失败，尝试在文提供的百度云中的data目录中下载# 检查文件
md5sum *.gz
# 改原始文件为只读，防止被修改
chmod u-w *

2. fastqc质量评估

# 质控所有gz压缩的原始数据，t启动多线程，一般与文件数量一致
fastqc *.gz -t 4
# 显示所有网页版质量评估报告文件，可下载本地或用firefox查看
ll *.html

3. Trimmomatic去接头和低质量序列

下载Illumina双端接头序列

curl -O -L http://dib-training.ucdavis.edu.s3.amazonaws.com/mRNAseq-semi-2015-03-04/TruSeq2-PE.fa

运行Trimmomatics

# 调用for循环批处理文件
for filename in *_1.fastq.gz
do# 提取双端公共文件名，并输出检验
base=$(basename $filename _1.fastq.gz)
echo $base# 运行去接头程序
TrimmomaticPE -threads 9 \${base}_1.fastq.gz \${base}_2.fastq.gz \${base}_1.qc.fq.gz ${base}_s1_se \${base}_2.qc.fq.gz ${base}_s2_se \ILLUMINACLIP:TruSeq2-PE.fa:2:40:15 \LEADING:2 TRAILING:2 \SLIDINGWINDOW:4:2 \MINLEN:25 
done

宏基因组拼接前必须去干净接头，防止引入人造序列对结果影响

4. 质控后再评估

fastqc *.qc.fq.gz -t 4
# 查看再次质控结果，与之前的比较试试
ll *.qc_fastqc.html

图1. 比较质控前后第一个样品右端接头污染水平。上图质控前接头污染水平近10%，质控后接近0.

评估报告的结果非常多，自己多读读，不懂上fastqc官网看帮助。

5. MultiQC多样品报告汇总(可选)

需要python3.5

# 激活Pythone3环境
. ~/py3/bin/activate
# 安装包
pip install git+https://github.com/ewels/MultiQC.git
# 生成多样品报告
multiqc . #

虽然是可选步骤，但对于多样品还是非常有意义的。可以方便比较，节省时间。

图2. 多样品质控前后比较。图像还是交互式的，鼠标悬停可显示样品名。

6. K-mer过滤

https://2017-cicese-metagenomics.readthedocs.io/en/latest/kmer_trimming.html

如果我们绘制样品k-mer丰度的柱状图，你会注意到存在大量的unqiue K-mers，即使测序质量很高，但它们也是由测序错误导致的。

图3. 序列末端低质量区有极高复杂度的kmer

本节继续在Python3下运行

# 对质控前后的数据统计单端丰度距离
abundance-dist-single.py -M 1e9 -k 21 SRR1976948_1.fastq.gz SRR1976948_1.fastq.gz.distabundance-dist-single.py -M 1e9 -k 21 SRR1976948_1.qc.fq.gz SRR1976948_1.qc.fq.gz.dist# 只对高覆盖度中的低丰度kmer剪切(更可能是测序错误)；低覆盖度保留
interleave-reads.py SRR1976948_1.qc.fq.gz SRR1976948_2.qc.fq.gz | trim-low-abund.py -V -M 8e9 -C 3 -Z 10 - -o SRR1976948.trim.fq

图4. kmer过滤原理：
只对高覆盖度中的低丰度kmer剪切(更可能是测序错误)；低覆盖度保

为什么要进行k-mer剪切

如果不做这步也是可以的。但会增加下游组装的工作量，本步可使结果更准确，并增加下游拼接速度，以及内存消耗。

unique-kmers.py SRR1976948_1.qc.fq.gz SRR1976948_2.qc.fq.gz
unique-kmers.py SRR1976948.trim.fq

结果如下：

# 质控后的32-mers数据
Estimated number of unique 32-mers in SRR1976948_1.qc.fq.gz: 65344914
Estimated number of unique 32-mers in SRR1976948_2.qc.fq.gz: 85395776
Total estimated number of unique 32-mers: 112758982# k-mer剪切后的数据
Estimated number of unique 32-mers in SRR1976948.trim.fq: 101285633
Total estimated number of unique 32-mers: 101285633