bowtie2+Samtools+GATK使用方法简介

1) Phase I: Raw data processing (对于mapping结果的处理)

 

(1)   通过mapping得到原始的mapping结果

GATK要求输入文件是BAM格式的，我们使用bowtie2进行mapping，bowtie2速度上比较有优势，而且输出的结果就是SAM格式的，用samtools转换成BAM格式即可，以下为bowtie2进行mapping的整个流程代码大致如下：

以S. cerevisiae为例，先下载fasta格式的参考基因组：

##Step1：用bowtie2-build生成索引文件，前缀名为s_cerevisiae

./bowtie2-build   s_cerevisiae.fasta    s_cerevisiae

##Step2： 这个和具体情况关系很大，参数自己弄懂后再决定。

基本命令：

bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} -S [<hit>]

必须参数：

-x <bt2-idx> 由bowtie2-build所生成的索引文件的前缀。首先 在当前目录搜寻，然后
在环境变量 BOWTIE2_INDEXES 中制定的文件夹中搜寻。
-1 <m1> 双末端测寻对应的文件1。可以为多个文件，并用逗号分开；多个文件必须和 -2 
<m2> 中制定的文件一一对应。比如:"-1 flyA_1.fq,flyB_1.fq -2 flyA_2.fq,flyB
_2.fq". 测序文件中的reads的长度可以不一样。
-2 <m2> 双末端测寻对应的文件2.
-U <r> 非双末端测寻对应的文件。可以为多个文件，并用逗号分开。测序文件中的reads的
长度可以不一样。
-S <hit> 所生成的SAM格式的文件前缀。默认是输入到标准输出。

以下是可选项：

输入

-q 输入的文件为FASTQ格式文件，此项为默认值。
-qseq 输入的文件为QSEQ格式文件。
-f 输入的文件为FASTA格式文件。选择此项时，表示--ignore-quals也被选择了。
-r 输入的文件中，每一行代表一条序列，没有序列名和测序质量等。选择此项时，表示--
ignore-quals也被选择了。
-c 后直接为比对的reads序列，而不是包含序列的文件名。序列间用逗号隔开。选择此项时，
表示—ignore-quals也被选择了。
-s/--skip <int> input的reads中，跳过前<int>个reads或者pairs。
-u/--qupto <int> 只比对前<int>个reads或者pairs（在跳过前<int>个reads或者
pairs后）。Default: no limit.
-5/--trim5 <int> 剪掉5'端<int>长度的碱基，再用于比对。(default: 0).
-3/--trim3 <int> 剪掉3'端<int>长度的碱基，再用于比对。(default: 0).
--phred33 输入的碱基质量等于ASCII码值加上33. 在最近的illumina pipiline中
得以运用。
--phred64 输入的碱基质量等于ASCII码值加上64.
--solexa-quals 将Solexa的碱基质量转换为Phred。在老的GA Pipeline版本中得以
运用。Default: off.
 --int-quals 输入文件中的碱基质量为用“ ”分隔的数值，而不是ASCII码。比如 40 40 
30 40...。Default: off.

                
## step3：把SAM格式转换成BAM格式
samtools view -bS sample01.sam -o sample01.bam  
## step4：对bam文件进行排序

如果没有readgroup需要先添加readgroup,命令是：

java  -jar AddOrReplaceReadGroups.jar I=sample.bam O=sample_addGroup.bam SORT_ORDER=coordinate CREATE_INDEX=true RGPL=illumina RGID=184 RGSM=sample184 RGLB=bar RGPU=pu184 VALIDATION_STRINGENCY=LENIENT

排序的命令为：

java -Xmx4g -jar   picard-tools-1.70/SortSam.jar \
            INPUT=sample01.bam \
            OUTPUT=sample01.sort.bam \
            SORT_ORDER=coordinate
这一步对bam文件进行排序的时候用的是picard里面的SortSam这个工具，而不是直接用是samtools来进行sort，因为samtools sort排完序以后不会更新BAM文件的头文件，SO这个tag显示的可能还是unsorted，下游程序处理的时候可能会报错；SORT_ORDER这边设置成coordinate表示按照染色体位置进行排序。

##step5： reference 是.fasta格式的需建立.dict文件，基本命令如下：

java -jar CreateSequenceDictionary.jar R= Homo_sapiens_assembly18.fasta O= Homo_sapiens_assembly18.dict

##step6：index模板基因组(reference 基因和比对基因组)

建立reference的index： samtools faidx Homo_sapiens_assembly19.fasta 

建立比对基因的基因模板：samtools index sorted.bam

(2)   处理原始的mapping结果

GATK提供了几种不同层次的处理方法，这边只讲一下最好也是最耗时的一种方式，其他的几种在此基础上调整一些步骤即可实现。

Best: multi-sample realignment with known sites and recalibration

同时考虑多个sample来进行本地重排以及校正，下面简单的介绍一下涉及到的几个步骤。

注: 下面的有些步骤虽然提到要合并文件，但有的工具本身可以指定多个输入文件的话实际上是不需要专门先去合并的，得到的结果就直接是一个合并好的文件，具体情况请根据实际碰到的为准。

 

for each sample

    lanes.bam <- merged lane.bam's for sample 

    第一步是把同一个sample测了多个lane的结果进行合并，合并可以用picard里面的MergeSamFiles工具来实现。

     例如我们要把lane1.bam, lane2.bam, lane3.bam, lane4.bam合并成sample.bam，可以用

## 合并每个sample的bam文件
java -Xmx4g -jar   picard-tools-1.70/MergeSamFiles.jar \
        INPUT=lane1.bam \
        INPUT=lane2.bam \
        INPUT=lane3.bam \
        INPUT=lane4.bam \
        OUTPUT=sample.bam

因为我们这边每个sample只有一个lane，所以这一步跳过。

 

    dedup.bam <- MarkDuplicates(lanes.bam)

    这一步是对每个sample的一些重复序列进行一些处理，这些重复的序列可能是因为PCR扩增的时候引入的一些引物序列，容易干扰下游结果，这个操作可以通过picard里面的MarkDuplicates工具实现，MarkDuplicates处理完成不会删除这些reads而是加上一个标签，带有这些标签的reads在后续处理的时候会被跳过。(注：dedup这一步只要在library层面上进行就可以了，例如一个sample如果建了多个库的话，对每个库进行dedup即可，不需要把所有库合成一个sample再进行dedup操作)

     以sample01为例，命令行代码如下：

## 用picard里面的MarkDuplicates工具去除PCR duplicates，这一步会生成
## 两个文件，一个dedup好的bam文件以及一个纯文本文件，这边起名叫
## sample01.dedup.metrics，里面包含duplicates的一些统计信息
java –Xmx4g -jar picard-tools-1.70/MarkDuplicates.jar \
        MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 \
        INPUT= sample01.sort.bam \
        OUTPUT= sample01.dedup.bam \
        METRICS_FILE= sample01.dedup.metrics

    realigned.bam <- realign(dedup.bam) [with known sites included if available]

这一步是sample层次的一个本地重排，因为在indel周围出现mapping错误的概率很高，而通过在indel周围进行一些重排可以大大减少由于indel而导致的很多SNP的假阳性，如果已经有一些可靠位点的信息的话 (例如dbsnp的数据)，这一步可以只在一些已知的可靠位点周围进行，这样得到的结果比较可靠也比较节省时间；而缺少这些数据的情况下则需要对所有位点都进行这样一个操作。

(题外话：之前版本v3的时候realign是在cohort层次进行的，不过相当耗时耗力，所以到这个版本就废弃了。。。(当然也有一些特例可能需要综合考虑才能得到更好的结果)

realign这一步用的是GATK里面的RealignerTargetCreator和IndelRealigner工具，如果已经有可靠的INDEL位点(如dbsnp数据等等，需要是VCF格式的)信息，可以通过-known参数来进行设置，让realign操作主要集中在这些位点附近，实际情况下可能很多物种并没有这样的参考数据，所以需要对所有INDEL位点进行这样的一个操作：

## 对INDEL周围进行realignment，这个操作有两个步骤，第一步生成需要进行
## realignment的位置的信息，第二步才是对这些位置进行realignment，最后生
## 成一个realin好的BAM文件sample01.realn.bam
java -Xmx4g -jar GenomeAnalysisTK.jar \
            -R ref.fasta \
            -T RealignerTargetCreator \
            -o sample01.realn.intervals \
            -I sample01.dedup.bam
java -Xmx4g -jar GenomeAnalysisTK.jar \
        -R ref.fasta \
        -T IndelRealigner \
        -targetIntervals sample01.realn.intervals \
        -o sample01.realn.bam \
        -I sample01.dedup.bam

     recal.bam <- recal(realigned.bam)

     sample.bam <- realigned.bam

最后一步是对realign得到的结果进行校正，因为mapping得到的结果的分值往往不能反应真实情况下的分值，而只有尽量消除这种偏态才能保证得到的结果的可靠性与正确性，所以如果情况允许的条件下都需要进行这样的一个校正步骤。这样最后我们就得到一个对于每个sample而言最好的处理结果。

校正这一步对低覆盖度的数据比较重要，可以消除很多假阳性；对于高覆盖度(10x~)的数据可以不做这一步，这一步最重要的是有已知的可靠位点作为参考， human因为研究的比较多，所以有很多这样的数据(包括dbsnp的数据等等)，而其他的物种则可能缺少这些数据，这个时候我们可以考虑先做一次初步的variants calling，然后筛选出我们认为最可靠的那些位点作为参考位点，用这部分的数据再来进行校正，下面提供一种用来生成参考位点的数据的方法，仅供参考：
## step1: variants calling by GATK
java -jar GenomeAnalysisTK.jar \
    -R ref.fasta \
    -T UnifiedGenotyper \
    -I sample01.realn.bam \
    -o sample01.gatk.raw.vcf \
    -stand_call_conf 30.0 \
    -stand_emit_conf 0 \
    -glm BOTH \
    -rf BadCigar
这边有一个-rf参数，是用来过滤掉不符合要求的reads，这边是把包含错误的Cigar字符串的reads给排除掉，关于Cigar字符串可以参考关于sam文件的说明(The SAM Format Speci?cation)，sam文件的第六行就是这边的Cigar字符串，-rf的其他参数可以参考GATK网站Read filters下面的条目http://www.broadinstitute.org/gatk/gatkdocs/
## step2: variants calling by samtools
samtools mpileup -DSugf ref.fasta sample01.realn.bam | \
    bcftools view -Ncvg - \
    > sample01.samtools.raw.vcf
## step3: 选取GATK和samtools一致的结果
java -Xmx4g -jar GenomeAnalysisTK.jar \
    -R ref.fasta \
    -T SelectVariants \
    --variant sample01.gatk.raw.vcf \
    --concordance sample01.samtools.raw.vcf \
    -o sample01.concordance.raw.vcf
## step4：筛选上面得到的结果，这边filter用到的几个标准可以根据实际情况
## 灵活更改，对QUAL值的筛选用的是$MEANQUAL，表示所有QUAL值的平均
## 值，linux底下这个值可以通过第一条命令行得到
## 计算平均值
MEANQUAL=`awk '{ if ($1 !~ /#/) { total += $6; count++ } }
                END { print total/count }' sample01.concordance.raw.vcf`
## 筛选
java -Xmx4g -jar GenomeAnalysisTK.jar \
    -R ref.fasta \
    -T VariantFiltration \
    --filterExpression " QD < 20.0 || ReadPosRankSum < -8.0 || FS > 10.0 || QUAL < $MEANQUAL" \
    --filterName LowQualFilter \
    --missingValuesInExpressio nsShouldEvaluateAsFailin g \
    --variant sample01.concordance.raw.vcf \
    --logging_level ERROR \
    -o sample01.concordance.flt.vcf
## step5：提取通过筛选标准的位点到结果文件中
grep -v "Filter" sample01.concordance.flt.vcf > sample01.confidence.raw.vcf

最后这个sample01.confidence.raw.vcf的结果就可以当作可靠的参考位点用于后面的分析。

校正这个步骤如果使用的是GATK 1.x版本的话需要使用CountCovariates和TableRecalibration等工具，2.x版本取消了这两个工具而改成BaseRecalibrator，下面都给出命令行代码：

## GATK 1.x版本
## step1：生成校正要用的文件sample01.recal_da ta.pre.csv
java -Xmx4g -jar GenomeAnalysisTK.jar \
    -R ref.fasta \
    -T CountCovariates \
    -cov ReadGroupCovariate \
    -cov QualityScoreCovariate \
    -cov CycleCovariate \
    -cov DinucCovariate \
    -knownSites sample01.confidence.raw.vcf \
    -I sample01.realn.bam \
    -recalFile sample01.recal_da ta.pre.csv
## step2：这一步和最后两步也可以跳过，主要作用是生成一系列的统计图表，
## 在做完校正之后可以再次生成上一步的文件，这样可以比较清楚的看到校正
## 前后的差别
java -Xmx4g -jar AnalyzeCovariates.jar \
    -recalFile sample01.recal_da ta.pre.csv \
    -outputDir sample01.graphs.pre \
    -ignoreQ 5
## step3：进行校正，生成校正后的文件sample01.recal.bam
java -Xmx4g -jar GenomeAnalysisTK.jar \
    -R ref.fasta \
    -T TableRecalibration \
    -I sample01.realn.bam \
    -recalFile sample01.recal_da ta.pre.csv \
    -o sample01.recal.bam
## step4：再次生成校正数据sample01.recal_da ta.post.csv
java -Xmx4g -jar GenomeAnalysisTK.jar \
    -R ref.fasta \
    -T CountCovariates \
    -cov ReadGroupCovariate \
    -cov QualityScoreCovariate \
    -cov CycleCovariate \
    -cov DinucCovariate \
    -knownSites sample01.confidence.raw.vcf \
    -I sample01.recal.bam \
    -recalFile sample01.recal_da ta.post.csv
## step5：生成统计图表
java -Xmx4g -jar AnalyzeCovariates.jar \
    -recalFile sample01.recal_da ta.post.csv \
    -outputDir sample01.graphs.post \
    -ignoreQ 5

 

2.x版本改动最大的就是淘汰了原来用于生成校正文件以及进行校正的两个工具：CountCovariates和TableRecalibration，取而代之的是新版本里面的BaseRecalibrator与PrintReads这两个工具，主要影响的就是原来的step1以及step3，这边就只给出这两步的代码：

## step1：生成校正要用的文件sample01.recal_da ta.grp
java -Xmx4g -jar GenomeAnalysisTK.jar \
    -R ref.fasta \
  -T BaseRecalibrator \
  -I sample01.realn.bam \
  -knownSites sample01.confidence.raw.vcf \
  -o sample01.recal_da ta.grp
这一步会生成好几个文件，包括关于分值的一些统计报告。另外这边也有-cov这个参数，可以像上面一样设置，不管设不设ReadGroupCovariate和QualityScoreCovariate两个是一定会有的。
## step3：进行校正，生成校正后的文件sample01.recal.bam
java -jar GenomeAnalysisTK.jar \
  -R ref.fasta \
  -T PrintReads \
  -I sample01.realn.bam \
  -BQSR sample01.recal_da ta.grp \
  -o sample01.recal.bam

下面是GATK网站上给出的几张校正前后的比较图：