20160420-序列比对前的准备工作

20160420-序列比对前的准备工作

Hello,大家好!1周不见,甚是想念啊!

在之前的文章中,我们已经为大家介绍了高通量测序技术的原理,储存格式,使用FastQC做质量控制的方法。那么在使用FastQC之后,如果我们发现了一些问题(序列质量不高,),那么我们该使用什么样的工具,去解决这些问题呢?这就是今天我们要干的事情。

测试数据:

为了方便大家学习高通量测序分析,我为大家准备了一个测试数据,是双端测序的fastq文件,每个文件包含50W条reads。以后我们所有教程中展示的结果全部都来自于这1组测试数据。下载地址如下。

fastx Toolkit

fastx Toolkit是包含处理fastq/fasta文件的一系列的工具,它是基于java开发的,我们高通量测序最常用到的是使用这个软件进行reads的裁剪(trim)。它的安装可以根据官网的指南进行:

在这里我们就不再赘述。下面我们对这个工具箱中的工具,具体的命令进行一个简要的介绍。

FASTQ-to-FASTA
说明:这个命令主要是用来转换FASTA格式与FASTQ格式$ fastq_to_fasta [-h] [-r] [-n] [-v] [-z] [-i] [-o]
[-h] = 获得帮助信息.
[-r] = 使用序号去代替fastq文件中原来的reads名.
[-n] = 如果fastq中有N,保留.(默认是有N的序列删除)
[-v] = 报告reads的总数
[-z] = 调用GZip软件,输出的文件自动经过压缩.
[-i] = 输入文件,可以是fastq/fasta格式.
[-o] =输出路径,如果不设置会直接输出到屏幕.

FASTX Statistics
说明:主要统计一下序列的基本信息,如GC含量什么的,很少用,基本使用FastQC代替
$ fastx_quality_stats 
[-h] [-i INFILE] [-o OUTFILE]
[-h] = 获得帮助信息
[-i] = FASTA/Q格式的输入文件
[-o] = 输出路径,如果不设置会直接输出到屏幕.

FASTA/Q Clipper
说明:主要是进行reads的过滤和adapter的裁剪
$ fastx_clipper [-h] [-a ADAPTER] [-i INFILE] [-o OUTFILE]
#参数很多,我们慢慢来看!
[-h] = 获得帮助信息.
[-a ADAPTER] = Adapter序列信息. 默认的是CCTTAAGG
[-l N] = 如果1条reads小于N就抛弃,默认5.
[-d N] = 保留adapter并保留后面的Nbp,如果设置-d 0等于没有用这个参数.
[-c] = 只保留包含adapter的序列
[-C] = 只保留不包含adapter的序列
[-k] = 报告adater的序列信息
[-n] = 如果reads中有N,保留reads.(默认是有N的序列删除)
[-v] = 报告序列总数
[-z] = 调用GZip软件,输出的文件自动经过压缩.
[-D]= Debug output.

FASTA/Q Trimmer
说明:这个是我最常用的工具,可以快速切序列
$ fastx_trimmer [-h] [-f N] [-l N] [-z] [-v] [-i INFILE] [-o OUTFILE]
[-h] = 获得帮助信息.
[-f N] = 序列中从第几个碱基开始保留. 默认是1.
[-l N] = 序列最后保留到多少个碱基,默认是整条序列全部保留.
[-z] = 调用GZip软件,输出的文件自动经过压缩.

备注:其实Fastx工具箱中还有很多其他的工具,但是都不是特别常用,所以我也就不进行介绍了。如果有感兴趣的同学,请直接到他们的官网上去阅读说明文档。

cutadapt软件

这个cutadapt软件是最常用的去adapter的工具。它是基于Python编写的一个Python包,Python包安装的办法在我们的公众号里面已经有比较详细的介绍,我们就不赘述。这里为大家附上cutadapt程序包的下载地址。


成功安装cutadapt以后,在命令行模式下输入:cutadapt就会出现类似于下面的提示,里面主要是cutadapt的基本介绍信息和最简单的用法。

一些最基本的用法
# cutadapt的功能特别强大,相对应的参数真的特别说,有几十个参数,我们平时只会用到很少的几个,我在这里为大家介绍一下。

# 最基本的形式,可以去掉3‘端的adapter序列
$ cutadapt -a AACCGGTT -o output.fastq input.fastq

# 可以直接输入或者输出压缩文件,不需要修改参数,输出文件的后面加上.gz
$ cutadapt -a AACCGGTT -o output.fastq.gz input.fastq.gz

# 假如去掉3‘端的adapter AAAAAAA 和5’端的adapter TTTTTTT
$ cutadapt -a AAAAAAA -g TTTTTTT -o output.fastq input.fastq

# cutadapt也可以用来进行reads的cut
# 去掉最前面的5bp
$ cutadapt -u 5 -o trimmed.fastq input_reads.fastq

# 进行reads测序质量的过滤
# cutadapt软件可以使用-q参数进行reads质量的过滤。基本原理就是,一般reads头和尾会因为测序仪状态或者是反应时间的问题造成测序质量差,比较粗略的一个过滤办法就是-q进行过滤。需要特别说明的是,这里的-q对应的数字和phred值是不一样的,它是软件根据一定的算法计算出来的。
# 3‘端进行一个简单的过滤,--quality-base=33是指序列使用的是phred33计分系统
$ cutadapt -q 10 --quality-base=33 -o output.fastq input.fastq 
# 3‘端 5’端都进行过滤,3'的阈值是10,5‘的阈值是15
$ cutadapt -q 10,15 --quality-base=33 -o output.fastq input.fastq 


Reads 长度的过滤
[--minimum-length N or -m N]
# 当序列长度小于N的时候,reads扔掉

[--too-short-output FILE]
# 上面参数获得的这些序列不是直接扔掉,而是输出到一个文件中

[--maximum-length N or -M N]
# 当序列长度大于N的时候,reads扔掉

[--too-long-output FILE]
# 上面参数获得的这些序列不是直接扔掉,而是输出到一个文件中


Paired-Reads的裁剪(trim)
# 现在很多的测序都是双端测序,那么从测序原理上来说,一对reads来自于1簇反应,所以一起进行adapter的trim可能效果更好。cutadapt自然也提供了这样的功能。
$ cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq
# -a是第1个文件的adapter序列
# -A是第2个文件的adapter序列
# -o是第1个输出文件
# -p是第2个输出文件

1个例子

其实在实际中,我们从公司拿到的数据大多数已经进行过cutadapt,我们其实更关注的是reads的trim。我们利用我们的测试数据进行一下实际处理说明。

我们首先先要用fastqc对test1.fastq与test2.fastq进行一下质量评估,评估的主要结果如下:

read1的测序质量图


read2的测序质量图


我们从上面两张图可以明显看出,read1的测序质量明显好于read2,一般我们确定要trim多少bp的时候都是按phred20这个标准进行评估。比如,对于我们测试数据,read1就不需要trim,read2需要保留1-85bp。对应的fastx_trimmer的命令如下:

fastx_trimmer -i test_data_2.fastq -o test_data_2_trim.fastq -f 1 -l 85

好啦!本期的内容就介绍到这里!

另外欢迎大家关注我们的微信公众号:高通量测序技术


——————————————————————————

另外欢迎各位参加我们的知乎Live:

1. 知乎Live:如何快速入门生物信息学 (涉及内容:测序原理,生物信息学发展历史,软件的安装与调试,入门路线图,介绍了RNA-Seq的分析流程并给出实践代码);

如何快速入门生物信息学


2. 知乎Live: 生信进阶第1课-重复Nature文章 (涉及内容:肺癌相关研究现状,RNA-Seq单细胞测序,RNA-Seq的建库方法,RNA-Seq的分析流程细节,相关生信图的绘制);

生信进阶第1课-重复Nature文章


3. 知乎Live:生信进阶第2课-基因组序列

(涉及内容:介绍基因组的序列结构,hg19与hg38的区别,ENCODE计划,常用的表观组学实验原理ChIP-Seq,Hi-C等,ChIP-Seq的标准处理流程,绘图原理)

生信进阶第2课-基因组序列


4. 知乎Live:不用编程怎么做生物信息学

(涉及内容:介绍生物信息学入门的几个层次,从命令行到图形界面再到命令行,绘制生物进化树,图形界面分析平台,使用图形界面处理RNA-Seq数据,使用图形界面分析ChIP-Seq数据,UCSC genome browser,WashU genome browser)

不用编程怎么做生物信息学

编辑于 2017-09-22

文章被以下专栏收录

    主要介绍以下内容: 传统生信技术,如序列比对,进化树构建,基因结构预测等; 常用数据库教程,如PubMed,NCBI中的序列信息,Ensembl,UCSC Genome Browser,Uniport 等等; 高通量测序技术,如测序方法,比对算法,软件用法,分析流程详解,数据格式等等; Python编程,对于常见生信数据的处理方法; R可视化,将分析好的数据进行可视化,前期主要介绍R的基本用法,后期结合实际例子进行常用的R可视化专题,如手把手教你写自己的heatmap函数分析差异基因等等; Deep Learning部分,介绍Deep Learning在生信中的应用; 癌症基因组学部分。