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}"""
        }
    }
}


Bpipe.run { setup + modifyBAMHeader + resortBAM + createBED + annotateVars }