gistfile1.groovy

about title: "Variant filtering from amplicon data"

PICARD_HOME="/n/HSPH/local/share/java/picard/"
GATK_HOME="/n/HSPH/local/share/java/gatk/"
SNPEFF_HOME="/n/HSPH/local/share/java/snpeff/"
SNPEFF_CONF="/n/HSPH/local/share/java/snpeff/snpEff.config"
HG19="/n/scratch00/hsph/biodata/genomes/Hsapiens/hg19/seq/hg19.fa"
HG19CHR='../chrTargets.bed'
TARGETS="../exonTargetsSorted.bed"
DBSNP132="/n/scratch00/hsph/biodata/genomes/Hsapiens/hg19/variation/dbsnp_132.vcf"

/*
 * Store original VCF file for later use, change
 * pointer to correct BAM file
 */
setup = {
    RAWVCF="$input2"
    exec """sed -i 's/bowtie/bwa/g' ${input2}""" 
    forward input1
}

/* 
 * Add missing BAM header information
 */
@Filter("header")
modifyBAMHeader = {
    exec """java -Xmx2g -jar $PICARD_HOME/AddOrReplaceReadGroups.jar INPUT=${input} OUTPUT=${output} RGID=1002 RGCN=CCCB RGLB=LungCancerCCCB RGPL=ILLUMINA RGPU=NA RGSM=`basename ${input}`"""
}

/* 
 * Re-sort and index reads based on contig order of the reference genome
 */
@Filter("sorted")
resortBAM = {
    exec """java -Xmx2g -jar $PICARD_HOME/ReorderSam.jar INPUT=$input1 OUTPUT=${output} R=$HG19"""
    exec """samtools index ${output}"""
}

/*
 * Convert BAM file to BED, filter out reads outside target genome for TEQC
 */
@Transform("bed")
createBED = {
    exec """bedtools intersect -wa -bed -abam ${input} -b $HG19CHR > ${output}"""
}

/* 
 * Using the [GATK VariantAnnotator](http://www.broadinstitute.org/gsa/wiki/index.php/VariantAnnotator) 
 * and limiting annotation to already called VCF files (the `-L` option), 
 * stick to the standard annotation:
 */
annotateVars = {
    from(["bam", "vcf"]) {
        transform("vcf") {
            exec """java -Xmx2g -jar $GATK_HOME/GenomeAnalysisTK.jar -T VariantAnnotator -I ${input.bam} -R $HG19 --variant:VCF ${input.vcf} -L:VCF ${input.vcf} --dbsnp $DBSNP132 -A ReadDepthAndAllelicFractionBySample -o ${output}"""
        }
    }
}

/*
 * Filter VCF file based on coverage, mapping quality, haplotype frequency
 * Somewhat convoluted way to check that the first called variant has a 
 * allelic frequency >= 0.25
 */
@Filter("filter")
flagVariants = { 
    exec """java -Xmx2g -jar $GATK_HOME/GenomeAnalysisTK.jar -R $HG19 -T VariantFiltration -o ${output} --variant ${input} --filterExpression 'DP < 40' --filterName 'DepthOfCoverage' --filterExpression 'vc.getGenotype(0).getAttribute("FA").split(",").get(0) < 0.25' --filterName 'AlleleFrequency'"""
}

/* 
 * Only retain variants passing all filters
 */
@Filter("pass")
keepPass = {
     exec """java -Xmx2g -jar $GATK_HOME/GenomeAnalysisTK.jar -R $HG19 -T SelectVariants --variant ${input} -o ${output} -ef"""
}


/*
 * Remove variants outside of targeted regions
 */
@Filter("target")
keepTarget = {
    exec """bedtools intersect -wa -header -a ${input} -b $TARGETS > ${output}"""
}


/*
 * Add [SnpEff information](http://www.broadinstitute.org/gsa/wiki/index.php/Adding_Genomic_Annotations_Using_SnpEff_and_VariantAnnotator) (V2.05, GrCh37.64):
 */

@Filter("snpeff")
addSnpEff = {
    exec """java -Xmx4G -jar $SNPEFF_HOME/snpEff.jar eff -c $SNPEFF_CONF -v -onlyCoding true -i vcf -o vcf GRCh37.64 ${input} > ${output}"""
}


/* 
 * VAT doesn't handle the annotated files well. Use bedtools to filter out non-final variants
 * from the raw vcf file
 */
@Filter("VAT")
createVAT = {
    exec """bedtools intersect -wa -header -a $RAWVCF -b ${input} > ${output}"""
}

/*
 * Limit VAT files to SNPs
 */
@Filter("snp")
subsetVAT = {
     exec """java -Xmx2g -jar $GATK_HOME/GenomeAnalysisTK.jar -T SelectVariants -R $HG19 --variant ${input} --selectTypeToInclude 'SNP' --out ${output}"""
}



/*
 * Knitr report using TEQC to provide coverage information, QC
 */
reportTEQC = {
    from("bed") {
    exec """Rscript ../../scripts/teqcWrap.R -r `pwd`/`basename $input` -t ../exonTargetsSorted.bed"""
    }
}

Bpipe.run { setup + modifyBAMHeader + resortBAM + createBED + annotateVars + flagVariants + keepPass + keepTarget + addSnpEff + createVAT + subsetVAT + reportTEQC }