analysis

#!/bin/sh

# software
VCFHACKS=~/bin/vcfhacks/
# EXPRESSION=~/projects/bridge_study/vcfs/work/expression.py
CSQ_QUERY=~/projects/bridge_study/vcfs/work/csq_query.py

# dirs
WORK=~/projects/bridge_study/vcfs/work/
AITMAN=${WORK}/aitman/
CAMBRIDGE=${WORK}/cambridge/

# files
VCF_AITMAN=${AITMAN}/vep.var.eds_09-06-16.aitman.ts99pt9.vcf_snp1pc_evs1pc_exac1pc.caddscored
VCF_CAMBRIDGE=${CAMBRIDGE}/vep.var._16-06-16.cambridge.ts99pt9.vcf_snp1pc_evs1pc_exac1pc.caddscored
EDS_BED=~/projects/bridge_study/updated_reportable_genes_GRCh37.bed
HELICAL_BED=~/projects/bridge_study/helical_domains.bed


######################################################
# sort_index()
#
# Sorts and Indexes the parsed vcf files.
#
######################################################

sort_index(){
    for vcf in $@; do
        vcf-sort $vcf | bgzip -c > ${vcf}.gz
        tabix -p vcf ${vcf}.gz
    done
}


####################################################
# AIMS:
# A. get all known pathogenic variants
#     A1: Variants in genes of interest.
#     A2: Filter on GT.Depth, GT.GQ and QUAL fields
#     A3. All known causitive variants
#     A4. LoF variants
#     A5. GLY-X-Y variants
# B. get all VUS variants
####################################################


for VCF in $VCF_AITMAN $VCF_CAMBRIDGE; do

    # A1: Filter for variants in genes of interest only
    getVariantsByLocation.pl -b $EDS_BED -i ${VCF}.vcf -o ${VCF}.eds_genes.vcf

    # A2: Convert samples with GT Depth < 20 or GQ > 20 to no calls and only keep variants
    #     with QUAL > 20
    filterGts.pl -d 20 -i ${VCF}.eds_genes.vcf -o ${VCF}.eds_genes.gt20.vcf
    getFunctionalVariants.pl \
          --input ${VCF}.eds_genes.gt20.vcf \
      --gq 20 \
      --var_qual 20 \
          --clinvar disable \
          --output ${VCF}.eds_genes.gt20.gq20.dp20.vcf

    # A3: ALL KNOWN CAUSATIVE VARIANTS
        # get Clinically Significant variants from ClinVarPathogenic and AS_CLNSIG fields
    echo "\n\nIdentifying clinically significant variants\n\n"
    getFunctionalVariants.pl \
          --input ${VCF}.eds_genes.gt20.gq20.dp20.vcf \
          --clinvar all \
          --output ${VCF}.eds_genes.gt20.gq20.dp20.clin_sig.vcf

    # A4: LOF VARIANTS
    echo "\n\nIdentifying LOF variants\n\n"
    getFunctionalVariants.pl \
      --input ${VCF}.eds_genes.gt20.gq20.dp20.vcf \
      --classes splice_acceptor_variant splice_donor_variant frameshift_variant stop_gained \
      --clinvar disable \
      --output ${VCF}.eds_genes.gt20.gq20.dp20.lof.vcf

    # A5: GLY-X-Y VARIANTS
    # filter for variants in helical domains
    # NOTE: getFunctionalVariants.pl --eval_filter won't work with "symbol eq 'COL1A2'"
    echo "\n\nIdentifying GLYXY variants\n\n"
    $CSQ_QUERY \
      --vcf ${VCF}.eds_genes.gt20.gq20.dp20.vcf \
      --helical_domains $HELICAL_BED \
      --out ${VCF}.eds_genes.gt20.gq20.dp20.helical_domain.vcf

    # A6: merge all filtered VCF files together
    # concatenate and sort all filtered VCF from A2-A4
    echo "\n\nConcatenating VCFs\n\n"
    DIRECTORY=`dirname $VCF`
    vcfs=`find $DIRECTORY -name "*.eds_genes.gt20.gq20.dp20*.vcf"`
    echo "\n\n\n$vcfs\n\n\n"
    sort_index $vcfs
    compressed_vcfs=`find $DIRECTORY -name "*.eds_genes.gt20.gq20.dp20*.vcf.gz" | xargs`
    vcf-concat $compressed_vcfs  > ${VCF}.eds_genes.gt20.gq20.dp20.pathogenic.vcf
    echo "\n\nSorting VCFs"
    sortVcf.pl -i ${VCF}.eds_genes.gt20.gq20.dp20.pathogenic.vcf -o ${VCF}.eds_genes.gt20.gq20.dp20.pathogenic.sorted.vcf
    # remove duplicate variants
    echo "\n\nRemoving Duplicates\n\n"
    awk -F '\t' '/^#/ {print;prev="";next;} {key=sprintf("%s\t%s\t%s",$1,$2,$3,$4,$5,$6);if(key==prev) next;print;prev=key;}' ${VCF}.eds_genes.gt20.gq20.dp20.pathogenic.sorted.vcf > ${VCF}.eds_genes.gt20.gq20.dp20.pathogenic.sorted.uniq.vcf
done


# Merge the aitman and cambridge pathogenic vcfs
AITMAN_VCF_PATHOGENIC=${VCF_AITMAN}.eds_genes.gt20.gq20.dp20.pathogenic.sorted.uniq.vcf
CAMBRIDGE_VCF_PATHOGENIC=${VCF_CAMBRIDGE}.eds_genes.gt20.gq20.dp20.pathogenic.sorted.uniq.vcf
sort_index $AITMAN_VCF_PATHOGENIC $CAMBRIDGE_VCF_PATHOGENIC
echo "\n\n\nMerging Pathogenic VCFS\n\n\n"
vcf-merge  ${AITMAN_VCF_PATHOGENIC}.gz ${CAMBRIDGE_VCF_PATHOGENIC}.gz > ${WORK}/eds_bridge.vep.var.ts99pt9.vcf_snp1pc_evs1pc_exac1pc.caddscored.eds_genes.gt20.gq20.dp20.pathogenic.sorted.uniq.combined.vcf