最近有多個(gè)朋友問(wèn)我關(guān)于二代重測(cè)序數(shù)據(jù)的分析問(wèn)題（clean data 到vcf），我寫(xiě)了一個(gè)完整的步驟如下，若你是多個(gè)樣品進(jìn)行分析，最好是寫(xiě)
循環(huán)語(yǔ)句，這里我僅舉一個(gè)例子“sample1”作為樣品進(jìn)行腳本編寫(xiě)。
############? 中國(guó)科學(xué)院植物研究所? ?#####################
############? 孟慶霖? ?#########
##2021-6-10 參考基因組索引構(gòu)建
#第一步 bwa
cd /datapool/.../raw/ref/
bwa index ref.fa
#第二步samtools生成.fai文件
samtools faidx ref.fa > ref.fa.fai
#第三步 gatk 生成.dict文件
gatk CreateSequenceDictionary \
-R ref.fa \
-O ref.dict
#到此，參考基因組的index構(gòu)建完成。
###############################################
##下面開(kāi)始fq（clean data）文件到vcf文件的生成
#1生成bam
bwa mem -M -t 4 -R "@RG\tID:sample1\tSM:sample1\tPL:illumina" \
ref.fa sample1_1.fq.gz sample1_2.fq.gz > | \
samtools view -bS - > sample1.bam
#2對(duì)bam進(jìn)行sort
gatk --java-options "-Xmx15g" SortSam \
-I sample1.bam \
--SORT_ORDER coordinate \
-O sample1.sort.bam
##這一步只是說(shuō)明可以對(duì)sam進(jìn)行sort，如果你是bam，就直接sort就行
##gatk也可以對(duì)sam進(jìn)行sort
gatk --java-options "-Xmx15g" SortSam \
-I sample1.sam \
--SORT_ORDER coordinate \
-O sample1.sort.sam
#3對(duì)bam文件標(biāo)記pcr重復(fù)
##若是簡(jiǎn)化基因組測(cè)序，則不要進(jìn)行這一步的運(yùn)算
gatk --java-options "-Xmx15g" \
MarkDuplicates \
-I sample1.sort.bam \
--METRICS_FILE sample1.metrics.txt \
-O sample1.sort.md.bam
#4對(duì)3獲得bam進(jìn)行index
samtools index -@ 4 sample1.sort.md.bam
#5gatk call vcf變異
gatk --java-options "-Xmx15g" \
HaplotypeCaller \
-R ref.fa \
--emit-ref-confidence GVCF \
-I sample1.sort.md.bam \
-O sample1.g.vcf
gatk --java-options "-Xmx15g" \
GenotypeGVCFs \
-R ref.fa \
-V sample1.g.vcf \
-O sample1.raw.vcf
#6合并多個(gè)樣品的vcf文件，其-I順序即樣品在vcf中的排列順序
gatk --java-options "-Xmx15g" \
MergeVcfs \
-I sample1.raw.vcf \
-I sample2.raw.vcf \
-I sample3.raw.vcf \
...
-O all.raw.vcf
#7對(duì)vcf進(jìn)行硬過(guò)濾
gatk --java-options "-Xmx15g" \
VariantFiltration \
-R ref.fa \
-V all.raw.vcf \
--filter-expression \
"AC<2||QD<2.0||FS>60.0||MQ<40.0||MQRankSum < -12.5||ReadPosRankSum < -8.0 " \
--filter-name "ref.fa-name" \
-O all.raw.H.vcf? #H=硬過(guò)濾“hard”
grep "#" all.raw.H.vcf > all.raw.HP.vcf
grep "PASS" all.raw.H.vcf >> all.raw.HP.vcf
gatk --java-options "-Xmx15g" \
IndexFeatureFile \
-F all.raw.HP.vcf
#8將indel和snp分開(kāi)
gatk --java-options "-Xmx15g" \
-R ref.fa --STRICT false \
-I all.raw.HP.vcf \
--INDEL_OUTPUT all.raw.HP.indel.vcf \
--SNP_OUTPUT all.raw.HP.snp.vcf
最后：如果你運(yùn)行過(guò)程出錯(cuò)了=>首先，核對(duì)是否完全按照我編寫(xiě)的腳本進(jìn)行；其次，查看是否路徑設(shè)對(duì)了。
