Broad Institute视频笔记 Introduction
本篇笔记是这个系列视频的第4个视频,我没有把所有的视频PPT放上面,因为有些PPT内容很空洞,所以只记录了其中一些。有兴趣的同学可以观看英文原版视频:here
上一个视频里讲了数据的预处理过程,完成了预处理后,我们拿到手的bam文件里有了重新校准的碱基质量、标记了duplicates,这时就可以准备call variants了。这个视频里我们主要是介绍call 几种不同的variants。再介绍一下variants caller输出文件。
大家都知道人类基因组大概有30亿个位点。其中有99.5%的DNA是一样的。只有其中很小一部分的variants在个体之间是不同的。
那么我们怎么才能把某一个特定的样品里的variants表示出来呢?在我们的分析中,参考基因组是很重要的,要注意你使用的基因组文件的版本,因为不同版本之间会有一些variants的不同。最好是使用最新版本的参考基因组来处理你的数据。
有几种不同的variants。这里我们讨论的大部分的variants是short variants,比如说SNPs和Indels。Indels的大小大约小于50bp,大于50bp的indels在haplotype caller是不会对其回应的。另一种variant是拷贝数的variants,还有一种就是移位。对于上一个视频里讲了数据的预处理过程,完成了预处理后,我们拿到手的bam文件里有了重新校准的碱基质量、标记了duplicates,这时就可以准备call variants了。这个视频里我们主要是介绍call 几种不同的variants。再介绍一下variants caller输出文件。
大家都知道人类基因组大概有30亿个位点。其中有99.5%的DNA是一样的。只有其中很小一部分的variants在个体之间是不同的。
那么我们怎么才能把某一个特定的样品里的variants表示出来呢?在我们的分析中,参考基因组是很重要的,要注意你使用的基因组文件的版本,因为不同版本之间会有一些variants的不同。最好是使用最新版本的参考基因组来处理你的数据。
有几种不同的variants。这里我们讨论的大部分的variants是short variants,比如说SNPs和Indels。Indels的大小大约小于50bp,大于50bp的indels在haplotype caller是不会对其回应的。另一种variant是拷贝数的variants,还有一种就是移位。另外主讲人提到有一个tool叫PathSeq可以用来检测你的DNA样品中是否有pathogen的污染。
上面这张图显示的是在IGV里,导入bam文件后,variants的展示形式。灰色条表示与参考基因组比对上的碱基。
上面这张PPT展示的是VCF文件的格式。(这里主讲人讲的太差劲了,哼哼唧唧半天。。。我干脆就直接打开了上面PPT里放的链接,自己看原版的vcf格式说明)vcf文件就是记录variants的一种文件。vcf文件分为hearder部分和records部分。header部分都是以“##”开头,之后用“键=值”这样的形式来表示的。header的第一行,就是VCF format版本号。INFO部分的说明对应的是Records部分里的第8列,在第8列里有多少种情况,这里就有多少行INFO;同样的FORMAT信息有3行,对应的是records部分第9列,在第9列里有多少种情况就有多少行FORMAT。比如你看到上面的第9列里有GT,GQ和DP,如果你想知道它们分别代表什么意思,就去上面的header部分的FORMAT3行里的“description”去找。
而record部分里,每一列分别代表:
(1)#CHROM:这个variant所在的染色体编号。
(2)POS:variant的位点在该染色体的位置。这里要注意的是,可以有多个variants有相同的POS。
(3)ID:如果这个variant是dbSNP数据库里的,就用rs加数字来表示。如果不存在dbSNP数据库里,用“.”表示,说明是一个新的variant。
(4)REF:参考基因组。每个碱基必须是A、C、G、T、N(大小写不敏感)中的一个。这里允许出现多个碱基。之前POS字段中的值指的是本列字符串中第一个碱基的位置。对于简单插入和缺失,字符串必须包含插入/缺失之前的碱基,除非该variant发生在contig上的第1个位置,在此情况下,必须包括该variant后的碱基。
(5)ALT:alternate base。这些等位基因不需要在所有样本中被call。由A、C、G、T、N、*(大小写不敏感)组成的基本字符串。
(6)QUAL:质量。Phred格式的质量分数。该值越高,代表这个位置的存在variant的可能性越大。如果前面的ALT列里是“.”,则这个值是−10log10(没有variant的可能性);如果前面的ALT列里不是“.”,这个值也是−10log10(有variant的可能性)。
(7)FILTER:如果这列显示的是PASS,说明这个位点call出来一个variant。
(8)INFO:一些附加信息。DP指的是测序深度,也就是说有多少条reads支持你的这个variant,比如上面这个图DP=14,说明在这个variant位点上有14条reads支持它的变异。(具体每一项都是什么,请参考官方网站:http://samtools.github.io/hts-specs/VCFv4.2.pdf),如果不想看英文原版的,这里有两篇文章可以参考:看懂变异记录结果文件(VCF)和vcf格式文件基础知识与编辑操作详解。
而对于somatic variant的分析和体细胞的有很大不同。比如肿瘤样品的call variants,因为我们在取肿瘤样品的时候,经常会取到周围的正常组织,所以这时我们数据的信噪比就会小得多,另外由于肿瘤的异质性,样品之间的variants也会有一些不一样。
为了降低背景噪音,我们可以通过增加样品数量,也可以通过与正常组织进行对照(大量的正常组织对照),这样就可以知道哪些variants是来自germline的。比如下面这张图:
在图中可以看到有些variants在germline里也存在,所以你可以推测哪些variants可能是肿瘤样品中特异存在的。
然后我们要说的就是coverage changes,覆盖度的变化。coverage profile的变化决定了是否发生了variant。所以你看到在某一个特定的区域由更多的reads,说明这是一个copy number variant的信号。
上面这张图就是对于体细胞的copy number variant分析的流程。在这个流程中比较重要的就是从去噪音到标准化coverage。对于外显子的数据就尤为重要。然后比较重要的步骤就是segmentation,这一步是决定你的events的边界在哪里。
视频里提到的所有的piplines(流程)代码都可以在官方github上找到:https://github.com/gatk-workflows