Merge pull request #6 from harvardinformatics/bam_resources

Bam resources

Author: tsackton

Snakefile_bam2vcf_gatk

@@ -18,7 +18,7 @@ bamDir = config["bamsForGatk"]
 gvcfDir = gatkDir + config["gvcfDir"]
 dbDir = gatkDir + config["dbDir"]
 vcfDir = gatkDir + config["vcfDir_gatk"]
-intDir = config["intDir"]
+intDir = config["intDir"] 
 maxIntervalLen = int(config["maxIntervalLen"])
 maxBpPerList = int(config["maxBpPerList"])
 maxIntervalsPerList = int(config["maxIntervalsPerList"])
@@ -27,7 +27,6 @@ minNmer = config["minNmer"]
 refBaseName = helperFun.getRefBaseName(config["ref"])
 
 # grab all samples for R1 to get list of names, no need to look at R2 which should have identical names
-#SAMPLES = ["ERR1013163"]
 SAMPLES = helperFun.getBamSampleNames(bamDir, bam_suffix)
 
 if not os.path.isdir(config["gatkDir"]):

Snakefile_fastq2bam

@@ -11,7 +11,6 @@ res_config = yaml.load(open("resources.yaml"))
 # rename variables from config file for clarity downstream
 fastq_suffix1 = config["fastq_suffix1"]
 fastq_suffix2 = config["fastq_suffix2"]
-fastqDir = config["fastqDir"]
 
 # this is where Snakemake output will go, specify with baseDir in config.yml
 fastqFilterDir = config["fastq2bamDir"] + config["fastqFilterDir"]

Snakefile_intervals

@@ -29,7 +29,7 @@ SAMPLES = helperFun.getBamSampleNames(bamDir, bam_suffix)
 if not os.path.isdir(intDir + "gatkLists"):
     os.system("mkdir -p " + intDir + "gatkLists")
 
-_, sample_dict, _ = helperFun.create_sample_dict("samples.csv")
+_, sample_dict, _ = helperFun.create_sample_dict(config["samples"])
 
 GENOMES = {sample_dict[k]["refGenome"]: sample_dict[k]["Organism"] for k in sample_dict.keys()}
 print(GENOMES)

config.yaml

@@ -3,11 +3,12 @@
 ##############################
 
 samples: "samples.csv"         # name of the sample metadata CSV 
-ref: "data/BHduck/genome/MU014702.1.fa"   # location of the reference genome
+ref: "data/BHduck/genome/MU014702.1.fa"   # location of the reference genome 
 spp: "hetAtr"                             # species name/code for the final VCF output file name
+genome: "MU014702.1"                      # name of the reference genome without the path or file extension
 
 # If using the fastq -> BAM workflow change these, otherwise ignore them
-fastqDir: "data/BHduck/fastq/"      # location of fastq files; must be followed by a "/"
+#fastqDir: "data/BHduck/fastq/"      # location of fastq files; must be followed by a "/"
 fastq_suffix1: "_1.fastq.gz"            # the suffix for the forward reads that follows all the sample names, e.g. "_1.fastq.gz" for sample "sampleName_1.fastq.gz"
 fastq_suffix2: "_2.fastq.gz"            # the suffix for the reverse reads that follows all the sample names, e.g. "_2.fastq.gz" for sample "sampleName_2.fastq.gz"
 

resources.yaml

@@ -28,6 +28,8 @@ dedup:
 # calculate BAM summaries with samtools and picard
 bam_sumstats: 
     mem: 9000 
+merge_bams:
+    mem: 9000
 
 ###
 # BAM -> VCF workflow

rules/bam2vcf_gatk.smk

@@ -148,7 +148,7 @@ rule gatherVcfs:
 rule vcftools:
     input:
         vcf = config["gatkDir"] + config['spp'] + "_final.vcf.gz",
-        int = intDir + "intervals_fb.bed"
+        int = intDir + config["genome"] + "_intervals_fb.bed"
     output: 
         missing = gatkDir + "missing_data_per_ind.txt",
         SNPsPerInt = gatkDir + "SNP_per_interval.txt"

rules/fastq2bam.smk

@@ -120,8 +120,11 @@ rule merge_bams:
         bai = bamDir + "{sample}_sorted.bam.bai"
     conda:
         "../envs/fastq2bam.yml"
+    resources:
+        mem_mb = lambda wildcards, attempt: attempt * res_config['merge_bams']['mem']
     shell:
         "samtools merge {output.bam} {input} && samtools index {output.bam}"
+
 rule dedup:
     input: 
         bamDir + "{sample}_sorted.bam",