最近有多个朋友问我关于二代重测序数据的分析问题（clean data 到vcf），我写了一个完整的步骤如下，若你是多个样品进行分析，最好是写
循环语句，这里我仅举一个例子“sample1”作为样品进行脚本编写。
############ 中国科学院植物研究所 #####################
############ 孟庆霖 #########
##2021-6-10 参考基因组索引构建
#第一步 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构建完成。
###############################################
##下面开始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对bam进行sort
gatk --java-options "-Xmx15g" SortSam \
-I sample1.bam \
--SORT_ORDER coordinate \
-O sample1.sort.bam
##这一步只是说明可以对sam进行sort，如果你是bam，就直接sort就行
##gatk也可以对sam进行sort
gatk --java-options "-Xmx15g" SortSam \
-I sample1.sam \
--SORT_ORDER coordinate \
-O sample1.sort.sam
#3对bam文件标记pcr重复
##若是简化基因组测序，则不要进行这一步的运算
gatk --java-options "-Xmx15g" \
MarkDuplicates \
-I sample1.sort.bam \
--METRICS_FILE sample1.metrics.txt \
-O sample1.sort.md.bam
#4对3获得bam进行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合并多个样品的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对vcf进行硬过滤
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=硬过滤“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分开
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
最后：如果你运行过程出错了=>首先，核对是否完全按照我编写的脚本进行；其次，查看是否路径设对了。
网友评论