WGS--bwa+gatk+snpEff流程
2021-03-31 本文已影响0人
一只烟酒僧
一、step0.0.creat_dirlist.sh
#需要的文件
#测序数据、基因组+索引、已知的SNP位点文件
#!/bin/bash
#创建文件系统
cd $wk_dir
mkdir fastq #存储原始数据及各个阶段的数据
mkdir code #用于存储代码
mkdir index #存储索引数据
mkdir genome #存储基因组fasta数据
mkdir res #存储最终结果数据
mkdir temp #存储中间的temp数据
mkdir otherfile #存储其它过程中间数据
#将构建好的index和基因组分别放到index和genome中
二、step0.environment
#!/bin/bash
wk_dir=/media/whq/282A932A2A92F3D2/282A932A2A92F3D2/WHQ/xxx
code_dir=/media/whq/282A932A2A92F3D2/WHQ/xxx/code
fastq_dir=$wk_dir/fastq
index=$wk_dir/index/gatk_human_hg38.fasta
genome=$wk_dir/genome/Homo_sapiens_assembly38.fasta
known_site_dir=$wk_dir/index/gatk_site
dir_prefix=FDSW* #这里是因为我的文件都是FDSW开头的
fastq1_suffix=-1r_1.clean.fq.gz
fastq2_suffix=-1r_2.clean.fq.gz
三、step1.bwa+sort+markdup+BQSR
#!/bin/bash
source ./step0.environment
#BWA+sort
cd $wk_dir/fastq
for id in `ls -d $dir_prefix`
do
{
input1=$id/$id$fastq1_suffix;input2=$id/$id$fastq2_suffix;sample=$id;index=$index;bwa mem -t 15 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $index $input1 $input2|samtools sort -@ 10 -o $id/$id.sorted.bam -
} 2>$wk_dir/temp/$id'.bwa.sorted.errlog' &
done
wait
#mark_dup+BQSR
cd $wk_dir/fastq
for id in `ls -d $dir_prefix`
do
{
input=$id/$id'.sorted.bam';output=$id/$id'.sorted.markdup.bam';gatk MarkDuplicates -I $input -O $output -M $id/$id'.markdup.matrix'
} 2>$wk_dir/temp/$id'.markdup.errlog' &
done
wait
for id in `ls -d $dir_prefix`
do
{
input=$id/$id'.sorted.markdup.bam' ;output=$id/$id'.bqsr.table';ref=$genome;know_site1=$known_site_dir/1000G_phase1.snps.high_confidence.hg38.vcf.gz;know_site2=$known_site_dir/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz;know_site3=$known_site_dir/dbsnp_138.hg38.vcf.gz;gatk BaseRecalibrator -I $input -O $output -R $ref --known-sites $know_site1 --known-sites $know_site2 --known-sites $know_site3
} 2>$wk_dir/temp/$id'.BQSR.1.errlog' &
done
wait
for id in `ls -d $dir_prefix`
do
{
input=$id/$id'.sorted.markdup.bam';output=$id/$id'.sorted.markdup.BQSR.bam';ref=$genome;bqsr=$id/$id'.bqsr.table';gatk ApplyBQSR -R $ref -I $input -O $output -bqsr $bqsr
} 2>$wk_dir/temp/$id'.BQSR.2.errlog' &
done
四、step2.call_SNP+indel.sh
这里使用了HaplotypeCaller -ERC GVCF 先产生GVCF文件，然后CombineGVCFs+GenotypeGVCFs
#!/bin/bash
source ./step0.environment
#call SNP+INDEL
cd $wk_dir/fastq
for id in `ls -d $dir_prefix`
do
{
input=$id/$id'.sorted.markdup.BQSR.bam';output=$id/$id'_raw.gvcf';dbsnp=$known_site_dir/dbsnp_138.hg38.vcf.gz;ref=$genome;gatk HaplotypeCaller -R $ref -I $input -O $output --dbsnp $dbsnp -ERC GVCF
} 2>$wk_dir/temp/$id'.call.snpindel.errlog' &
done
wait
V1=FDSW210010721/FDSW210010721_raw.gvcf
V2=FDSW210010722/FDSW210010722_raw.gvcf
V3=FDSW210010723/FDSW210010723_raw.gvcf
V4=FDSW210010724/FDSW210010724_raw.gvcf
#合并gvcf
gatk CombineGVCFs \
-R $genome \
-V $V1 \
-V $V2 \
-V $V3 \
-V $V4 \
-O combined.g.vcf
#joint genotyping
gatk GenotypeGVCFs \
-R $genome \
-V combined.g.vcf \
-G StandardAnnotation \
-D $known_site_dir/dbsnp_138.hg38.vcf.gz \
-A Coverage \
-O raw_variants.vcf
五、step3.VQSR.sh
#!/bin/bash
source ./step0.environment
cd $wk_dir/fastq
VQSR_dir=/media/whq/282A932A2A92F3D2/282A932A2A92F3D2/WHQ/xxx/index/VQSR_file
# call SNP
time gatk VariantRecalibrator \
-R $genome \
-V raw_variants.vcf \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 $VQSR_dir/hapmap_3.3.hg38.vcf.gz \
-resource:omini,known=false,training=true,truth=false,prior=12.0 $VQSR_dir/1000G_omni2.5.hg38.vcf.gz \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 $VQSR_dir/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 $VQSR_dir/dbsnp_138.hg38.vcf.gz \
-an DP \
-an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode SNP -tranche 100.0 \
-tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
-O raw_variants.vcf.snp.recal \
--tranches-file raw_variants.vcf.snp.tranches \
--rscript-file raw_variants.vcf.snp.plots.R
time gatk ApplyVQSR \
-R $genome \
-V raw_variants.vcf \
--truth-sensitivity-filter-level 99.0 \
--tranches-file raw_variants.vcf.snp.tranches \
--recal-file raw_variants.vcf.snp.recal \
-mode SNP \
-O raw_variants.snp.VQSR.vcf.gz
#call INDEL 这里使用raw_variants.snp.VQSR.vcf.gz作为输入，可以直接在上一步的SNP输出文件中写入INDEL的filter信息！
gatk VariantRecalibrator \
-R $genome \
-V raw_variants.snp.VQSR.vcf.gz \
-resource:mills,known=true,training=true,truth=true,prior=12.0 $VQSR_dir/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode INDEL \
--max-gaussians 6 \
--rscript-file raw_variants.snp.indel.plots.R \
--tranches-file raw_variants.snp.indel.tranches \
-O raw_variants.snp.indel.recal
time gatk ApplyVQSR \
-R $genome \
-V raw_variants.snp.VQSR.vcf.gz \
--truth-sensitivity-filter-level 99.0 \
--tranches-file raw_variants.snp.indel.tranches \
--recal-file raw_variants.snp.indel.recal \
-mode INDEL \
-O raw_variants.snp.indel.VQSR.vcf.gz
六、step4.snpEff+gatk+ANNOVAR+VEP.sh
VEP:https://github.com/Ensembl/ensembl-vep#install
VEP配置及报错处理：https://www.cnblogs.com/afeiyuanda/p/13168166.html
SNPeff:
ANNOVAR:
#!/bin/bash
source ./step0.environment
cd $wk_dir/fastq
snpeff=/home/whq/miniconda3/share/snpeff-5.0-0/snpEff.jar
#snpEff annotation
input=raw_variants.snp.indel.VQSR.vcf.gz
output1=raw_variants.snp.indel.VQSR.vcf.snpeff.gz
java -Xmx8g -jar $snpeff GRCh38.99 $input >$output1
#gatk annotation
#ANNOVAR annotation
七、selectVariation.sh
#!/bin/bash
#注意！这里的vcf文件不能是压缩格式的！
source ./step0.environment
cd $wk_dir/fastq
gatk SelectVariants \
-R $genome \
-V raw_variants.snp.indel.VQSR.snpeff.vcf \
--select-type-to-include SNP \
-O raw_variants.snp.VQSR.snpeff.vcf
gatk SelectVariants \
-R $genome \
-V raw_variants.snp.indel.VQSR.snpeff.vcf \
--select-type-to-include INDEL \
-O raw_variants.indel.VQSR.snpeff.vcf
八、step6.statistics.density.snp.indel.sh
使用bcftools拆分VCF文件：https://www.cnblogs.com/muuyouzhi/p/9303389.html
#!/bin/bash
#注意！这里的vcf文件不能是压缩格式的！
source ./step0.environment
cd $wk_dir/fastq
#拆分总VCF为各样本独立的VCF
ls -d FDS*|while read id;do output=$id/${id}.snp.VQSR.snpeff.vcf;bcftools view -s $id raw_variants.snp.VQSR.snpeff.vcf -O v -o $output;done
ls -d FDS*|while read id;do output=$id/${id}.indel.VQSR.snpeff.vcf;bcftools view -s $id raw_variants.indel.VQSR.snpeff.vcf -O v -o $output;done
#计算SNP及INDEL在各样品中的密度 --SNPdensity是窗口宽度
ls -d FDS*|while read id;do input=$id/${id}.snp.VQSR.snpeff.vcf ;output=$id/${id}.snp.VQSR.snpeff.density;less $input |awk '$10!~/0\/0|\.\/\.|0\|0/{print $0}'|vcftools --vcf - --SNPdensity 25000 --out $output;done
ls -d FDS*|while read id;do input=$id/${id}.indel.VQSR.snpeff.vcf ;output=$id/${id}.indel.VQSR.snpeff.density;less $input |awk '$10!~/0\/0|\.\/\.|0\|0/{print $0}'|vcftools --vcf - --SNPdensity 25000 --out $output;done
九、step7.statistics.distribution.snp.indel.sh
