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 标记重复比对序列并去除(见第三部分)。
网友评论