基因組重測序SNP_calling

1 介紹

基因組重測序是對已知基因組序列的物種進行不同個體的基因組測序，并在此基礎(chǔ)上對個體或群體進行差異性分析。

基于全基因組重測序技術(shù)，人們期望快速進行資源普查篩選，尋找到大量遺傳變異，實現(xiàn)遺傳進化分析及重要性狀候選基因的預(yù)測。隨著測序成本降低和擁有參考基因組序列物種增多，全基因組重測序成為動植物育種和群體進化研究迅速有效的方法。

? SNP 強調(diào)群體在這個位點有多樣性。

? Genotype 強調(diào)個體在這個位點是哪種基因型。SNV屬于genotype的范疇，表示這個genotype 是個突變，是第一無二的。

gatk分析流程
gatk相關(guān)參數(shù)

2 基因組SNP calling分析流程

? 首先利用bwa,samtools,picard對參考基因組ref.fa建立索引，隨后使用基因組比對軟件bwa將fastq文件比對回參考基因組，生成sam文件(這里使用例子數(shù)據(jù)由于小于10M，所以使用bwa的-is參數(shù)，當(dāng)參考基因組很大的時候，使用bwtsw參數(shù))

? 使用picard對sam文件排序后生成bam文件(這里如果使用samtools的話會刪除重復(fù)位點，而picard只會標(biāo)記重復(fù)并不會刪除)。之后使用samtools對bam文件建立索引，利用已經(jīng)注釋過的indel.vcf文件對indel周圍進行realignment，這個操作有兩個步驟，第一步生成需要進行realignment的位置的信息，第二步 才是對這些位置進行realignment，最后生成一個重排好的BAM。

? 接下來對mapping得到的bam文件進行矯正，消除偏態(tài)，去除假陽性，這里主要針對測序質(zhì)量不是特別好的數(shù)據(jù)，可以根據(jù)fastqc結(jié)果對測序結(jié)果進行解讀，如果比較好的話可以跳過。

? 使用HaplotypeCaller發(fā)現(xiàn)變異并生成gvcf文件。隨后使用CombineGVCFs將不同文件整合，GenotypeGVCFs函數(shù)輸出vcf文件。

? 下一步使用VariantRecalibrator校準(zhǔn)，并根據(jù)可視化結(jié)果確定最佳參數(shù)，確定參數(shù)之后再運行一次VariantRecalibrator得出結(jié)果

VQSR(Variant Quality Score Recalibration)

? 如何確定參數(shù)：(1) an注釋信息選擇:從圖中可以看出，這個模擬結(jié)果可以很好的將真實的變異位點和假陽性變異位點分開，形成了明顯的界限， 也就是說，如果一個變異位點的這兩個注釋值，只要有一個落在了界限之外，就會被過濾掉。最主要的是要 看右邊兩個圖片，只要能很好的區(qū)分開novel和known以及filtered和retained就可以。其實在如何選擇注釋值 存在一定得主觀性，因此，在做VariantRecalibrator時可以做兩次，第一次盡可能的多的選擇這些注釋值， 第一遍跑完之后，選擇幾個區(qū)分好的，再做一次VariantRecalibrator，然后再做ApplyRecalibration。 (2) tranche值的設(shè)定:如果這個值設(shè)定的比較高的話，那么最后留下來的變異位點就會多，但同時假陽性的位點也會相應(yīng)增加； run兩遍VariantRecalibrator，第一遍的時候多寫幾個閾值，看哪個閾值結(jié)果好，選擇一個最好的閾值，再run一遍 VariantRecalibrator。 區(qū)分標(biāo)準(zhǔn)： 1. 看結(jié)果中已知變異位點與新發(fā)現(xiàn)變異位點之間的比例，這個比例不要太大，大多數(shù)新發(fā)現(xiàn)的變異都是假陽性，如果太多的 話，可能假陽性的比例就比較大； 2. 看保留的變異數(shù)目，這個就要根據(jù)具體的需求進行選擇了。 – 3. 看TI/TV值，對于人類全基因組，這個值應(yīng)該在2.15左右，對于外顯子組，這個值應(yīng)該在3.2左右，不要太小或太大，越接近 這個數(shù)值越好，這個值如果太小，說明可能存在比較多的假陽性。 千人中所選擇的tranche值是99，僅供參考。

? 最后使用annovar對vcf進行注釋并可視化結(jié)果

3 代碼展示

3.1建立索引

對于基因組

#首先建立參考基因組索引
##bwa
bwa index -a is ./ref/TAIR9_chr1.random.fasta
##samtools
samtools faidx ./ref/TAIR9_chr1.random.fasta
##picard
java   -jar   /software/picard-tools-1.119/CreateSequenceDictionary.jar \ 
REFERENCE=./ref/TAIR9_chr1.random.fasta\ 
OUTPUT=./ref/TAIR9_chr1.random.dict

對于轉(zhuǎn)錄組

###STAR建立索引
STAR -runThreadN 8 -runMode genomeGenerate \
        -genomeDir ./star_index/ \
        -genomeFastaFiles ./genome/chrX.fa \
        -sjdbGTFfile ./genes/chrX.gtf

3.2 比對回參考基因組

基因組使用bwa

#bwa將測序文件比對回參考基因組
bwa mem -t 4 -M -R "@RG\tID:Sample12\tLB:Sample12\tPL:Illumina\tPU:Sample12\tSM:Sample12" \ 
./ref/TAIR9_chr1.random.fasta \ 
./reads/vars_20x/Sample12_1.fq\ 
./reads/vars_20x/Sample12_2.fq\ 
>./reads/vars_20x/Sample12.bwa.sam
##這里引號中表示對SAM 頭文件添加標(biāo)簽

轉(zhuǎn)錄組使用STAR 2-pass，第一步比對

STAR -runThreadN 8 -genomeDir ./star_index/ \
        -readFilesIn ./samples/ERR188044_chrX_1.fastq.gz ./samples/ERR188044_chrX_2.fastq.gz \
        -readFilesCommand zcat \
        -outFileNamePrefix ./star_1pass/ERR188044

根據(jù)比對結(jié)果重新建立索引

STAR -runThreadN 8 -runMode genomeGenerate \
        -genomeDir ./star_index_2pass/ \
        -genomeFastaFiles ./genome/chrX.fa \
        -sjdbFileChrStartEnd ./star_1pass/ERR188044SJ.out.tab

進行第二次比對

STAR -runThreadN 8 -genomeDir ./star_index_2pass/ \
        -readFilesIn ./samples/ERR188044_chrX_1.fastq.gz ./samples/ERR188044_chrX_2.fastq.gz \
        -readFilesCommand zcat \
        -outFileNamePrefix ./star_2pass/ERR188044

3.3 sort并轉(zhuǎn)換為bam文件

基因組使用SortSam.jar

#picard排序并轉(zhuǎn)換為bam文件
java   -jar   /software/picard-tools-1.119/SortSam.jar \ 
VALIDATION_STRINGENCY=SILENT \ 
INPUT=./reads/vars_20x/Sample12.bwa.sam \ 
OUTPUT=./reads/vars_20x/Sample12.bwa.sort.bam \ 
SORT_ORDER=coordinate

轉(zhuǎn)錄組使用AddOrReplaceReadGroups.jar

java -jar picard.jar AddOrReplaceReadGroups \
        I=./star_2pass/ERR188044Aligned.out.sam \
        O=./star_2pass/ERR188044_rg_added_sorted.bam \
        SO=coordinate \
        RGID=ERR188044 \
        RGLB=rna \
        RGPL=illumina \
        RGPU=hiseq \
        RGSM=ERR188044

3.4 標(biāo)記重復(fù)

#picard標(biāo)記重復(fù)
java -jar  /software/picard-tools-1.119/MarkDuplicates.jar \ 
MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 \ 
INPUT=./reads/vars_20x/Sample12.bwa.sort.bam \ 
OUTPUT=./reads/vars_20x/Sample12.bwa.sort.dedup.bam \ 
METRICS_FILE=./reads/vars_20x/Sample12.bwa.sort.dedup.metrics \ 
VALIDATION_STRINGENCY=LENIENT

3.5 對bam文件建立索引

#samtools對bam文件建立索引
/software/samtools-1.3/samtools index \ 
./reads/vars_20x/Sample12.bwa.sort.dedup.bam

由于 STAR 軟件使用的 MAPQ 標(biāo)準(zhǔn)與 GATK 不同，而且比對時會有 reads 的片段落到內(nèi)含子區(qū)間，需要進行一步 MAPQ 同步和 reads 剪切，使用 GATK 專為 RNA-seq 應(yīng)用開發(fā)的工具 SplitNCigarReads 進行操縱，它會將落在內(nèi)含子區(qū)間的 reads 片段直接切除，并對 MAPQ 進行調(diào)整。DNA 測序的重測序應(yīng)用中也有序列比對軟件的 MAPQ 與 GATK 無法直接對接的情況，需要進行調(diào)整。

java -jar GenomeAnalysisTK.jar -T SplitNCigarReads \
        -R ./genome/chrX.fa \
        -I ./star_2pass/ERR188044_dedup.bam \
        -o ./star_2pass/ERR188044_dedup_split.bam \
        -rf ReassignOneMappingQuality \
        -RMQF 255 \
        -RMQT 60 \
        -U ALLOW_N_CIGAR_READS

3.6 重排bam文件(選做)

之后就是可選的 Indel Realignment，對已知的 indel 區(qū)域附近的 reads 重新比對，可以稍微提高 indel 檢測的真陽性率，如果對于某些物種沒有已知indel區(qū)域的話官網(wǎng)推薦可以跳過這步直接生成vcf文件，之后使用第一次生成的vcf文件作為已知indel再走一遍流程。

#重排
##第一步首先生成位置信息
java   -jar  /software/GenomeAnalysisTK.jar \ 
-R   ./ref/TAIR9_chr1.random.fasta \ 
-T   RealignerTargetCreator \ 
-o   ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.intervals \ 
-I   ./reads/vars_20x/Sample12.bwa.sort.dedup.bam \ 
-known   ./variants/TAIR9_chr1.random.indel.vcf

##根據(jù)上一步位置信息，生成重排好的bam文件
java     -jar    /software/GenomeAnalysisTK.jar \ 
-R   ./ref/TAIR9_chr1.random.fasta \ 
-T    IndelRealigner \ 
-targetIntervals ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.intervals \ 
-o    ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.bam \ 
-I    ./reads/vars_20x/Sample12.bwa.sort.dedup.bam \ 
-known   ./variants/TAIR9_chr1.random.indel.vcf \ 
--consensusDeterminationModel KNOWNS_ONLY \ 
-LOD 0.4

3.7 BQSR矯正bam文件(選做)

主要是針對測序質(zhì)量不是特別高的數(shù)據(jù)

#bam文件矯正，選做，看測序質(zhì)量
java   -jar   /software/GenomeAnalysisTK.jar \ 
-T   BaseRecalibrator \ 
-R  ./ref/TAIR9_chr1.random.fasta \ 
-I   ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.bam \ 
-knownSites ./variants/TAIR9_chr1.random.vcf \ 
-o   ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.recal.table

java   -jar   /software/GenomeAnalysisTK.jar \ 
-T   PrintReads \ 
-R   ./ref/TAIR9_chr1.random.fasta \ 
-I   ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.bam \ 
-BQSR   ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.recal.table \ 
-o   ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.recal.bam

3.8 變異檢測，使用HaplotypeCaller 

#call variation發(fā)現(xiàn)變異并生成gvcf文件
java     -jar     /software/GenomeAnalysisTK.jar \ 
-R ./ref/TAIR9_chr1.random.fasta \
-T HaplotypeCaller \
-I  ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.recal.bam \
-nct 2 \ 
--emitRefConfidence GVCF \ 
--variant_index_type LINEAR \ 
--variant_index_parameter 128000 \ 
-o    ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.recal.hc.gvcf

3.9 整合gvcf文件并轉(zhuǎn)化為vcf文件

#整合gvcf文件
java    -jar    /software/GenomeAnalysisTK.jar \ 
-R     ./ref/TAIR9_chr1.random.fasta \ 
-T     CombineGVCFs \ 
-o     ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.gvcf \ 
-V     ./reads/vars_20x/Sample12.bwa.sort.dedup.realn_known.recal.hc.gvcf \ 
-V     ./reads/vars_20x/Sample13.bwa.sort.dedup.realn_known.recal.hc.gvcf

#gvcf文件轉(zhuǎn)化為vcf文件
java    -jar    /software/GenomeAnalysisTK.jar \ 
-R    ./ref/TAIR9_chr1.random.fasta \ 
-T    GenotypeGVCFs \ 
-nt 8 \ 
-V   ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.gvcf \ 
-o    ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.vcf

3.10 VariantRecalibrator校準(zhǔn)SNP

java   -jar   /software/GenomeAnalysisTK.jar \ 
-T    VariantRecalibrator \ 
-R   ./ref/TAIR9_chr1.random.fasta \ 
-input   ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.vcf \ 
-resource:original,known=false,training=true,truth=true,prior=15.0 \ 
./variants/TAIR9_chr1.random.snp.vcf \ 
-an DP  -an QD  -an FS  -an MQRankSum -an ReadPosRankSum \ 
-mode  SNP \ 
--maxGaussians 2  --maxNegativeGaussians 1  --minNumBadVariants 1000 \ 
-tranche  100.0   -tranche  99.9  -tranche  99.0  -tranche  90.0 \ 
-recalFile 
./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.SNP.recal \ 
-tranchesFile
./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.SNP.tranches\ 
-rscriptFile ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.SNP_plots.R

3.11 最終輸出

java    -jar   /software/GenomeAnalysisTK.jar \ 
-T    ApplyRecalibration \ 
-R    ./ref/TAIR9_chr1.random.fasta \ 
-input    ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.vcf \ 
-mode    SNP \ 
--ts_filter_level 99.0 \ 
-recalFile ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.SNP.recal\ 
-tranchesFile ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.SNP.tranches\ 
-o    ./variants/Samples/Samples.bwa.sort.dedup.realn_known.recal.hc.recal_SNP.vcf

3.12 后續(xù)分析

隨后可以使用annovar對vcf文件進行注釋，可以使用convert2annovar.pl進行格式轉(zhuǎn)換。比如將VCF和pileup格式的文件轉(zhuǎn)換為annovar的輸入格式。

#對gtf文件進行注釋
/software/UCSC/gtfToGenePred -genePredExt test.gtf test_refGene.txt
#進行注釋
convert2annovar.pl -format pileup variant.pileup -outfile variant.query
convert2annovar.pl -format vcf4 test.vcf > test.avinput
#輸出注釋文件
annotate_variation.pl -hgvs -buildver test_refGene.txt --geneanno test.avinput refdir/   
-outfile annotate.snp
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容