视频地址:https://www.youtube.com/watch?v=dgzkXLfFAf8&list=PLjiXAZO27elDHGlQwfd06r7coiFtpPkvI&index=3
这一个视频主要是讲数据预处理的过程。主要分成三步:
(1)比对
(2)标记duplicates
(3)碱基质量分数重新校准
我们为什么要进行数据的预处理呢?(用主讲人的话就是:garbage-in和garbage-out的过程)我们拿到的数据受到技术偏差和人为因素的影响,会产生一些duplicates,所以在做突变体calling之前你要对你的数据进行清理。
第一步就是比对。对于DNA测序来说,主讲人推荐使用BWA进行比对。对于RNA测序,推荐使用STAR进行比对:
在比对的过程中,有的reads能很好的识别来源,比如下面的region1和region2。也许会有一些local variants,但是仍然可以分辨出是出自于哪一个区域。但是有些时候,你的样品里会有重复的区域,这种duplicates会使得情况变得复杂。
在大多数的情况下,你的参考基因组是非常特异的,在比对的时候就很容易把reads比对到特异的基因组上。但是如果你的参考基因组本身就有很多重复区域,在比对的时候,软件就会很“困惑”,不知道要把reads比对到哪里去。比如下图中的2A和2B。这时我们就需要双端测序。
因为每一对read的长度是已知的,所以就可以确定到底哪一对read才是真正比对到基因组上的。
比对后生成的文件是SAM文件,以及它的二进制文件BAM文件。这些文件的格式有两部分,一个是头文件,里面是一些信息,比如版本、参考序列、和一些read group的信息;另一部分是record信息,对于每一个read来讲都是唯一的,其中包括read名字,flowcell信息,lane信息,read的方向,比对到正链还是负链,read的位置,比对质量。并且还有mate信息,它告诉你read比对到哪里,以及insert的大小。之后是read的序列,质量分数,最后是metadata,是一些read group的信息。
(这时有人提问,问的问题没有收音进来,只有回答:在一般情况下,在数据前期处理中,你不需要做任何trimming,你可以mark这些接头,之后这些接头会被忽略掉。传统的操作你还会做一些质量trimming,这也是完全没有必要做的,因为variants callers会把碱基质量进行计算,把其作为variant caller的模型。)
下面讲的是关于CIGAR字符。CIGAR字符很好的提供了每一条read的结构信息。下面是个例子:这个CIGAR字符所传递的信息就是1个soft clip,3个比对上的碱基,1个缺失,2个比对上,1个insertion,1个比对上。
对于RNA-seq,就是一个特殊的情况了。因为你的read可能被内含子被分开,如果你用DNA测序比对软件进行比对,你会得到大量的deletions的reads。所以需要另外的比对软件,比如说STAR。
下面讲第二步,标记duplicates。在RNA-seq数据里,你也需要标记这些duplicates。但是一般不选择去除这些duplicates,因为有时你会需要重新进行分析过程,如果直接去除了duplicates,你可能会丢失一些信息。还有在分析不同的批次的数据时,一定要使用相同的比对软件进行处理。
duplicates 的来源主要有两个:一个是文库制备过程中产生的,另一个是测序过程中产生的。
下面简单的说一下软件是如何“决定”是否把一条read标记为duplicate。一般是看一条Read的unclipped的5'端和read的方向(两个方向都会进行比对),并不是把read的全长进行比对来判断的。
下面是两张IGV的截图,分别显示和隐藏了duplicates的reads。软件只是把duplicates进行标记而已,并没有把它们删掉。这也方便让你在IGV里分别查看全部reads和部分reads。
第三步是重新校准碱基质量分数。为什么要做这一步呢?variant calling高度依赖你的数据里碱基质量。但是在测序过程中会存在一些测序仪系统偏差,你需要把这些从你的数据中清理掉。这个软件叫做BQSR,它可以识别测序仪的系统错误,BQSR会使用机器学习产生一个经验错误模型,然后把这个模型应用在样品里所有的reads上,从而对质量分数进行校准。
下面这张图展示的是,在进行重新校准质量分数后,你的比对文件里会出现新的质量分数。主讲人说他们通常把原始的质量分数去掉,来节省储存空间。