一步一步学会撰写SNP Calling 流程（二）
3.1 DNASeq 数据分析SNP/Indel 策略
策略一： BWA + Samtools /Picard + GATK
官宣：GATK 分析流程
1) BWA 比对及 Samtools 转化为 bam 格式
BWA 是李恒大神写的一款比对软件。
BWA MEM比对模块是有一定适用范围的：它是专门为长read比对设计的，目的是为了解决，第三代测序技术这种能够产生长达几十kb甚至几Mbp的read情况。一般只有当read长度≥70bp的时候，才推荐使用，如果比这个要小，建议使用BWA ALN模块。
## ======== Step 2 bwa alignment ========
rule bwaAlign:
input:
R1 = WORKDIR + "Step1.fastqFilter/{sample}/{sample}.R1.fq.gz",
R2 = WORKDIR + "Step1.fastqFilter/{sample}/{sample}.R2.fq.gz"
output:
bam = WORKDIR + "Step2.bwaAlign/{sample}.bam"
log:
WORKDIR + "logs/Step2.bwaAlign/{sample}.align.logs"
benchmark:
WORKDIR + "benchmark/Step2.bwaAlign/{sample}.benchmark"
threads:
8
params:
"-k 19 -M -Y"
shell:
"bwa mem -t {threads} {params} {BWA_INDEX} {input.R1} {input.R2} |samtools view -@ {threads} -Sb - >{output.bam} 2> {log}"
2) Picard 数据处理
由于在 BWA 比对过程中我们并没有设置 bam 的头信息。因此我们使用 Picard 添加 bam 头文件并排序。添加的头信息是以 @RG 开头的一行进行存储的，其主要功能是将比对的 reads 进行分组，不同的组之间的测序过程被认为是相互独立的，分组信息为后续 Mark Duplicate 提供一定的依据。主要有以下几个部分：
1. ID，这是Read Group的分组ID，一般设置为测序的lane ID (不同lane之间的测序过程认为是独立的)，下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；
2. PL，指的是所用的测序平台，这个信息不要随便写！特别是当我们需要使用GATK进行后续分析的时候，更是如此！这是一个很多新手都容易忽视的一个地方。
3. SM，样本ID，同样非常重要，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的 lane 分布测出来，这个时候SM名字就是可以用于区分这些样本。
4. LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的 fq 文件名中找到，如果上面的 lane ID 足够用于区分的话，也可以不用设置LB；
除了以上这四个之外，还可以自定义添加其他的信息，不过如无特殊的需要，对于序列比对而言，这4个就足够了。这些信息设置好之后，在RG字符串中要用制表符（\t）将它们分开。
在GATK中，PL只允许被设置为：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS 或 UNKNOWN 这几个信息。基本上就是目前市场上存在着的测序平台，当然，如果实在不知道，那么必须设置为 UNKNOWN，名字方面不区分大小写。如果你在分析的时候这里没设置正确，那么在后续使用GATK过程中可能会碰到类似如下的错误：
ERROR MESSAGE: The platform (xx) associated with read group GATKSAMReadGroupRecord @RG:xx is not a recognized platform.
如果你的数据是CG测序的那么记得不要写成CG！而要写COMPLETE。
例如我们自己的一个 Illumina Nova 测序数据（截取）所示：
Illumina 测序数据示意
由 @ST 行可知：
| 描述 | 标识 |
|---|---|
| the unique instrument name | ST-E00578 |
| the run id | 294 |
| the flowcell id | HLKK3CCXY |
| flowcell lane | 6 |
| Tile number within the flowcell lane | 1101 |
| 'x'-coordinate of the cluster within the tile | 3853 |
| 'y'-coordinate of the cluster within the tile | 2047 |
| the member of a pair, 1 or 2 (paired-end or mate-pair reads only) | 1 |
| Y if the read is filtered, N otherwise | N |
| 0 when none of the control bits are on, otherwise it is an even number | 0 |
| index sequence | CAGCGTTA |
Snakemake 分析代码如下：
## ======== Step 3.1 AddOrReplaceReadGroups ========
rule AddOrReplaceReadGroups:
input:
bam = WORKDIR + "Step2.bwaAlign/{sample}.bam"
output:
temp(WORKDIR + "Step3.Picard/{sample}.add_rg_sort.bam")
log:
WORKDIR + "logs/Step3.Picard/{sample}.add_rg.log"
params:
rg = "SO=coordinate RGID={sample} RGLB=lib RGPL=illumina RGPU=hiseq RGSM={sample}"
shell:
"picard AddOrReplaceReadGroups I={input} O={output} {params.rg}"
ID：输入 Reads ID号（字符串类型）
LB：Reads 文库名
PL：测序平台 ( Illunima 或 Solid )
PU：测序平台下级单位名称（run的名称）
SM：样本名称
标记好 Header 信息后，使用 GATK MarkDuplicates 标记重复比对序列并去除（见第三部分）。
