GATK Best Practices for RNA-Seq

前言

GATK是由Broad institute开发的用于变异检测(例如: SNP, Indel等)高通量测序数据. 目前GATK已经被广泛应用于人类以及哺乳动物从WGS (Whole Genome Sequencing) 数据中检测遗传变异。本文讲展示如何使用GATK官方提供的Best Practice对群体RNA-Seq数据进行变异检测。

所需要的工具

• STAR
• sambamba
• GATK (v3.7 or v3.8)
• picard

详细步骤及解释

1、为参考基因组创建index

• 根据GATK官方提供的Best Practice,在众多软件中STAR拥有最高的Sensitivity和Specificity.

  STAR --runThreadN 20 --runMode genomeGenerate --genomeDir ~/genomes/Maize/maize_genome_4_35/STAR_idx/ --genomeFastaFiles ~/genomes/Maize/maize_genome_4_35/Zea_mays_AGPv4_dna_sm_chromosome.fasta --sjdbGTFfile ~/genomes/Maize/maize_genome_4_35/maize_v4.35.gtf --sjdbOverhang 100
  

• 上述所用参数解释

  --runThreadN: 指定线程数目
  --runMode：指定STAR子程序genomeGenerate
  --genomeDir: 指定参考基因组index路径
  --genomeFastaFiles: 参考基因组的fasta文件
  --sjdbGTFfile: 参考基因组GTF注释文件
  --sjdbOverhang: 该参数用来指定用于拼接junctions数据库时junctions周围基因组序列的长度。理想情况下, 这个值应该是reads长度减1。例如，对于Illumina双端reads（100bp）,这个值应该设置为99。但是大多数情况下，设置成100也能达到理想的效果。
  

2、STAR 2-pass alignment

• 2-pass alignment是指先使用指定参数比对到参考基因组，然后使用上一步生成的junction信息(SI.out.tab)重新创建基因组index，然后使进行第二轮比对

  #指定参考基因组index
  star_ref="~/genomes/Maize/maize_genome_4_35/STAR_idx"
  ref="~/Zea_mays_AGPv4_dna_sm_chromosome.fasta"
  #SRR accession
  bn="/your/path/SRR1234567"
  opdir=$bn"_Results"
  mkdir $opdir
  
  #Path to gatk and picard tools
  gatk="~/software/GATK_3.7/GenomeAnalysisTK.jar"
  picard="~/software/picard/build/libs/picard.jar"
  
  # time用于统计该行命令所需要的时间
  time STAR --outFileNamePrefix $opdir/$bn --outSAMtype BAM SortedByCoordinate --outFilterMultimapNmax 1 --outSAMstrandField intronMotif --genomeDir $star_ref --runThreadN 20  --readFilesIn $bn"_1_fastp_filter.fastq" $bn"_2_fastp_filter.fastq" --twopassMode Basic
  
  #sambamba index (~1min)
  echo -e "["$(date)"]\tIndexing.."
  time sambamba index -p -t 20 $opdir/$bn"Aligned.sortedByCoord.out.bam"
  

• 上述所用参数解释

  --outFileNamePrefix: 输出文件名前缀
  --outSAMtype: 指定输出SAM类型，该处使用排序过后的BAM (BAM SortedByCoordinate)
  --outFilterMultimapNmax: 指定输出唯一匹配的reads, 1代表uniquely mapped reads
  --outSAMstrandField:  对于非链特异性的RNA-Seq，Cufflinks/Cuffdiff需要指定此参数
  --genomeDir: 参考基因组index文件
  --runThreadN: 指定线程数目
  --readFilesIn: FASTQ文件
  --twopassMode: 指定2-pass比对模式,即先使用上述参数比对到参考基因组，然后使用上一步生成的junction信息(SI.out.tab)重新创建基因组index，然后使用上述参数进行第二轮比对
  

3、Add read groups

• 替换BAM文件中的read groups. GATK对于read groups的解释入下：There is no formal definition of what is a read group, but in practice, this term refers to a set of reads that were generated from a single run of a sequencing instrument.

  time java -XX:ParallelGCThreads=20 -jar $picard AddOrReplaceReadGroups I=$opdir/$bn"Aligned.sortedByCoord.out.bam" O=$opdir/$bn"Aligned.rg_added.sortedByCoord.out.bam" SO=coordinate RGID=id RGLB=library RGPL=platform RGPU=machine RGSM=sampleName 2>$opdir/$bn.AddReadGroups.log
  

注意!!!:上述命令中RGSM=sampleName参数，在对群体数据进行变异检测的时候，一定要保证每个样本名字不一样,否则后续生成的GVCF将无法合并，在实际操作中可以将RGSM设置成SRR accession或者直接设置为样本名称.

4、Remove duplicates

• 在这里重复reads是指来源于同一个DNA片段的reads，来源于建库过程（例如：建库过程中的PCR

  time java -XX:ParallelGCThreads=20 -jar $picard MarkDuplicates I=$opdir/$bn"Aligned.rg_added.sortedByCoord.out.bam" O=$opdir/$bn"Aligned.rg_added.dedupped.bam" CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT M=$opdir/output.metrics 2>$opdir/$bn.MarkDuplicates.log
  

5、SplitNCigarReads

• SplitNCigarReads是专为RNA-Seq数据设计的工具。主要功能见下图：

  SplitNCigarReads

  更多解释参见这里

  time java -XX:ParallelGCThreads=20 -d64 -jar $gatk -T SplitNCigarReads -allowPotentiallyMisencodedQuals -R $ref -I $opdir/$bn"Aligned.rg_added.dedupped.bam" -o $opdir/$bn"_split.bam" -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS 2>$opdir/$bn.SplitNCigarReads.log
  

6、Indel realignment (optional)

• 在分割完reads之后，对indel重新比对一定程度上可以挽回一些reads，但是对整体SNP结果影响不大，考虑到此步骤耗时（～30mins），因此在实际操作中这个步骤可以省略。

  #Create targets for indel realignment (2mins)
  echo -e "["$(date)"]\tCreating targets for indel realignment.."
  time java -d64 -jar $gatk -T RealignerTargetCreator -allowPotentiallyMisencodedQuals -R $ref -I $opdir/$bn"_split.bam" -o $opdir/$bn".intervals" -nt 20
  
  #Perform indel realignment (25mins)
  echo -e "["$(date)"]\tPerforming Indel Realignment.."
  time java -d64 -jar $gatk -T IndelRealigner -allowPotentiallyMisencodedQuals -R $ref -I $opdir/$bn"_split.bam" -targetIntervals $opdir/$bn".intervals" -o $opdir/$bn"_processed.bam"
  

7、Run HaplotypeCaller

• HaplotypeCaller既可以找SNP又可以找Indel

  #Run HaplotypeCaller
  echo -e "["$(date)"]\tRunning HaplotypeCaller.."
  time java -Xmx120g -XX:ParallelGCThreads=20 -d64 -jar $gatk -T HaplotypeCaller -allowPotentiallyMisencodedQuals -R $ref -I $opdir/$bn"_split.bam" -dontUseSoftClippedBases -stand_call_conf 20.0 -ERC GVCF -nct 20 -o $opdir/$bn".g.vcf" 2>$opdir/$bn.HaplotypeCaller.log
  

注意!!!: 上述命令中ERC参数用来指定输出文件格式为GVCF，所谓GVCF就是VCF的变种，即输出每一个位置的信息（不局限于变异）。

8、Combine GVCF

• 在所有样本都跑完2-7步骤之后，对于每一个样本生成的GVCF文件，以下命令将其合并成包含多个样本的GVCF

  java -jar $gatk -R $ref -T CombineGVCFs --variant sample1.g.vcf --variant sample2.g.vcf -o merge.g.vcf
  

注意!!!: 以上命令有几个样本就有几个variant参数, 同时一定要保证每个GVCF文件最后一列的行名不一样,否则将无法合并。

9、将GVCF转换为VCF

• 使用GenotypeGVCFs工具将GVCF转换为VCF

  java -jar $gatk -T GenotypeGVCFs -R $ref --variant merge.g.vcf -o output.vcf  -stand_call_conf 20
  

10、过滤

• 尽管GATK强烈推荐使用VQSR对检测到的变异进行质量控制。VQSR通过机器学习的方法利用多个不同的数据特征训练一个高斯混合模型对变异数据进行质控，然而不幸的是使用VQSR需要具备以下两个条件：

  • 需要一个精心准备的已知变异集，它将作为训练质控模型的真集。比如，对于人类来说，有Hapmap、OMNI，1000G和dbsnp等际性项目的数据可以作为高质量的已知变异集。GATK的bundle主要就是对这四个数据集做了精心的处理和选择，然后把它们作为VQSR时的真集位点。
  • 要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练的时候会因为可用的变异位点数目不足而无法进行。

  然而对于大多数物种都没有这样的数据集，这时候就需要执行硬筛选

  所谓硬筛选, 就是使用VQSR使用的数据指标进行“一刀切”:

  #使用SelectVariants，选出只有两个等位基因的SNP
  time java -jar GenomeAnalysisTK.jar -T SelectVariants  -R reference.fasta -V input.vcf -o snp.vcf -selectType SNP -restrictAllelesTo BIALLELIC
  time java -XX:ParallelGCThreads=20 -d64 -jar $gatk -T VariantFiltration -R $ref -V snp.vcf -window 35 -cluster 3 -filterName FS -filter "FS > 60.0" -filterName QD -filter "QD < 2.0" -o snp_filtered.vcf
  
  #使用SelectVariants，选出Indel
  time java -jar GenomeAnalysisTK.jar -T SelectVariants  -R reference.fasta -V input.vcf -o indel.vcf -selectType INDEL
  time java -XX:ParallelGCThreads=20 -d64 -jar $gatk -T VariantFiltration -R $ref -V indel.vcf -window 35 -cluster 3 -filterName FS -filter "FS > 60.0" -filterName QD -filter "QD < 2.0" -o indel_filtered.vcf 
  

  关于GATK执行硬筛选更详细的解释，可以查看这里。

VCF文件的压缩

• 由于生成的gVCF以及VCF文件都相对较大，尤其是在处理群体数据的时候。所以为了节省硬盘空间，可以使用工具bZIP对gVCF进行压缩。

#压缩
bgzip test.gVCF
#解压
tabix -p vcf test.gVCF.gz

• bZIP的安装可以参考这里