WGS--bwa+gatk+snpEff流程

作者: 一只烟酒僧 | 来源:发表于2021-03-31 20:41 被阅读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

网友评论

本文标题：WGS--bwa+gatk+snpEff流程

本文链接：https://www.haomeiwen.com/subject/depbqltx.html

延伸阅读

深度阅读

您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

WGS--bwa+gatk+snpEff流程

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

重测序分析

生信分析工具包