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

本文英文原版见下方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

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

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

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

需要python3.5

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

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

image
image

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

6. K-mer过滤

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

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

image
图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

image
图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) 发贴,并转发链接入群,问题及解答方便检索,造福后人。
image

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

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

发布了542 篇原创文章 · 获赞 346 · 访问量 137万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 编程工作室 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览