前言
GATK是由Broad institute開發(fā)的用于變異檢測(例如: SNP, Indel等)高通量測序數(shù)據(jù). 目前GATK已經(jīng)被廣泛應(yīng)用于人類以及哺乳動物從WGS (Whole Genome Sequencing) 數(shù)據(jù)中檢測遺傳變異。本文講展示如何使用GATK官方提供的Best Practice對群體RNA-Seq數(shù)據(jù)進(jìn)行變異檢測。
所需要的工具
詳細(xì)步驟及解釋
1、為參考基因組創(chuàng)建index
-
根據(jù)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 -
上述所用參數(shù)解釋
--runThreadN: 指定線程數(shù)目 --runMode：指定STAR子程序genomeGenerate --genomeDir: 指定參考基因組index路徑 --genomeFastaFiles: 參考基因組的fasta文件 --sjdbGTFfile: 參考基因組GTF注釋文件 --sjdbOverhang: 該參數(shù)用來指定用于拼接junctions數(shù)據(jù)庫時(shí)junctions周圍基因組序列的長度。理想情況下, 這個(gè)值應(yīng)該是reads長度減1。例如，對于Illumina雙端reads（100bp）,這個(gè)值應(yīng)該設(shè)置為99。但是大多數(shù)情況下，設(shè)置成100也能達(dá)到理想的效果。
2、STAR 2-pass alignment
-
2-pass alignment是指先使用指定參數(shù)比對到參考基因組，然后使用上一步生成的junction信息(SI.out.tab)重新創(chuàng)建基因組index，然后使進(jìn)行第二輪比對
#指定參考基因組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用于統(tǒng)計(jì)該行命令所需要的時(shí)間 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" -
上述所用參數(shù)解釋
--outFileNamePrefix: 輸出文件名前綴 --outSAMtype: 指定輸出SAM類型，該處使用排序過后的BAM (BAM SortedByCoordinate) --outFilterMultimapNmax: 指定輸出唯一匹配的reads, 1代表uniquely mapped reads --outSAMstrandField: 對于非鏈特異性的RNA-Seq，Cufflinks/Cuffdiff需要指定此參數(shù) --genomeDir: 參考基因組index文件 --runThreadN: 指定線程數(shù)目 --readFilesIn: FASTQ文件 --twopassMode: 指定2-pass比對模式,即先使用上述參數(shù)比對到參考基因組，然后使用上一步生成的junction信息(SI.out.tab)重新創(chuàng)建基因組index，然后使用上述參數(shù)進(jìn)行第二輪比對
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參數(shù)，在對群體數(shù)據(jù)進(jìn)行變異檢測的時(shí)候，一定要保證每個(gè)樣本名字不一樣,否則后續(xù)生成的GVCF將無法合并，在實(shí)際操作中可以將RGSM設(shè)置成SRR accession或者直接設(shè)置為樣本名稱.
4、Remove duplicates
-
在這里重復(fù)reads是指來源于同一個(gè)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數(shù)據(jù)設(shè)計(jì)的工具。主要功能見下圖：
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結(jié)果影響不大，考慮到此步驟耗時(shí)（～30mins），因此在實(shí)際操作中這個(gè)步驟可以省略。
#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參數(shù)用來指定輸出文件格式為GVCF，所謂GVCF就是VCF的變種，即輸出每一個(gè)位置的信息（不局限于變異）。
8、Combine GVCF
- 在所有樣本都跑完2-7步驟之后，對于每一個(gè)樣本生成的GVCF文件，以下命令將其合并成包含多個(gè)樣本的GVCF
java -jar $gatk -R $ref -T CombineGVCFs --variant sample1.g.vcf --variant sample2.g.vcf -o merge.g.vcf
注意!!!: 以上命令有幾個(gè)樣本就有幾個(gè)variant參數(shù), 同時(shí)一定要保證每個(gè)GVCF文件最后一列的行名不一樣,否則將無法合并。
9、將GVCF轉(zhuǎn)換為VCF
- 使用GenotypeGVCFs工具將GVCF轉(zhuǎn)換為VCF
java -jar $gatk -T GenotypeGVCFs -R $ref --variant merge.g.vcf -o output.vcf -stand_call_conf 20
10、過濾
-
盡管GATK強(qiáng)烈推薦使用VQSR對檢測到的變異進(jìn)行質(zhì)量控制。VQSR通過機(jī)器學(xué)習(xí)的方法利用多個(gè)不同的數(shù)據(jù)特征訓(xùn)練一個(gè)高斯混合模型對變異數(shù)據(jù)進(jìn)行質(zhì)控，然而不幸的是使用VQSR需要具備以下兩個(gè)條件：
- 需要一個(gè)精心準(zhǔn)備的已知變異集，它將作為訓(xùn)練質(zhì)控模型的真集。比如，對于人類來說，有Hapmap、OMNI，1000G和dbsnp等際性項(xiàng)目的數(shù)據(jù)可以作為高質(zhì)量的已知變異集。GATK的bundle主要就是對這四個(gè)數(shù)據(jù)集做了精心的處理和選擇，然后把它們作為VQSR時(shí)的真集位點(diǎn)。
- 要求新檢測的結(jié)果中有足夠多的變異，不然VQSR在進(jìn)行模型訓(xùn)練的時(shí)候會因?yàn)榭捎玫淖儺愇稽c(diǎn)數(shù)目不足而無法進(jìn)行。
然而對于大多數(shù)物種都沒有這樣的數(shù)據(jù)集，這時(shí)候就需要執(zhí)行硬篩選
所謂硬篩選, 就是使用VQSR使用的數(shù)據(jù)指標(biāo)進(jìn)行“一刀切”:
#使用SelectVariants，選出只有兩個(gè)等位基因的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關(guān)于GATK執(zhí)行硬篩選更詳細(xì)的解釋，可以查看這里。
VCF文件的壓縮
- 由于生成的gVCF以及VCF文件都相對較大，尤其是在處理群體數(shù)據(jù)的時(shí)候。所以為了節(jié)省硬盤空間，可以使用工具bZIP對gVCF進(jìn)行壓縮。
#壓縮
bgzip test.gVCF
#解壓
tabix -p vcf test.gVCF.gz
- bZIP的安裝可以參考這里
