1 介紹
基因組重測(cè)序是對(duì)已知基因組序列的物種進(jìn)行不同個(gè)體的基因組測(cè)序城瞎，并在此基礎(chǔ)上對(duì)個(gè)體或群體進(jìn)行差異性分析。
基于全基因組重測(cè)序技術(shù)疾瓮，人們期望快速進(jìn)行資源普查篩選脖镀，尋找到大量遺傳變異，實(shí)現(xiàn)遺傳進(jìn)化分析及重要性狀候選基因的預(yù)測(cè)狼电。隨著測(cè)序成本降低和擁有參考基因組序列物種增多蜒灰，全基因組重測(cè)序成為動(dòng)植物育種和群體進(jìn)化研究迅速有效的方法。
? SNP 強(qiáng)調(diào)群體在這個(gè)位點(diǎn)有多樣性肩碟。
? Genotype 強(qiáng)調(diào)個(gè)體在這個(gè)位點(diǎn)是哪種基因型强窖。SNV屬于genotype的范疇，表示這個(gè)genotype 是個(gè)突變腾务，是第一無(wú)二的毕骡。
2 基因組SNP calling分析流程
? 首先利用bwa,samtools,picard對(duì)參考基因組ref.fa建立索引削饵，隨后使用基因組比對(duì)軟件bwa將fastq文件比對(duì)回參考基因組岩瘦，生成sam文件(這里使用例子數(shù)據(jù)由于小于10M，所以使用bwa的-is參數(shù)窿撬，當(dāng)參考基因組很大的時(shí)候启昧，使用bwtsw參數(shù))
? 使用picard對(duì)sam文件排序后生成bam文件(這里如果使用samtools的話會(huì)刪除重復(fù)位點(diǎn)，而picard只會(huì)標(biāo)記重復(fù)并不會(huì)刪除)劈伴。之后使用samtools對(duì)bam文件建立索引密末，利用已經(jīng)注釋過(guò)的indel.vcf文件對(duì)indel周圍進(jìn)行realignment，這個(gè)操作有兩個(gè)步驟跛璧，第一步生成需要進(jìn)行realignment的位置的信息严里，第二步 才是對(duì)這些位置進(jìn)行realignment，最后生成一個(gè)重排好的BAM追城。
? 接下來(lái)對(duì)mapping得到的bam文件進(jìn)行矯正刹碾，消除偏態(tài)，去除假陽(yáng)性座柱，這里主要針對(duì)測(cè)序質(zhì)量不是特別好的數(shù)據(jù)迷帜，可以根據(jù)fastqc結(jié)果對(duì)測(cè)序結(jié)果進(jìn)行解讀，如果比較好的話可以跳過(guò)色洞。
? 使用HaplotypeCaller發(fā)現(xiàn)變異并生成gvcf文件戏锹。隨后使用CombineGVCFs將不同文件整合，GenotypeGVCFs函數(shù)輸出vcf文件火诸。
? 下一步使用VariantRecalibrator校準(zhǔn)锦针，并根據(jù)可視化結(jié)果確定最佳參數(shù)，確定參數(shù)之后再運(yùn)行一次VariantRecalibrator得出結(jié)果
? 如何確定參數(shù)：(1) an注釋信息選擇:從圖中可以看出，這個(gè)模擬結(jié)果可以很好的將真實(shí)的變異位點(diǎn)和假陽(yáng)性變異位點(diǎn)分開(kāi)伞插，形成了明顯的界限割粮， 也就是說(shuō)，如果一個(gè)變異位點(diǎn)的這兩個(gè)注釋值媚污，只要有一個(gè)落在了界限之外舀瓢，就會(huì)被過(guò)濾掉。最主要的是要 看右邊兩個(gè)圖片耗美，只要能很好的區(qū)分開(kāi)novel和known以及filtered和retained就可以京髓。其實(shí)在如何選擇注釋值 存在一定得主觀性，因此商架，在做VariantRecalibrator時(shí)可以做兩次堰怨，第一次盡可能的多的選擇這些注釋值， 第一遍跑完之后蛇摸，選擇幾個(gè)區(qū)分好的备图，再做一次VariantRecalibrator，然后再做ApplyRecalibration赶袄。 (2) tranche值的設(shè)定:如果這個(gè)值設(shè)定的比較高的話揽涮，那么最后留下來(lái)的變異位點(diǎn)就會(huì)多，但同時(shí)假陽(yáng)性的位點(diǎn)也會(huì)相應(yīng)增加饿肺； run兩遍VariantRecalibrator蒋困，第一遍的時(shí)候多寫(xiě)幾個(gè)閾值，看哪個(gè)閾值結(jié)果好敬辣，選擇一個(gè)最好的閾值雪标，再run一遍 VariantRecalibrator。 區(qū)分標(biāo)準(zhǔn)： 1. 看結(jié)果中已知變異位點(diǎn)與新發(fā)現(xiàn)變異位點(diǎn)之間的比例溉跃，這個(gè)比例不要太大村刨，大多數(shù)新發(fā)現(xiàn)的變異都是假陽(yáng)性，如果太多的 話撰茎，可能假陽(yáng)性的比例就比較大嵌牺； 2. 看保留的變異數(shù)目，這個(gè)就要根據(jù)具體的需求進(jìn)行選擇了乾吻。 – 3. 看TI/TV值髓梅，對(duì)于人類全基因組，這個(gè)值應(yīng)該在2.15左右绎签，對(duì)于外顯子組枯饿，這個(gè)值應(yīng)該在3.2左右，不要太小或太大诡必，越接近 這個(gè)數(shù)值越好奢方，這個(gè)值如果太小搔扁，說(shuō)明可能存在比較多的假陽(yáng)性。 千人中所選擇的tranche值是99蟋字，僅供參考稿蹲。
? 最后使用annovar對(duì)vcf進(jìn)行注釋并可視化結(jié)果
3 代碼展示
3.1建立索引
對(duì)于基因組
#首先建立參考基因組索引
##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
對(duì)于轉(zhuǎn)錄組
###STAR建立索引
STAR -runThreadN 8 -runMode genomeGenerate \
-genomeDir ./star_index/ \
-genomeFastaFiles ./genome/chrX.fa \
-sjdbGTFfile ./genes/chrX.gtf
3.2 比對(duì)回參考基因組
基因組使用bwa
#bwa將測(cè)序文件比對(duì)回參考基因組
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
##這里引號(hào)中表示對(duì)SAM 頭文件添加標(biāo)簽
轉(zhuǎn)錄組使用STAR 2-pass，第一步比對(duì)
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ù)比對(duì)結(jié)果重新建立索引
STAR -runThreadN 8 -runMode genomeGenerate \
-genomeDir ./star_index_2pass/ \
-genomeFastaFiles ./genome/chrX.fa \
-sjdbFileChrStartEnd ./star_1pass/ERR188044SJ.out.tab
進(jìn)行第二次比對(duì)
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 對(duì)bam文件建立索引
#samtools對(duì)bam文件建立索引
/software/samtools-1.3/samtools index \
./reads/vars_20x/Sample12.bwa.sort.dedup.bam
由于 STAR 軟件使用的 MAPQ 標(biāo)準(zhǔn)與 GATK 不同鹊奖，而且比對(duì)時(shí)會(huì)有 reads 的片段落到內(nèi)含子區(qū)間苛聘，需要進(jìn)行一步 MAPQ 同步和 reads 剪切，使用 GATK 專為 RNA-seq 應(yīng)用開(kāi)發(fā)的工具 SplitNCigarReads 進(jìn)行操縱忠聚，它會(huì)將落在內(nèi)含子區(qū)間的 reads 片段直接切除设哗，并對(duì) MAPQ 進(jìn)行調(diào)整。DNA 測(cè)序的重測(cè)序應(yīng)用中也有序列比對(duì)軟件的 MAPQ 與 GATK 無(wú)法直接對(duì)接的情況两蟀，需要進(jìn)行調(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，對(duì)已知的 indel 區(qū)域附近的 reads 重新比對(duì)赂毯，可以稍微提高 indel 檢測(cè)的真陽(yáng)性率战虏，如果對(duì)于某些物種沒(méi)有已知indel區(qū)域的話官網(wǎng)推薦可以跳過(guò)這步直接生成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文件(選做)
主要是針對(duì)測(cè)序質(zhì)量不是特別高的數(shù)據(jù)
#bam文件矯正，選做遣鼓，看測(cè)序質(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 變異檢測(cè)啸盏，使用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對(duì)vcf文件進(jìn)行注釋重贺，可以使用convert2annovar.pl進(jìn)行格式轉(zhuǎn)換骑祟。比如將VCF和pileup格式的文件轉(zhuǎn)換為annovar的輸入格式。
#對(duì)gtf文件進(jìn)行注釋
/software/UCSC/gtfToGenePred -genePredExt test.gtf test_refGene.txt
#進(jìn)行注釋
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
