首页 > 数据库 > > GATK流程--利用Pegasus : bwa算法简介
一
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
原理简述:
主要为了在内存消耗低的情况下实现序列比对(自然,动态规划是不行的)
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
它的意义参见 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
四、利用Pegasus来跑GATK

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