WGS--bwa+gatk+snpEff流程

一舍沙、step0.0.creat_dirlist.sh

#需要的文件
#測序數(shù)據(jù)详羡、基因組+索引耕渴、已知的SNP位點文件
#!/bin/bash
#創(chuàng)建文件系統(tǒng)
cd $wk_dir
mkdir fastq  #存儲原始數(shù)據(jù)及各個階段的數(shù)據(jù)
mkdir code #用于存儲代碼
mkdir index  #存儲索引數(shù)據(jù)
mkdir genome  #存儲基因組fasta數(shù)據(jù)
mkdir res  #存儲最終結(jié)果數(shù)據(jù)
mkdir temp  #存儲中間的temp數(shù)據(jù)
mkdir otherfile  #存儲其它過程中間數(shù)據(jù)
#將構(gòu)建好的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 先產(chǎn)生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