二代测序方式分为三种:
- single read 单端测序
- paired-end read 双端测序
- mate-pair read 配对测序
- single read 单端测序
虽然有三种方式,但是大多数数据(包括本课题)为双端测序,至于三者之间的差别、优缺点以及适用场合可以自行搜索。
接下来讲的内容都是基于双端测序的。
测序仪原始下机的数据我们称为raw data,二代测序是将DNA片段打断了再测的,每个测序片段我们称为read,质控完的数据称为clean data。既然是双端测序,那么文件就是成对出现的,分别记录reads两端的信息:一般的命名是*.fq1.gz、*.fq2.gz(’ * ‘ 表示通配符),这是一个fastq文件,通常以fq或fastq作为后缀,具体内容下方会介绍。这里我给大家展示我们课题的一个真实数据,给大家看看大小以及命名方式。
我们可以看到,首先我把我的用户名打码了,呵呵,开个玩笑。
首先这个样本总共有五个文件,这就是公司给的原始数据,我原封不动的上传到了我的服务器上,也没改名也没怎样,我在图片上画了五个蓝色的框,下面是分别代表的意思:
1 | W2018001:是样本的名字,其实我一开始提交给公司的名字是2018001,为什么多了个W我也不知道啊。 |
- 我们可以看到,这个30X的WGS的数据量还是很大的,所以为什么他会是压缩文件,不同公司给的命名可能不一样,不过大体不会相差太多,具体命名的含义也可以问问公司的。
- 这里给到的是一个样本2对文件,实际情况是,可能有多对或者只有一对文件,对于这样的情况,我们的处理方式一般是看作一个文件处理,就是merge一下,有见过有些课题组是在这一步直接merge,我的处理方式是后期bam文件再merge。
- 好了,拿到这么一个数据,我们做的第一件事情应该是校验数据的完整性!!!这点很重要,即使一般情况下,传输都不会出错,但这一步还是要做的,怎么去做呢,就看这个截图里的第一个文件MD5.txt,其实这就是一个文本文件,打开了看看是这样的
- 就是这样的四行信息,也就是一个文件名对应一个校验码,最简单的校验方式,linux系统自带的MD5校验命令:
1 | $ md5sum -c MD5.txt |
- 这个还是需要一点点时间的,大家要耐心等待。展示一下我这里的结果吧(其实在上传完数据的时候我就校验过)
好了,校验完了,没问题,那么就给大家介绍一下文件的格式。
文件是什么呢,其实就是文本文件,可以直接查看,以下是示例:
一条read一条记录,一条记录占四行,第一行是注释,第二行是序列,就是我们所说的ATCG碱基序列,第三行是‘+’,第四行是对应的每个碱基的测序质量,也就是fastq中的“q”。每条记录之间是没有空格的。
贴一下关于第一行的解释
EAS139 | The unique instrument name |
---|---|
136 | Run ID |
FC706VJ | Flowcell ID |
2 | Flowcell lane |
2104 | Tile number within the flowcell lane |
15343 | ‘x’-coordinate of the cluster within the tile |
197393 | ‘y’-coordinate of the cluster within the tile |
1 | Member of a pair, 1 or 2 (paired-end or mate-pair reads only) |
Y | Y if the read fails filter (read is bad), N otherwise |
18 | 0 when none of the control bits are on, otherwise it is an even number |
ATCACG | Index sequence |
第一个是测序机器ID,独一无二的,你可以去illumina官网查到的(我没查过)
第二个是这个机器跑的次数,一般机器都是有使用寿命的,所以次数越多肯定结果越差,那么一般在什么区间合适呢,我的老师给的建议是200-9999,为啥还有下限呢,那是因为,机器刚开始跑的那几次,需要磨合嘛,不是很准。看来这个展示数据不是很好啊,吓得我赶紧去看了下我的数据,还好,在212次(懒得贴图了,hexo弄个图比简书还麻烦),在范围内。
倒数第三个,表示数据有没有被过滤过,一般我们的原始数据肯定是N的,要是给的原始数据是Y,你就得好好问问公司原因了。
其余具体的我就不解释了,有问题可以评论交流。
第二行要注意的是,可能有N的出现,那是因为有些碱基没被测出来。
第四行是ASCII码表示的碱基质量,这样就能保证质量是用一个字符表示的,和碱基一一对应。公式如下:
P=10^-[(n-33)/10]举个栗子:比如“?”在ASCII码表上对应的编号是63,那么n就是63,减去33以后就是30,也就是我们说的Q30了,所以
Q = n - 33
,P算出来就是0.001,这个就是错误率,反过来,准确率就是99.9%,Q30就是准确率99.9%,同理Q20就是准确率99%。ok,介绍完格式,介绍两款质控的软件,fastqc和fastp
1 | fastqc: |
- 以上是我的安装方式,fastqc是一个很常用的软件,我想只要不是小白都用过。对于ubuntu用户请用我提供的方式进行安装,用
apt-get install fastqc
可能在使用的时候出现如下报错
- 至于fastp,是用来清洗质量不好的reads的,其实类似的软件很多,包括trimmomatic等,之所以选择这个是因为用过那么多软件后,发现fastp无论在安装还是使用上都很顺手,仅此而已。
- 好了,现在让我们来使用fastqc对raw data的质量进行分析并查看
1 | $ fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN |
- 结果如下,要放在早几年,估计是很难见到这么干净的原始数据了,为此我还专门请教了一些老师,这么干净的数据的可信程度。得益于illumina公司对仪器的升级改善,数据的质量也是越来越好了。
但是即便如此,我们还是要对raw data进行清洗的,数据的质量直接决定了后面分析的准确性,是一切分析的前提,我们一直都在说QC,QC也是基础,但是它的重要性不容忽视。下面是我要查看的质控的内容:
• read各个位置的碱基质量值分布
• 碱基的总体质量值分布
• read各个位置上碱基分布比例,目的是为了分析碱基的分离程度
• GC含量分布 • read各位置的N含量
• read是否还包含测序的接头序列下面是我的对于clean data的指标:
1.将含有接头的reads对去除,或者cut接头
2.测序错误率分布检查,一般情况下,每个碱基位置的测序错误率都应该低于1%
3.GC含量,理论上A=T,C=G,前几个碱基可能会有波动,GC含量整体在41.5%左右
4.平均质量值分布:Q30平均比例≥80%,平均错误率≤0.1%为了达到这样的目标,我的筛选是条件是:
去除含接头(adapter)的一对reads;可以酌情考虑去除接头保留reads
当某一端read中的N含量超过该条read长度比例10%,去除该对reads
当某一端read中低质量(Q≤20)碱基数超过该read长度比例的50%时,去除该对reads
1 | $ fastp -i 1.fq.gz -I 2.fq.gz \ |
别的不多说,解释一下我用到的几个参数:
-c:对overlap区域进行纠错,适用于PE
-q:设置低质量的标准,默认是15,多数公司这里设为5
-u:低质量碱基所占比例,默认40代表40%,只要有一条read不满足条件就成对丢掉
-n:过滤N碱基过多的reads,15代表个数,因为一般PE150的reads长度是150
-w:线程数,这里要注意!范围是1-16,建议是8,亲测,8和16差不多
-5 -3:根据质量值来截取reads,分别对应5‘端和3’端,得到reads长度可能不等具体请参考官网说明
这个时候再看clean data结果,就无需fastqc了,因为fastp也会生成一份报告。
- 看到结果,其实可以发现,数据还是保留了大部分的,所以完全不用担心cut太狠,导致数据不够。我的原则是,在数据量面前,更应该保证的是数据的准确性,即使会有一点损失,也是值得的,毕竟research更注重的也是准确性,相对于降低假阴性,我更倾向于提高真阳性。
- 那么得到clean data以后就可以进行下一步mapping了,mapping内容下次再与大家分享。
- 水平有限,要是存在什么错误请指出,可发送邮箱至shiyuant@outlook.com!请大家多多批评指正,相互交流,共同成长,谢谢!!!