ITPub博客

首页 > 数据库 > > GATK流程--利用Pegasus : bwa算法简介

GATK流程--利用Pegasus : bwa算法简介

作者:飘絮絮絮丶 时间:2022-04-18 11:18:44 0 删除 编辑

一

Pegasus是一个Workflow Management System。

利用它,设计一个DAG图作为一连串任务的流程,方便的进行并行、串行,并可从错误中恢复。

二、更多信息

https://pegasus.isi.edu/about/

三、GATK

illumina NGS的原理 https://blog.csdn.net/u010608296/article/details/88831797

Gatk Best Practice:https://software.broadinstitute.org/gatk/best-practices/workflow?id=11165

 ① 将Reads比对到参考基因组上:BWA

BWA :Burrows-Wheeler Aligner

原理简述:

主要为了在内存消耗低的情况下实现序列比对(自然,动态规划是不行的)

Burrows-Wheeler数据压缩:

一种很自然的压缩方法是:ABAABADBADAD => 6A3B3D, 但是DNA序列不能这么压缩,因为DNA序列的信息蕴含在了碱基之间的顺序中。

现在将序列“panamabananas”做如图的Burrows-Wheeler Transform,获得的矩阵按照首字母排序,可以发现最后一列在原始字串中是位于第一列之前。可以看到,最后一列有多个a连续出现,这是因为,在原始序列中 =,‘an’这个2-mer多次出现(在DNA序列中,这种情况很多),因此,我们对最后一列进行压缩:SMNPB2N5A$A可以起到压缩作用(注意,ATCG组成的超长DNA中,这种压缩效果更显著)。由最后一列,经过首字母排序,可以推导第一列。因此,在比对时,我们只用保存参考基因组的Burrows-Wheeler Matrix 的最后一列的压缩形式。

 右图:由于英文文章中,‘and’的出现次数很多,所以他的BWT变换最后一列会出现连续的a(n、d亦然)。

 

通过第一列与最后一列,我们可以推导出整一个序列:最后一列在原始字串中是位于第一列之前,所以,从最后一列的$开始,看第一列的字符,然后再于最后一列找刚才的字符,再看第一列。。。

 „

在将read比对到参考基因组上时,我们从read的最后一位于BWMatrix最后一列开始寻找,结果如右图所示:在最开始的时候,我们为每一行算出一个SuffixArray(Text)作为定位参考。

 

现实中的比对是有错配的。在从第一列到最后一列的对应搜寻中,就可以相应放宽要求。过程如下图,在第一次搜寻,引入一个错配,最底下的情况,与第二次搜寻,又引入了一个错配,那么即停止并丢弃。

 

参看Phillip Compeau & Pavel Pevzner的Bioinformatics Algorithms: An Active Learning Approach

 

 

1
bwa mem ref.fa read1.fq read2.fq > aln-pe.sam

SN: Reference sequence name. The SN tags and all individual AN names in all @SQ lines must be distinct. The value of this field is used in the alignment records in RNAME and RNEXT fields. Regular expression: [:rname:∧*=][:rname:]*

LN: Reference sequence name. The SN tags and all individual AN names in all @SQ lines must be distinct. The value of this field is used in the alignment records in RNAME and RNEXT fields. Regular expression: [:rname:∧*=][:rname:]*

 

这步生成sam文件,然后要利用samtools将sam转为bam

SAM:Sequence Alignment Map
BAM:Binary Alignment Map
https://samtools.github.io/hts-specs/SAMv1.pdf

 

如果是将fastq文件分成多个来bwa的,现在要用picard将他们合并

1
java -jar picard.jar MergeSamFiles <argument...>

 

 

②Mark Duplicates:MarkDuplicatesSpark

  • Reads are sorted into coorrdinate-order
  • Mark the Duplicates, causing the marked pairs to be ignored by default during the variant discovery process.

   它的意义参见 https://www.jianshu.com/p/cad6a359a368:SNP的发现

 View Code

 

❗️在这一步之前要加上头信息:利用picard的AddOrReplaceReadGroups

  这个信息随便填也行,可以直接用示例

1
2
3
4
5
6
7
8
java -jar picard.jar AddOrReplaceReadGroups \
    I=input.bam \
    O=output.bam \
    RGID=4 \
    RGLB=lib1 \
    RGPL=ILLUMINA \
    RGPU=unit1 \
    RGSM=20

Mark Duplicates:

1
2
3
gatk MarkDuplicatesSpark \
            -I input.bam \
            -O marked_duplicates.bam

 

③ Base (Quality Score) Recalibration

  • 这一步的目的是利用机器学习方法纠正由测序仪给出的质量分数(Base Quality Score)
  • 这个误差来自:文库的制备,芯片生产,测序仪
  • 校准方法:
    • 收集covariate measurement from all base calls in datasets,建立统计模型。
      • The initial statistics collection can be parallelized by scattering across genomic coordinates, typically by chromosome or batches of chromosomes but this can be broken down further to boost throughput if needed.
      • Then the per-region statistics must be gathered into a single genome-wide model of covariation; this cannot be parallelized but it is computationally trivial, and therefore not a bottleneck.
    • 利用这个模型对碱基质量进行校准 (可以并行) 

 四、利用Pegasus来跑GATK

来自 “ ITPUB博客 ” ,链接:https://blog.itpub.net/69946831/viewspace-3125796/ ,如需转载,请注明出处,否则将追究法律责任。

请登录后发表评论 登录
全部评论

注册时间:2019-09-09

  • 博文量
    2646
  • 访问量
    42259