使用BWA+samtools+gatk4.0进行SNP calling

首先下载数据
eg：

$ wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz ###ref
###批量下载NGS数据，首先下载RunInfo
#1
prefetch SRR12424925
#2
wget -c https://sra-download.ncbi.nlm.nih.gov/traces/sra62/SRR/012133/SRR12424925(可以写个循环)#eg for a in list; do ***; done

01.align

首先定义相关用到的变量，然后使用baw创建索引之后进行map。

#!/usr/bin/bash
COV=/public/home/*/data/test
REF=$COV/ref/CM3.6.1_pseudomol.fa
SRR1=m1-115
SRR2=m1-116
BAM1=$COV/align/$SRR1.bam
BAM2=$COV/align/$SRR2.bam
R1_1=$COV/fq/${SRR1}.1.fq.gz
R1_2=$COV/fq/${SRR1}.2.fq.gz
R2_1=$COV/fq/${SRR2}.1.fq.gz
R2_2=$COV/fq/${SRR2}.2.fq.gz
VARI=$COV/variant
##bwa比对，samtools排序并构建索引，为了之后更快调用比对文件
mkdir -p $COV/align && cd $COV/align
bwa index $REF
samtools faidx $REF
bwa mem -R "@RG\tID:$SRR1/tSM:$SRR1\tLB:WGS\tPL:Illumina" $REF $R1_1 $R1_2 |samtools sort > $BAM1
bwa mem -R "@RG\tID:$SRR2/tSM:$SRR2\tLB:WGS\tPL:Illumina" $REF $R2_1 $R2_2 |samtools sort > $BAM2
samtools index $BAM1
samtools index $BAM2
mkdir -p $VARI

02.snp calling

对已经map好的文件使用gatk进行snp的calling，新版本gatk4.0已经整合相关软件的所用功能

##SNP calling
#构建.dict文件
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" CreateSequenceDictionary \
    -R $REF \
    -O $COV/ref/CM3.6.1_pseudomol.dict
#标记PCR重复reads
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" MarkDuplicates \
    -I $BAM1 -O ${SRR1}_MarkDup.bam \
    -M ${SRR1}.metrics
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" MarkDuplicates \
    -I $BAM2 -O ${SRR2}_MarkDup.bam \
    -M ${SRR2}.metrics
#用samtools对标记后的文件创建索引
samtools index ${SRR1}_MarkDup.bam
samtools index ${SRR2}_MarkDup.bam

#-bamout：将一整套经过gatk程序重新组装的单倍体基因型（haplotypes）输出到文件
#-stand-call-conf :低于这个数字的变异位点被忽略，可以设成标准30（默认是10）
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" HaplotypeCaller \
    -R $REF \
    -I ${SRR1}_MarkDup.bam -O ${VARI}/${SRR1}_gatk.gvcf \
    --emit-ref-confidence GVCF \
    -stand-call-conf 30
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" HaplotypeCaller \
    -R $REF \
    -I ${SRR2}_MarkDup.bam -O ${VARI}/${SRR2}_gatk.gvcf \
    --emit-ref-confidence GVCF \
    -stand-call-conf 30
#通常，GATK的HaplotypeCaller流程跑完后，得到的是单个样本的gVCF文件，即全基因组层面的变异信息，根据需要进行样本间gVCF文件的合并，并转化为VCF文件：
#合并多样本gVCF文件
gatk --java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true' GenomicsDBImport \
    -L chr00 \
    -V ${VARI}/${SRR1}_gatk.gvcf \
    -V ${VARI}/${SRR2}_gatk.gvcf \
    --genomicsdb-workspace-path $VARI/vcfdb
#gVCF转为VCF
gatk GenotypeGVCFs \
    -R $REF \
    -V gendb://$VARI/vcfdb \
    -O ${VARI}/HaplotypeCaller_gatk.vcf
gatk MergeVcfs \
    -I *.vcf.gz...
    -O *.vcf.gz

03.snp_filter

###过滤结果 
# 使用SelectVariants，选出SNP
gatk SelectVariants \
    -select-type SNP \
    -V chr00.vcf.gz \
    -O chr00.snp.vcf.gz

# 为SNP作硬过滤
gatk VariantFiltration \
    -V chr00.snp.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "Filter" \
    -O chr00.snp.filter.vcf.gz

# 使用SelectVariants，选出Indel
gatk SelectVariants \
    -select-type INDEL \
    -V chr00.vcf.gz \
    -O chr00.indel.vcf.gz

# 为Indel作过滤
gatk VariantFiltration \
    -V chr00.indel.vcf.gz \
    --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "Filter" \
    -O chr00.indel.filter.vcf.gz

# 重新合并过滤后的SNP和Indel
gatk MergeVcfs \
    -I chr00.snp.filter.vcf.gz \
    -I chr00.indel.filter.vcf.gz \
    -O chr00.filter.vcf.gz

# 删除无用中间文件
rm -f 

最终获得相关vcf.gz文件

04.snp annotation

first tool

使用VEP完成变异的注释

#!/usr/bin/bash
VEP=/your_path_to/ensembl-vep/vep
$VEP --fasta CM3.6.1_pseudomol.fa \
--vcf --merged --fork 10 --hgvs --force_overwrite --everything \
--offline --dir_cache /your_path_to/ensembl-vep/cachedir \
-i *.vcf.gz \
-o *.VEP.vcf.gz

second tool

使用snpEFF进行变异的注释

编辑.../snpEff/snpEff.config, 在# Databases & Genomes后增加新的物种信息记录

# melon
melon.genome : melon

在.../snpEff/data/新建一个melon文件夹，存入相关参考数据

mkdir -p .../snpEff/data/melon
$ ll data
sequences.fa.gz #参考基因组
genes.gff.gz #注释文件GFF2/3

build建库

java -jar snpEff/snpEff.jar build -gff3 -v melon

使用snpeff进行注释

#所有注释信息
java -jar snpEff.jar ann XXX XXX.vcf.gz > snpeff.vcf
#去除上游、下游、UTR、基因间区的注释结果
java -jar snpEff.jar ann -no-utr -no-downstream -no-upstream -no-intergenic XXX XXX.vcf.gz > snpeff.vcf
#eg
java -jar snpEff/snpEff.jar ann -no-utr -no-downstream -no-upstream -no-intergenic melon melon.vcf.gz > melon_eff.vcf
#