2019年5月26日下午,无意中看到hanli0902的关于NGS分析的博文https://blog.csdn.net/hanli1992/article/details/82790386有很多需要学习的地方,在这里贴一些并就不理解之处做些笔记,仅供参考。
NGS基础——NGS分析
NGS 分析步骤
1. 质量分析
fastqc、multiqc、SolexaQA
测序数据的质量好坏会影响我们的下游分析。但不同的测序平台其测序错误率的图谱都是有差别的。因此,非常建议在我们分析测序数据之前先搞清楚如下两个地方:
原始数据是通过哪种测序平台产生的,它们的错误率分布是怎么样的,是否有一定的偏向性和局限性,是否会显著受GC含量的影响等;评估它们有可能影响哪些方面的分析;
其中fastqc的使用和结果分析在Day1的博文中有简要介绍,另外两个软件的使用方法因为不是今天的重点,若有时间再去查阅相关文章。(在此标记出来,可能在下一篇博文中介绍)
2. 数据过滤
去除接头序列和低质量序列,常用工具:SOAPnuke、cutadapt、untrimmed、fastp、sickle、seqtk、Trimmomatic、Trim Galore(合并了FastQC和cutadapt)同样标记,具体细节查阅后再介绍
3. 比对
使用BWA(exact match)软件把这些短序列和参考基因组进行对比, 确定短序列在基因组上的位置,生成*.sam文件。标记
4. 排序
用samtools等进行比对、排序、建索引,得到reads在参考序列中的位置及质量值,并进行数据格式的转换,得到*.bam文件。标记
5. Bam文件再处理
remove duplicates → indel realign → BQSR
remove duplicates
使用Picard软件(标记)把测序产生的冗余信息和噪声去掉,并对数据质量进行评价。
在制备文库的过程中,由于PCR扩增过程中会存在一些偏差,也就是说有的序列会被过量扩增。这样,在比对的时候,这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列,不能作为变异检测的证据,因此,要尽量去除这些由PCR扩增所形成的duplicates,这一步可以使用picard-tools来完成。去重复的过程是给这些序列设置一个flag以标志它们,方便GATK的识别。还可以设置 REMOVE_DUPLICATES=true来丢弃duplicated序列。对于是否选择标记或者删除,对结果应该没有什么影响,GATK官方流程里面给出的例子是仅做标记不删除。这里定义的重复序列是这样的:如果两条reads具有相同的长度而且比对到了基因组的同一位置,那么就认为这样的reads是由PCR扩增而来,就会被GATK标记。
indel realign
BWA采取exact match策略,在indel附近比对效果不好,需进行局部重新比对。一般来说,绝大部分需要进行重新比对的基因组区域,都是因为插入/缺失的存在,因为在indel附近的比对会出现大量的碱基错配,这些碱基的错配很容易被误认为SNP。还有,在比对过程中,比对算法对于每一条read的处理都是独立的,不可能同时把多条reads与参考基因组比对来排错。因此,即使有一些reads能够正确的比对到indel,但那些恰恰比对到indel开始或者结束位置的read也会有很高的比对错误率,这都是需要重新比对的。Local realignment就是将由indel导致错配的区域进行重新比对,将indel附近的比对错误率降到最低。可以使用GATK工具进行indel重新比对。
BQSR(Base Quality Score Recalibration)
BQSR是对bam文件里reads的碱基质量值进行重新校正,使最后输出的bam文件中reads中碱基的质量值能够更加接近真实的与参考基因组之间错配的概率。例如,在reads碱基质量值被校正之前,我们要保留质量值在Q25以上的碱基,但是实际上质量值在Q25的这些碱基的错误率在1%,也就是说质量值只有Q20,这样就会对后续的变异检测的可信度造成影响。还有,在边合成边测序的测序过程中,在reads末端碱基的错误率往往要比起始部位更高。另外,AC的质量值往往要低于TG。BQSR的就是要对这些质量值进行校正。使用GATK对重新比对的BAM文件做质量校准。
6. 获取突变
用GATK等得到**.vcf**文件,从中获取如外显子区域的数据,错义突变,热点突变,靶向药物相关的突变等感兴趣的突变信息。
7. 注释
使用Annovar对这些变异位点进行功能注释(标记), 得到一个易于理解的变异位点列表。
下面就一些格式文档进行说明,同样来自hanli0902的博文:https://blog.csdn.net/hanli1992/article/details/82790386
SAM
SAM/BAM格式专用于存储基于参考序列的比对序列,SAM(Sequence Alignment Map)是“序列比对映射”的首字母缩写,是带有比对信息的序列文件(即告诉你这个reads在染色体上的位置等),用于储存序列数据。
经过专门的比对软件,如BWA、BOWTIE2(标记)等,得到的SAM文件如下所示,需要研究的就是如下这几行。
Coor:坐标的简写,方便查看比对
ref:参考序列
r001/2表示paired end数据,
r003是嵌合read,
r004则是原序列打断后比对结果。
这里没太看懂,有时间再查一下
BAM
BAM是(SAM的)二进制格式,因为它是压缩的、所以数据量更小;因为它是有索引的,所以可以更快地访问它。
VCF(Variant Call Format)
VCF文件分为两部分内容:以“#”开头的注释部分;没有“#”开头的主体部分。
CHROM和POS
:代表参考序列名和variant的位置;如果是INDEL的话,位置是INDEL的第一个碱基位置。
ID
:variant的ID。比如在dbSNP中有该SNP的id,则会在此行给出;若没有,则用’.'表示其为一个novel variant。
REF
和 ALT
:参考序列的碱基 和 Variant的碱基。
QUAL
:Phred格式(Phred_scaled)的质量值,表 示在该位点存在variant的可能性;该值越高,则variant的可能性越大;计算方法:Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1,该位点为variant的概率为90%。
FILTER
:使用上一个QUAL值来进行过滤的话,是不够的。GATK能使用其它的方法来进行过滤,过滤结果中通过则该值为”PASS”;若variant不可靠,则该项不为”PASS”或”.”。
INFO
: 这一行是variant的详细信息。
FORMAT
和 NA12878
:这两行合起来提供了’NA12878′这个sample的基因型的信息。’NA12878′代表这该名称的样品,是由BAM文件中的@RG下的 SM 标签决定的。
GT
:样品的基因型(genotype)。两个数字中间用 /
分 开,这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele; 1 表示样品中variant的allele; 2表示有第二个variant的allele。因此: 0/0 表示sample中该位点为纯合的,和ref一致; 0/1 表示sample中该位点为杂合的,有ref和variant两个基因型; 1/1 表示sample中该位点为纯合的,和variant一致。
AD
和 DP
:AD(Allele Depth)为sample中每一种allele的reads覆盖度,在diploid中则是用逗号分割的两个值,前者对应ref基因型,后者对应variant基因型; DP(Depth)为sample中该位点的覆盖度。
GQ
:基因型的质量值(Genotype Quality)。Phred格式(Phred_scaled)的质量值,表示在该位点该基因型存在的可能性;该值越高,则Genotype的可能性越 大;计算方法:Phred值 = -10 * log (1-p) p为基因型存在的概率。
PL
:指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1),这三种基因型的概率总和为1。和之前不一致,该值越大,表明为该种基因型的可能 性越小。 Phred值 = -10 * log § p为基因型存在的概率。
身体不适,有点急躁,明天再来