生信小白学习日记Day7——WGS分析流程（picard）

2019年6月2日，周日，天气晴，pass 上午。开始学习NGS分析，继BWA比对和SAM文件排序转BAM后的流程。

NGS分析

step5 Mark Duplications

参考这篇：GATK使用方法详解。https://www.plob.org/article/7009.html
在制备文库的过程中，由于PCR扩增过程中会存在一些偏差，也就是说有的序列会被过量扩增。这样，在比对的时候，这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates，这一步可以使用picard-tools来完成。去重复的过程是给这些序列设置一个flag以标志它们，方便GATK的识别。还可以设置 REMOVE_DUPLICATES=true 来丢弃duplicated序列。对于是否选择标记或者删除，对结果应该没有什么影响，GATK官方流程里面给出的例子是仅做标记不删除。这里定义的重复序列是这样的：如果两条reads具有相同的长度而且比对到了基因组的同一位置，那么就认为这样的reads是由PCR扩增而来，就会被GATK标记。
e.g.：

java -jar picard-tools-1.96/MarkDuplicates.jar

REMOVE_DUPLICATES= false

MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000

INPUT=hg19.reorder.sort.addhead_03.bam

OUTPUT=hg19.reorder.sort.addhead.dedup_04.bam METRICS_FILE=hg19.reorder.sort.addhead.dedup_04.metrics

实验室常用脚本中还会加入：VALIDATION_STRINGENCY=LENIENT
去查了一下，不知道理解的对不对：在BWA 比对生成SAM文件时，将没有map到基因组上的read归到了ref以外的区域，其MAPQ值不为0，而Picard认为这些read是不应该出现的，所以会报错（“MAPQ should be 0 for unmapped read” or “CIGAR should have zero elements for unmapped read”）。如果想忽略报错的话，就使用这行代码。
BWA can produce SAM records that are marked as unmapped but have non-zero MAPQ and/or non-"*" CIGAR. Typically this is because BWA found an alignment for the read that hangs off the end of the reference sequence. Picard considers such input to be invalid. In general, this error can be suppressed in Picard programs by passing VALIDATION_STRINGENCY=LENIENT or VALIDATION_STRINGENCY=SILENT. For ValidateSamFile, you can pass the arguments IGNORE=INVALID_MAPPING_QUALITY IGNORE=INVALID_CIGAR.（https://sourceforge.net/p/picard/wiki/Main_Page/）
注意： dedup这一步只要在library层面上进行就可以了，例如一个sample如果建了多个库的话，对每个库进行dedup即可，不需要把所有库合成一个sample再进行dedup操作。其实并不能准确的定义被mask的reads到底是不是duplicates，重复序列的程度与测序深度和文库类型都有关系。最主要目的就是尽量减小文库构建时引入文库的PCR bias。

Step 6: Index Bam File

对上一步得到的结果生成索引文件, 可以用samtools完成，生成的索引后缀是bai。

e.g.

samtools index hg19.reorder.sort.addhead.dedup_04.bam

Step 7: Local realignment around indels

这一步的目的就是将比对到indel附近的reads进行局部重新比对，将比对的错误率降到最低。一般来说，绝大部分需要进行重新比对的基因组区域，都是因为插入/缺失的存在，因为在indel附近的比对会出现大量的碱基错配，这些碱基的错配很容易被误认为SNP。还有，在比对过程中，比对算法对于每一条read的处理都是独立的，不可能同时把多条reads与参考基因组比对来排错。因此，即使有一些reads能够正确的比对到indel，但那些恰恰比对到indel开始或者结束位置的read也会有很高的比对错误率，这都是需要重新比对的。Local realignment就是将由indel导致错配的区域进行重新比对，将indel附近的比对错误率降到最低。

主要分为两步