run:
    if len(input) == 1:
        shell("cp {input.run} {output.bam} ")
    else:
        shell("{samtools_app} merge --threads {threads} {output.bam} {input.run} {input.reseq} ")

shell:
    "{samtools_app} sort -@ {threads} {input.bam_file} -o {output.bam_sort}"

shell:
    "java -jar {picard_app} MarkDuplicates I={input.bam_file} "
    "O={output.bam_noDup} M={output.dup_file} REMOVE_DUPLICATES={params.rm_dup} "

shell:
    "{samtools_app} stats  {input.bam_file}  --threads {threads} > {output.txt_file}"

shell:
    "java -jar {picard_app} CollectWgsMetrics I={input.bam_file} "
    "O={output.collect} R={input.fasta_ref} "

shell:
    "{samtools_app} index {input.bam_file} {output.index_file}"

shell:
    "{gatk_app} BaseRecalibrator -I {input.bam_file} -R {input.fasta_ref} "
    " --known-sites {params.snps_1000} --known-sites {params.indels_Mills} "
    " --known-sites {params.hg38_indels} --known-sites {params.hg38_snps} "
    " --known-sites {params.snps_hapmap} -O {output.recal_table} "

shell:
    "{gatk_app} HaplotypeCaller "
    "-I {input.bam_file} -R {input.fasta_ref} -O {output.vcf_file} -ERC GVCF "

shell:
    "{samtools_app} view -hb -q 20 {input.bam_file} > {output.filter_bam}"

shell:
    "grep -w {params.chr} {input.interval_file} > {output.chr_file}"

run:
    iFiles = " ".join(["-V " + file for file in input.vcf_file ])
    shell("mkdir -p " + params.tempdir)
    shell("{gatk_app}  GenomicsDBImport " + iFiles + "  --genomicsdb-workspace-path {output.database} -L {input.interval_file} --tmp-dir {params.tempdir}")

shell:
    "{gatk_app} --java-options \"-Xmx4g\" GenotypeGVCFs -R {input.fasta_ref} -V gendb://{input.database} -O {output.cohort_vcf}"

run:
    iFiles = " ".join(["-I " + file for file in input.vcf_files ])
    shell("java -jar {picard_app} GatherVcfs " + iFiles + " -O {output.cohort} ")

shell:
    "{bcftools_app} norm -f {input.fasta_ref} -m -any -Oz -o {output.cohort} {input.vcf_file} "

shell:
    "{gatk_app} --java-options \"-Djava.io.tmpdir={params.tempdir} -Xms4G -Xmx4G -XX:ParallelGCThreads=2\" VariantRecalibrator \
      -tranche 100.0 -tranche 99.95 -tranche 99.9 -tranche 99.8 \
      -tranche 99.5 -tranche 99.0 -tranche 97.0 -tranche 96.0 \
      -tranche 95.0 -tranche 94.0 \
      -tranche 93.5 -tranche 93.0 -tranche 92.0 -tranche 91.0 -tranche 90.0 \
      -R {input.fasta_ref} \
      -V {input.cohort} \
      --resource:hapmap,known=false,training=true,truth=true,prior=15.0 {input.snps_hapmap} \
      --resource:omni,known=false,training=true,truth=false,prior=12.0 {input.snps_omni} \
      --resource:1000G,known=false,training=true,truth=false,prior=10.0 {input.snps_1000} \
      -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR  \
      -mode SNP -O {output.recal} --tranches-file {output.tranches} \
      --rscript-file {output.script} "

shell:
    "{gatk_app} --java-options \"-Djava.io.tmpdir={params.tempdir} -Xms4G -Xmx4G -XX:ParallelGCThreads=2\" VariantRecalibrator \
      -tranche 100.0 -tranche 99.95 -tranche 99.9 -tranche 99.8 \
      -tranche 99.5 -tranche 99.0 -tranche 97.0 -tranche 96.0 \
      -tranche 95.0 -tranche 94.0 -tranche 93.5 -tranche 93.0 \
      -tranche 92.0 -tranche 91.0 -tranche 90.0 \
      -R {input.fasta_ref} \
      -V {input.cohort} \
      --resource:mills,known=false,training=true,truth=true,prior=12.0 {input.indels_Mills} \
      --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 {input.hg38_snps} \
      -an QD -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \
      -mode INDEL -O {output.recal} --tranches-file {output.tranches} \
      --rscript-file {output.script} "

shell:
    "{gatk_app} --java-options \"-Djava.io.tmpdir={params.tempdir} \
      -Xms2G -Xmx2G -XX:ParallelGCThreads=2\" ApplyVQSR \
      -V {input.cohort} \
      --recal-file {input.recal} \
      -mode SNP \
      --tranches-file {input.tranches} \
      --truth-sensitivity-filter-level 99.8 \
      --create-output-variant-index true \
      -O {output.vcf}"

shell:
    "{gatk_app} --java-options \"-Djava.io.tmpdir={params.tempdir} \
      -Xms2G -Xmx2G -XX:ParallelGCThreads=2\" ApplyVQSR \
      -V {input.vcf} \
      -mode INDEL \
      --recal-file {input.recal} \
      --tranches-file {input.tranches} \
      --truth-sensitivity-filter-level 99.9 \
      --create-output-variant-index true \
      -O {output.vcf} "

shell:
    "{gatk_app} --java-options \"-Djava.io.tmpdir={params.tempdir} -Xms2G -Xmx2G -XX:ParallelGCThreads=2\" CalculateGenotypePosteriors "
    " -V {input.vcf} "
    " --supporting-callsets {input.wgs_1000GP} -O {output.vcf} "

shell:
    "{bcftools_app} +fill-tags {input.vcf} -Ou | {bcftools_app} view -t {params.chr_str} -Ou | {bcftools_app} annotate --set-id '%CHROM\_%POS\_%REF\_%ALT' -Oz -o {output.vcf} "

shell:
    "{bcftools_app} filter -S . -e 'FMT/DP<3 | FMT/GQ<20' -Oz -o {output.vcf} {input.vcf} "

shell:
    "{bcftools_app} filter -i 'FILTER=\"PASS\"' -Ou {input.vcf}| {bcftools_app} filter -i 'INFO/HWE > 1e-6' -Ou | {bcftools_app} filter -i 'F_MISSING < 0.15' -Oz -o {output.vcf}"

shell:
    "java -Xmx24g -jar {beagle_app} gt={input.vcf} nthreads={threads} out={params.prefix}  impute=false"

shell:
    "tabix -p vcf {input} "

shell:
    "{bcftools_app} +fill-tags  -Oz -o {output.vcf} {input.vcf}"

shell:
    "for chr in {{1..22},X}; do echo chr$\{chr\} $\{chr\} >> {params.chrs}; done"
    "{bcftools_app} annotate --rename-chrs {params.chrs} {input} -Oz -o {output} --set-id '%CHROM\:%POS\:%REF\:%ALT' "

shell:
    "{bcftools_app} filter -i 'INFO/MAF > {params.maf}' -Oz -o {output.vcf} {input.vcf} "

shell:
    "{bcftools_app} +dosage {input.vcf} -- -t GT > {output.table}"

shell:
    "{bcftools_app} query -f '%ID\n' {input.vcf}  > {output.table}"

shell:
    "Rscript {params.script} -g {input.genotype} -v {input.variants} -r {params.results}"