加载中…
正文 字体大小:

【原创】GATK使用方法详解(包含bwa使用)第一部分

(2014-03-03 11:07:29)
标签:

gatk

bwa

snp

indel

分类: 生物信息

由于新浪博客规定,每篇文章不可超过2万字符,因此分4篇发布。

一、使用GATK前须知事项:

1)对GATK的测试主要使用的是人类全基因组和外显子组的测序数据,而且全部是基于illumina数据格式,目前还没有提供其他格式文件(如Ion Torrent)或者实验设计(RNA-Seq)的分析方法。

2GATK是一个应用于前沿科学研究的软件,不断在更新和修正,因此,在使用GATK进行变异检测时,最好是下载最新的版本,目前的版本是2.8.12014-02-25)。下载网站:http://www.broadinstitute.org/gatk/download

3)在GATK使用过程中(见下面图),有些步骤需要用到已知变异信息,对于这些已知变异,GATK只提供了人类的已知变异信息,可以在GATKFTP站点下载(GATK resource bundle)。如果要研究的不是人类基因组,需要自行构建已知变异,GATK提供了详细的构建方法。

4GATK在进行BQSRVQSR的过程中会使用到R软件绘制一些图,因此,在运行GATK之前最好先检查一下是否正确安装了R和所需要的包,所需要的包大概包括ggplot2gplotsbitopscaToolscolorspacegdatagsalibreshapeRColorBrewer等。如果画图时出现错误,会提示需要安装的包的名称。

 

二、GATK的使用流程

【原创】GATK使用方法详解(包含bwa使用)第一部分

GATK最佳使用方案:3大步骤。原始数据的处理变异检测初步分析。

 

第一大步:原始数据的处理

 

1. 对原始下机fastq文件进行过滤和比对(mapping

对于Illumina下机数据推荐使用bwa进行mapping

 

Bwa比对步骤大致如下:

1)对参考基因组构建索引:

     例子:bwa index -a bwtsw hg19.fa。最后生成文件:hg19.fa.ambhg19.fa.annhg19.fa.bwthg19.fa.pachg19.fa.sa

     构建索引时需要注意的问题:bwa构建索引有两种算法,两种算法都是基于BWT的,这两种算法通过参数-a is -a bwtsw进行选择。其中-a bwtsw对于短的参考序列是不工作的,必须要大于等于10Mb-a is是默认参数,这个参数不适用于大的参考序列,必须要小于等于2G

2)寻找输入reads文件的SA坐标。

     对于pair end数据,每个reads文件单独做运算,single end数据就不用说了,只有一个文件。

     例子:pair end

bwa  aln  hg19.fa  read1.fq.gz  -l 30  -k 2  -t 4  -I  > read1.fq.gz.sai

bwa  aln  hg19.fa  read2.fq.gz  -l 30  -k 2  -t 4  -I  > read2.fq.gz.sai

single end

bwa  aln  hg19.fa  read.fq.gz  -l 30  -k 2  -t 4  -I  > read.fq.gz.sai

主要参数说明:

-o int:允许出现的最大gap数。

-e int:每个gap允许的最大长度。

-d int:不允许在3’端出现大于多少bpdeletion

-i int:不允许在reads两端出现大于多少bpindel

-l intRead前多少个碱基作为seed,如果设置的seed大于read长度,将无法继续,最

好设置在25-35,与-k 2 配合使用。

-k int:在seed中的最大编辑距离,使用默认2,与-l配合使用。

-t int:要使用的线程数。

-R int:此参数只应用于pair end中,当没有出现大于此值的最佳比对结果时,将会降低标

准再次进行比对。增加这个值可以提高配对比对的准确率,但是同时会消耗更长的

时间,默认是32

-I int:表示输入的文件格式为Illumina 1.3+数据格式。

-B int:设置标记序列。从5’端开始多少个碱基作为标记序列,当-B为正值时,在比对之

前会将每个read的标记序列剪切,并将此标记序列表示在BC SAM 标签里,对于

pair end数据,两端的标记序列会被连接。

-b :指定输入格式为bam格式。bwa  aln  hg19.fa  read.bam  > read.fq.gz.sai

 

3)生成sam格式的比对文件。如果一条read比对到多个位置,会随机选择一种。

     例子:single endbwa  samse  hg19.fa  read.fq.gz.sai  read.fq.gz  > read.fq.gz.sam

           参数:-n int:如果reads比对次数超过多少次,就不在XA标签显示。

                 -r str:定义头文件‘@RG\tID:foo\tSM:bar’,如果在此步骤不进行头文件定义,在

后续GATK分析中还是需要重新增加头文件。

           pair endbwa sampe -a 500 read1.fq.gz.sai read2.fq.gz.sai read1.fq.gz read2.fq.gz > read.sam

           参数:-a int:最大插入片段大小。

                 -o intpair endreads中其中之一所允许配对的最大次数,超过该次数,将被视为

single end。降低这个参数,可以加快运算速度,对于少于30bpread,建

议降低-o值。

                 -r str:定义头文件。同single end

                 -n int:每对reads输出到结果中的最多比对数。

 

对于最后得到的sam文件,将比对上的结果提取出来(awk即可处理),即可直接用于GATK的分析。

注意:由于GATK在下游的snpcalling时,是按染色体进行callsnp的。因此,在准备原始sam文件时,可以先按染色体将文件分开,这样会提高运行速度。但是当数据量不足时,可能会影响后续的VQSR分析,这是需要注意的。

 

2. sam文件进行进行重新排序(reorder

BWA生成的sam文件时按字典式排序法进行的排序(lexicographically)进行排序的(chr10chr11…chr19chr1chr20…chr22chr2chr3…chrMchrXchrY),但是GATK在进行callsnp的时候是按照染色体组型(karyotypic)进行的(chrMchr1chr2…chr22chrXchrY),因此要对原始sam文件进行reorder。可以使用picard-tools中的ReorderSam完成。

e.g

java -jar picard-tools-1.96/ReorderSam.jar

I=hg19.sam

O=hg19.reorder_00.sam

REFERENCE=hg19.fa

 

注意:

1. 这一步的头文件可以人工加上,同时要确保头文件中有的序号在下面序列中也有对应的。虽然在GATK网站上的说明chrM可以在最前也可以在最后,但是当把chrM放在最后时可能会出错。

2. 在进行排序之前,要先构建参考序列的索引。

   e.g. samtools faidx hg19.fa。最后生成的索引文件:hg19.fa.fai

3. 如果在上一步想把大文件切分成小文件的时候,头文件可以自己手工加上,之后运行这一步就好了。

 

 

3. sam文件转换成bam文件(bam是二进制文件,运算速度快)

这一步可使用samtools view完成。

e.g. samtools view -bS hg19.reorder_00.sam -o hg19.sam_01.bam

 

4. bam文件进行sort排序处理

这一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。可以使用picard-toolsSortSam完成。

e.g.

java -jar picard-tools-1.96/SortSam.jar

INPUT=hg19.sam_01.bam

OUTPUT=hg19.sam.sort_02.bam

SORT_ORDER=coordinate

 

5. bam文件进行加头(head)处理

GATK2.0以上版本将不再支持无头文件的变异检测。加头这一步可以在BWA比对的时候进行,通过-r参数的选择可以完成。如果在BWA比对期间没有选择-r参数,可以增加这一步骤。可使用picard-toolsAddOrReplaceReadGroups完成。

e.g.

java -jar picard-tools-1.96/AddOrReplaceReadGroups.jar

I=hg19.sam.sort_02.bam

O=hg19.reorder.sort.addhead_03.bam

ID=hg19ID

LB=hg19ID

PL=illumine

PU=hg19PU

SM=hg19

ID str:输入readsID号;LBread集文库名;PL:测序平台(illunimasolid);PU:测序平台下级单位名称(run的名称);SM:样本名称。

注意:这一步尽量不要手动加头,本人尝试过多次手工加头,虽然看起来与软件加的头是一样的,但是程序却无法运行。

 

6. Merge

如果一个样本分为多个lane进行测序,那么在进行下一步之前可以将每个lanebam文件合并。

e.g.

java -jar  picard-tools-1.70/MergeSamFiles.jar

INPUT=lane1.bam

INPUT=lane2.bam

INPUT=lane3.bam

INPUT=lane4.bam

……

INPUT=lane8.bam

OUTPUT=sample.bam

 

7. Duplicates Marking

在制备文库的过程中,由于PCR扩增过程中会存在一些偏差,也就是说有的序列会被过量扩增。这样,在比对的时候,这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列,不能作为变异检测的证据,因此,要尽量去除这些由PCR扩增所形成的duplicates,这一步可以使用picard-tools来完成。去重复的过程是给这些序列设置一个flag以标志它们,方便GATK的识别。还可以设置 REMOVE_DUPLICATES=true 来丢弃duplicated序列。对于是否选择标记或者删除,对结果应该没有什么影响,GATK官方流程里面给出的例子是仅做标记不删除。这里定义的重复序列是这样的:如果两条reads具有相同的长度而且比对到了基因组的同一位置,那么就认为这样的reads是由PCR扩增而来,就会被GATK标记。

e.g.

java -jar picard-tools-1.96/MarkDuplicates.jar

REMOVE_DUPLICATES= false

MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000

INPUT=hg19.reorder.sort.addhead_03.bam

OUTPUT=hg19.reorder.sort.addhead.dedup_04.bam METRICS_FILE=hg19.reorder.sort.addhead.dedup_04.metrics

注意: dedup这一步只要在library层面上进行就可以了,例如一个sample如果建了多个库的话,对每个库进行dedup即可,不需要把所有库合成一个sample再进行dedup操作。其实并不能准确的定义被maskreads到底是不是duplicates,重复序列的程度与测序深度和文库类型都有关系。最主要目的就是尽量减小文库构建时引入文库的PCR bias

 

8. 要对上一步得到的结果生成索引文件

可以用samtools完成,生成的索引后缀是bai

e.g.

samtools index hg19.reorder.sort.addhead.dedup_04.bam

 

作者:葡萄糖

欢迎加入群:分子进化与基因组学

0

阅读 评论 收藏 转载 喜欢 打印举报
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 电话:4006900000 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有