美文网首页生物信息学习NGS生物信息学与算法
一步一步学会撰写SNP Calling 流程(二)

一步一步学会撰写SNP Calling 流程(二)

作者: 正踪大米饭儿 | 来源:发表于2018-12-22 22:07 被阅读5次

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

    相关文章

      网友评论

        本文标题:一步一步学会撰写SNP Calling 流程(二)

        本文链接:https://www.haomeiwen.com/subject/kviwkqtx.html