宏基因组分析实战教程2. 数据质控fastqc, Trimmomatic, MultiQC

最新推荐文章于 2024-04-20 17:10:38 发布

刘永鑫Adam

最新推荐文章于 2024-04-20 17:10:38 发布

阅读量1.3w

点赞数 6

分类专栏：宏基因组宏基因组分析

宏基因组同时被 2 个专栏收录

76 篇文章 167 订阅

订阅专栏

宏基因组分析

18 篇文章 317 订阅

订阅专栏

本文英文原版见下方github链接，由中科院朱微金博士翻译、测试、并进行中文注释和补充，全网首发“宏基因组”公众号。

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

前情提要

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

测试数据

刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：
1. 下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
2. 预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
3. 数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
4. 方便结果预览和跳过问题步骤；服务器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 = 8888

EOF

# 
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.dist

abundance-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

结果只经过了简单的尾部过滤，k-mer的数量减少了10%以上，对下游分析的准确度和速度都非常有帮助。

按Kmer质控后的结果，感觉趣的再用fastqc评估一下，看看有什么变化？

接下来的文章来会介绍k-mer更大的用途，猜猜是什么？

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内五十位PI，六百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加创始人好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决推荐生信技能树-微生物组版块(http://www.biotrainee.com/forum-88-1.html) 发贴，并转发链接入群，问题及解答方便检索，造福后人。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA

刘永鑫Adam

关注

6
点赞
踩
63

收藏

觉得还不错? 一键收藏
1
评论
宏基因组分析实战教程2. 数据质控fastqc, Trimmomatic, MultiQC

本文英文原版见下方github链接，由中科院朱微金博士翻译、测试、并进行中文注释和补充，全网首发“宏基因组”公众号。https://2017-cicese-metagenomics.readthedocs.io/en/latest/toc.html前情提要如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章宏基因组分析理论教程微生物组入门圣经+宏基因组分析实操课程1.
复制链接

扫一扫