-We recommend installing snakemake into a conda environment called `snakemake`.
-## Setup
+snpArcher is a reproducible workflow optimized for nonmodel organisms and comparisons across datasets, built on the [Snakemake](https://snakemake.readthedocs.io/en/stable/index.html#) workflow management system. It provides a streamlined approach to dataset acquisition, variant calling, quality control, and downstream analysis.
-### Sample metadata sheet
+### Usage
+For usage instructions and complete documentation, please visit our [docs](https://snparcher.readthedocs.io/en/latest/).
-The workflow requires a comma seperated metadata sheet to run. The fields of the sheet are outlined below:
-| Field | Description |
-| ---- | -------------|
-| BioSample | The name of the sample. |
-| LibraryName | LibraryID for sample, **must be unique.** |
-| Run | The SRR for the sample, if applicable. If not, must be some **unique** value. |
-| Organism | The name of the organism. |
-| refGenome | Reference genome accession, if applicable. *See note* |
-| refPath | Path to local reference genome, if applicable. *See note* |
-| BioProject | If applicable. Otherwise any value is acceptable. |
-| fq1 | Optional. Path to read 1 for sample |
-| fq2 | Optional. Path to read 2 for sample |
+### Datasets generated by snpArcher
+A number of resequencing datasets have been run with snpArcher generating consistent variant calls, available via [Globus](https://www.globus.org/) in the [Comparative Population Genomics Data collection](https://app.globus.org/file-manager?origin_id=a6580c44-09fd-11ee-be16-195c41bc0be4&origin_path=%2F). Details of data processing are described [in our manuscript](https://www.biorxiv.org/content/10.1101/2023.06.22.546168v1). If you use any of these datasets in your projects, please cite both the [snpArcher paper](https://www.biorxiv.org/content/10.1101/2023.06.22.546168v1) and the original data producers.
-*Note:* refGenome is always required. refPath is optional, but when specified, a name for the assembly (in refGenome) must also be included.
-
-It is important to note that samples are proccessed together based on their `refGenome` metadata, so **all BioSamples that share a reference genome will ultimately end up in the same final vcf file.** If you are mapping multiple populations / species to a single reference genome, and want separate VCF files for each population / species, you will need to split your final vcf after the pipeline completes, or run multiiple indpendent sample sheets in different results directories.
-
-If your reads (and, optionally, your local reference genome) are stored in somewhere seperate of the workflow (e.g.: a scratch disk) then you can specify the path to your reads using the `fq1` and `fq2` fields, and the location of your reference genome fasta (*note: must be uncompressed*) in the `refPath` field.
-
-A python script `workflow/write_samples.py` is included to help write the sample sheet for you. This script has three required arguments:
-|Argument| Description|
-| ------ | ---------- |
-| `-s / --sample_list` | Path to a sample list. One sample per line |
-| `-f / --fastq_dir` | Path to directory containing ALL fastq files. It is assumed that each fastq file will contain the sample name uniquely. |
-| `-o / --org` | The organism name. |
-
-Additionally, one of the following options must also be specified:
-|Argument| Description|
-| ------ | ---------- |
-| `-r / --ref` | Path to reference fasta. |
-| `-a / --acc` | NCBI accession of reference. |
-
-### Workflow configuration
-
-The other file that needs to be updated can be found under `config/config.yml`:
-
-```
-##############################
-# Variables you need to change
-##############################
-
-samples: "config/samples.csv" # path to the sample metadata CSV
-resource_config: "config/resources.yaml" # path to resources yaml config
-final_prefix: "" # prefix for final output files
-intervals: True #Set to True if you want to perform variant calling using interval approach.
-sentieon: False #set to True if you want to use sentieon, False if you want GATK
-sentieon_lic: "" #set to path of sentieon license
-remote_reads: False # Set True if reads are in a Google Bucket seperate from --default-remote-prefix.
-remote_reads_prefix: "" # set to google bucket prefix where reads live
-bigtmp: "" #Set to a path with lots of free space to use for commands that require large amounts of temp space; defaults to system tmpdir if empty
-cov_filter: True #set to True if you want to include coverage thresholds in the callable sites bed file (default uses mappability only)
-
-```
-You are required to set a `final_prefix`: your final vcf file will be called {final_prefix}_final.vcf. To run out of the box, we recommend leaving `intervals` True and `sention` and `remote_reads` (for cloud operations) False. You will need to change `samples` to point to your sample sheet. You may need to set `bigtmp` to something depending on your system configuration, as some steps can requires 100+ Gb of temp space for large datasets; `tmp/` is a decent option, which will create a temp directory in the directory you launch snpArcher from.
-
-The default `cov_filter` is very loose, removing only regions of the genome with 0 coverage and excessively high coverage (>10000x). You can modify these parameters in the config, as needed. The default options are also optimized for low coverage (<10x) data. If you have >15x coverage, you probably want to use the high coverage parameters instead, which can also be changed in the config.
-
-Increasing the `num_gvcf_intervals` will scale the pipeline wider, so there will be more, shorter jobs. Decreasing `num_gvcf_intervals` will create fewer, longer jobs. The optimal setup will depend on your HPC system.
-
-### Run the workflow
-
-Execute the Snakemake workflow by running the command:
-`snakemake --use-conda --cores <# of cores to use>`
-
-If you are running on an HPC system, you'll need to include a profile for, e.g. `slurm` or `sge`:
-`snakemake --profile profiles/slurm`
-
-We include a slurm profile with snpArcher, but you will need to modify the queue names, at a minimum, for your institution. We welcome contributions of other profiles. The `run_pipeline.sh` script is an example of how you might launch a snakemake job using slurm, and assuming you have snakemake installed in a conda environment called snakemake, although you will need to modify the partition names and possibly the resource requests depending on your local environment (and we would recommend changing the out and err files to something more informative, e.g. {final_prefix}-sm-%j.out and {final_prefix}-sm-%j.err).
-
-## Run example data
-
-To run the test data, activate your snakemake conda environment and execute the following command:
-`snakemake -d .test/ecoli --use-conda --cores <# of cores to use>`
-
-## Output
-
-Output can be found in the `results` folder and includes everal key files:
-
-- `results/{ASSEMBLY_NAME}`
-
-The main output of the pipeline is a single VCF with genotype calls for every individual:
-
-- `results/{ASSEMBLY_NAME}/{final_prefix}_final.vcf.gz`
-
-By default, this file contains all SNPs and Indels identified and has the basic GATK filters applied. No filtering has been done on the VCF, so it will include all individuals from the sample sheet and all variants identified. The filters are applied as annotations within the VCF file.
-
-A very simple example for removing filtered sites and only retaining biallelic SNPs is, e.g.:
-
-```
-bcftools view -v snps -m2 -M2 -f .,PASS -e 'AF==1 | AF==0 | ALT="*" | TYPE~"indel" | ref="N"' {input.vcf} -O z -o {output.filtered}
-```
+### Citing snpArcher
+- Cade D Mirchandani, Allison J Shultz, Gregg W C Thomas, Sara J Smith, Mara Baylis, Brian Arnold, Russ Corbett-Detig, Erik Enbody, Timothy B Sackton, A fast, reproducible, high-throughput variant calling workflow for population genomics, Molecular Biology and Evolution, 2023;, msad270, https://doi.org/10.1093/molbev/msad270
+- Also, make sure to cite the tools you used within snpArcher.
diff --git a/config/config.yaml b/config/config.yaml
index 12a92620..b5dd2068 100644
--- a/config/config.yaml
+++ b/config/config.yaml
@@ -3,12 +3,11 @@
##############################
samples: "config/samples.csv" # path to the sample metadata CSV
-resource_config: "config/resources.yaml" # path to resources yaml config
final_prefix: "" # prefix for final output files
intervals: True #Set to True if you want to perform variant calling using interval approach.
sentieon: False #set to True if you want to use sentieon, False if you want GATK
sentieon_lic: "" #set to path of sentieon license
-remote_reads: False # Set True if reads are in a Google Bucket seperate from --default-remote-prefix.
+remote_reads: False # Set True if reads are in a location seperate from --default-remote-prefix.
remote_reads_prefix: "" # set to google bucket prefix where reads live
bigtmp: "" #Set to a path with lots of free space to use for commands that require large amounts of temp space; defaults to system tmpdir if empty
cov_filter: True #set to True if you want to include coverage thresholds in the callable sites bed file (default uses mappability only)
@@ -18,11 +17,15 @@ trackhub_email: ""
# Variables you *might* need to change
##############################
+# Set reference genome here if you would like to you use the same reference genome for all samples in sample sheet. See docs for more info.
+# refGenome:
+# refPath:
+
# Interval approach options, only applicable if intervals is True
minNmer: 500 # the minimum Nmer used to split up the genome; e.g. a value of 200 means only Nmers 200 or greater are used to define the boundaries of intervals. The minimum is 50.
num_gvcf_intervals: 50 # The maximum number of intervals to create for GVCF generation. Note: the actual number of intervals may be less than the specified value if the reference genome has very few gaps.
db_scatter_factor: 0.15 # Scatter factor for calculating number of intervals to create for genomics db generation. (scatter_factor * num_samples * num_gvcf_intervals) gives us number of db intervals to create. Reccomend <1
-
+ploidy: 2 # Ploidy for HaplotypeCaller and Sentieon Haplotyper
## Coverage options ##
## default pipeline is optimized for low coverage data - if using high coverage data (> 10x), uncomment high coverage options and comment out low coverage options
diff --git a/config/resources.yaml b/config/resources.yaml
deleted file mode 100644
index 05eeef12..00000000
--- a/config/resources.yaml
+++ /dev/null
@@ -1,121 +0,0 @@
-###
-# fastq2bam rules
-##
-
-# fastq download
-get_fastq_pe:
- threads: 8
- mem: 4000
-
-# index reference
-index_ref:
- mem: 10000
-
-# fastp program
-fastp:
- threads: 8
- mem: 4000
-
-# bwa mapping
-bwa_map:
- threads: 31
- mem: 15000
-# sort bam with picard's SortSam tool
-sort_bam:
- threads: 4
- mem_per_thread: 1000
-
-#merge bams
-merge_bams:
- mem: 9000
- threads: 2
-
-# mark duplicates with picard's MarkDuplicates
-dedup:
- threads: 31
- mem: 9000
-# calculate BAM summaries with samtools and picard
-bam_sumstats:
- mem: 9000
-merge_bams:
- mem: 9000
-bedtools:
- mem: 4000
-
-# Sentieon tools
-sentieon_map:
- machine_type: "n2d-standard-32"
- threads: 31
- mem: 15000
-
-sentieon_dedup:
- machine_type: "n2d-standard-32"
- threads: 31
- mem: 15000
-
-sentieon_haplotyper:
- machine_type: "n2d-standard-32"
- threads: 31
- mem: 15000
-
-sentieon_combine_gvcf:
- machine_type: "n2d-standard-32"
- threads: 31
- mem: 15000
- disk_mb: 2000000
-
-###
-# Intervals workflow
-###
-
-# preprocess genome, create intervals
-# picard's create CreateSequenceDictionary, samtools faidx
-process_ref:
- mem: 15000
-# custom python algo to create intervals
-create_intervals:
- mem: 5000
-
-## Callable sites workflow
-
-# genmap map
-genmap:
- threads: 10
- mem: 20000
-genmap_sort:
- mem: 4000
-compute_d4:
- mem: 4000
- threads: 4
-merge_d4:
- mem: 10000
-callable_bed:
- mem: 10000
-
-###
-# bam2vcf workflows
-###
-
-# gatk HaplotypeCaller
-bam2gvcf:
- mem: 9000
-# gatk GenomicsDBImport
-gvcf2DB:
- mem: 32000
-# gatk GenotypeGVCFs
-DB2vcf:
- mem: 32000
-# gatk filterVcfs
-filterVcfs:
- mem: 4000
-# gatk GatherVcfs
-gatherVcfs:
- mem: 9000
-# vcftools program
-vcftools:
- mem: 8000
-# plink
-plink:
- threads: 5
-admixture:
- mem: 4000
diff --git a/config/samples.csv b/config/samples.csv
deleted file mode 100644
index 19a8cd84..00000000
--- a/config/samples.csv
+++ /dev/null
@@ -1,2 +0,0 @@
-BioSample,LibraryName,refGenome,Run,Organism,BioProject
-testENA,EK7.12,GCA_000008865.2,SRR10058855,Escherichia coli,PRJNA563564
diff --git a/config/test.csv b/config/test.csv
deleted file mode 100644
index b35e9d55..00000000
--- a/config/test.csv
+++ /dev/null
@@ -1,3 +0,0 @@
-BioSample,LibraryName,refGenome,Run,Organism,BioProject,
-sample,sample1,genome1,1,test,x
-
diff --git a/config/test_coords.csv b/config/test_coords.csv
deleted file mode 100644
index 62489259..00000000
--- a/config/test_coords.csv
+++ /dev/null
@@ -1,247 +0,0 @@
-BioSample,LibraryName,refGenome,Organism,fq1,fq2,lat,long,Run,BioProject
-test_A01,test_A01,genome1,test,test_A01_R1.fastq.gz,test_A01_R2.fastq.gz,39.4433764,-122.8245446,test_A01,x
-test_A01,test_A01_S1_L001,genome1,test,test_A01_S1_L001_R1_001.fastq.gz,test_A01_S1_L001_R2_001.fastq.gz,39.4433764,-122.8245446,test_A01_S1_L001,x
-test_B01,test_B01,genome1,test,test_B01_R1.fastq.gz,test_B01_R2.fastq.gz,39.4433764,-122.8245446,test_B01,x
-test_B01,test_B01_S13_L001,genome1,test,test_B01_S13_L001_R1_001.fastq.gz,test_B01_S13_L001_R2_001.fastq.gz,39.4433764,-122.8245446,test_B01_S13_L001,x
-test_C01,test_C01_S25_L001,genome1,test,test_C01_S25_L001_R1_001.fastq.gz,test_C01_S25_L001_R2_001.fastq.gz,39.4433764,-122.8245446,test_C01_S25_L001,x
-test_C01,test_C01,genome1,test,test_C01_R1.fastq.gz,test_C01_R2.fastq.gz,39.4433764,-122.8245446,test_C01,x
-test_D01,test_D01,genome1,test,test_D01_R1.fastq.gz,test_D01_R2.fastq.gz,39.3759636,-122.866367,test_D01,x
-test_D01,test_D01_S37_L001,genome1,test,test_D01_S37_L001_R1_001.fastq.gz,test_D01_S37_L001_R2_001.fastq.gz,39.3759636,-122.866367,test_D01_S37_L001,x
-test_E01,test_E01,genome1,test,test_E01_R1.fastq.gz,test_E01_R2.fastq.gz,39.40491,-122.62053,test_E01,x
-test_E01,test_E01_S49_L001,genome1,test,test_E01_S49_L001_R1_001.fastq.gz,test_E01_S49_L001_R2_001.fastq.gz,39.40491,-122.62053,test_E01_S49_L001,x
-test_F01,test_F01_S61_L001,genome1,test,test_F01_S61_L001_R1_001.fastq.gz,test_F01_S61_L001_R2_001.fastq.gz,39.40491,-122.62053,test_F01_S61_L001,x
-test_F01,test_F01,genome1,test,test_F01_R1.fastq.gz,test_F01_R2.fastq.gz,39.40491,-122.62053,test_F01,x
-test_G01,test_G01,genome1,test,test_G01_R1.fastq.gz,test_G01_R2.fastq.gz,39.40491,-122.62053,test_G01,x
-test_G01,test_G01_S73_L001,genome1,test,test_G01_S73_L001_R1_001.fastq.gz,test_G01_S73_L001_R2_001.fastq.gz,39.40491,-122.62053,test_G01_S73_L001,x
-test_H01,test_H01,genome1,test,test_H01_R1.fastq.gz,test_H01_R2.fastq.gz,35.0604116,-119.6941782,test_H01,x
-test_H01,test_H01_S85_L001,genome1,test,test_H01_S85_L001_R1_001.fastq.gz,test_H01_S85_L001_R2_001.fastq.gz,35.0604116,-119.6941782,test_H01_S85_L001,x
-test_A02,test_A02,genome1,test,test_A02_R1.fastq.gz,test_A02_R2.fastq.gz,35.0604116,-119.6941782,test_A02,x
-test_A02,test_A02_S2_L001,genome1,test,test_A02_S2_L001_R1_001.fastq.gz,test_A02_S2_L001_R2_001.fastq.gz,35.0604116,-119.6941782,test_A02_S2_L001,x
-test_B02,test_B02,genome1,test,test_B02_R1.fastq.gz,test_B02_R2.fastq.gz,39.3759636,-122.866367,test_B02,x
-test_B02,test_B02_S14_L001,genome1,test,test_B02_S14_L001_R1_001.fastq.gz,test_B02_S14_L001_R2_001.fastq.gz,39.3759636,-122.866367,test_B02_S14_L001,x
-test_C02,test_C02_S26_L001,genome1,test,test_C02_S26_L001_R1_001.fastq.gz,test_C02_S26_L001_R2_001.fastq.gz,39.3759636,-122.866367,test_C02_S26_L001,x
-test_C02,test_C02,genome1,test,test_C02_R1.fastq.gz,test_C02_R2.fastq.gz,39.3759636,-122.866367,test_C02,x
-test_D02,test_D02,genome1,test,test_D02_R1.fastq.gz,test_D02_R2.fastq.gz,35.0604116,-119.6941782,test_D02,x
-test_D02,test_D02_S38_L001,genome1,test,test_D02_S38_L001_R1_001.fastq.gz,test_D02_S38_L001_R2_001.fastq.gz,35.0604116,-119.6941782,test_D02_S38_L001,x
-test_E02,test_E02,genome1,test,test_E02_R1.fastq.gz,test_E02_R2.fastq.gz,36.7042146,-120.8221395,test_E02,x
-test_E02,test_E02_S50_L001,genome1,test,test_E02_S50_L001_R1_001.fastq.gz,test_E02_S50_L001_R2_001.fastq.gz,36.7042146,-120.8221395,test_E02_S50_L001,x
-test_F02,test_F02_S62_L001,genome1,test,test_F02_S62_L001_R1_001.fastq.gz,test_F02_S62_L001_R2_001.fastq.gz,36.7042146,-120.8221395,test_F02_S62_L001,x
-test_F02,test_F02,genome1,test,test_F02_R1.fastq.gz,test_F02_R2.fastq.gz,36.7042146,-120.8221395,test_F02,x
-test_G02,test_G02,genome1,test,test_G02_R1.fastq.gz,test_G02_R2.fastq.gz,35.36144,-120.31198,test_G02,x
-test_G02,test_G02_S74_L001,genome1,test,test_G02_S74_L001_R1_001.fastq.gz,test_G02_S74_L001_R2_001.fastq.gz,35.36144,-120.31198,test_G02_S74_L001,x
-test_H02,test_H02_S86_L001,genome1,test,test_H02_S86_L001_R1_001.fastq.gz,test_H02_S86_L001_R2_001.fastq.gz,34.61753,-118.58407,test_H02_S86_L001,x
-test_H02,test_H02,genome1,test,test_H02_R1.fastq.gz,test_H02_R2.fastq.gz,34.61753,-118.58407,test_H02,x
-test_A03,test_A03_S3_L001,genome1,test,test_A03_S3_L001_R1_001.fastq.gz,test_A03_S3_L001_R2_001.fastq.gz,34.61753,-118.58407,test_A03_S3_L001,x
-test_A03,test_A03,genome1,test,test_A03_R1.fastq.gz,test_A03_R2.fastq.gz,34.61753,-118.58407,test_A03,x
-test_B03,test_B03_S15_L001,genome1,test,test_B03_S15_L001_R1_001.fastq.gz,test_B03_S15_L001_R2_001.fastq.gz,34.61753,-118.58407,test_B03_S15_L001,x
-test_B03,test_B03,genome1,test,test_B03_R1.fastq.gz,test_B03_R2.fastq.gz,34.61753,-118.58407,test_B03,x
-test_C03,test_C03,genome1,test,test_C03_R1.fastq.gz,test_C03_R2.fastq.gz,32.98273,-116.79362,test_C03,x
-test_C03,test_C03_S27_L001,genome1,test,test_C03_S27_L001_R1_001.fastq.gz,test_C03_S27_L001_R2_001.fastq.gz,32.98273,-116.79362,test_C03_S27_L001,x
-test_D03,test_D03_S39_L001,genome1,test,test_D03_S39_L001_R1_001.fastq.gz,test_D03_S39_L001_R2_001.fastq.gz,32.98273,-116.79362,test_D03_S39_L001,x
-test_D03,test_D03,genome1,test,test_D03_R1.fastq.gz,test_D03_R2.fastq.gz,32.98273,-116.79362,test_D03,x
-test_E03,test_E03_S51_L001,genome1,test,test_E03_S51_L001_R1_001.fastq.gz,test_E03_S51_L001_R2_001.fastq.gz,32.98273,-116.79362,test_E03_S51_L001,x
-test_E03,test_E03,genome1,test,test_E03_R1.fastq.gz,test_E03_R2.fastq.gz,32.98273,-116.79362,test_E03,x
-test_F03,test_F03,genome1,test,test_F03_R1.fastq.gz,test_F03_R2.fastq.gz,32.7120598,-116.2312235,test_F03,x
-test_F03,test_F03_S63_L001,genome1,test,test_F03_S63_L001_R1_001.fastq.gz,test_F03_S63_L001_R2_001.fastq.gz,32.7120598,-116.2312235,test_F03_S63_L001,x
-test_G03,test_G03_S75_L001,genome1,test,test_G03_S75_L001_R1_001.fastq.gz,test_G03_S75_L001_R2_001.fastq.gz,32.96822,-116.82805,test_G03_S75_L001,x
-test_G03,test_G03,genome1,test,test_G03_R1.fastq.gz,test_G03_R2.fastq.gz,32.96822,-116.82805,test_G03,x
-test_H03,test_H03,genome1,test,test_H03_R1.fastq.gz,test_H03_R2.fastq.gz,32.96822,-116.82805,test_H03,x
-test_H03,test_H03_S87_L001,genome1,test,test_H03_S87_L001_R1_001.fastq.gz,test_H03_S87_L001_R2_001.fastq.gz,32.96822,-116.82805,test_H03_S87_L001,x
-test_A04,test_A04,genome1,test,test_A04_R1.fastq.gz,test_A04_R2.fastq.gz,40.3160017,-122.8557539,test_A04,x
-test_A04,test_A04_S4_L001,genome1,test,test_A04_S4_L001_R1_001.fastq.gz,test_A04_S4_L001_R2_001.fastq.gz,40.3160017,-122.8557539,test_A04_S4_L001,x
-test_B04,test_B04_S16_L001,genome1,test,test_B04_S16_L001_R1_001.fastq.gz,test_B04_S16_L001_R2_001.fastq.gz,40.3160017,-122.8557539,test_B04_S16_L001,x
-test_B04,test_B04,genome1,test,test_B04_R1.fastq.gz,test_B04_R2.fastq.gz,40.3160017,-122.8557539,test_B04,x
-test_C04,test_C04,genome1,test,test_C04_R1.fastq.gz,test_C04_R2.fastq.gz,40.3160017,-122.8557539,test_C04,x
-test_C04,test_C04_S28_L001,genome1,test,test_C04_S28_L001_R1_001.fastq.gz,test_C04_S28_L001_R2_001.fastq.gz,40.3160017,-122.8557539,test_C04_S28_L001,x
-test_D04,test_D04_S40_L001,genome1,test,test_D04_S40_L001_R1_001.fastq.gz,test_D04_S40_L001_R2_001.fastq.gz,36.7042146,-120.8221395,test_D04_S40_L001,x
-test_D04,test_D04,genome1,test,test_D04_R1.fastq.gz,test_D04_R2.fastq.gz,36.7042146,-120.8221395,test_D04,x
-test_E04,test_E04_S52_L001,genome1,test,test_E04_S52_L001_R1_001.fastq.gz,test_E04_S52_L001_R2_001.fastq.gz,35.31353,-118.08886,test_E04_S52_L001,x
-test_E04,test_E04,genome1,test,test_E04_R1.fastq.gz,test_E04_R2.fastq.gz,35.31353,-118.08886,test_E04,x
-test_F04,test_F04,genome1,test,test_F04_R1.fastq.gz,test_F04_R2.fastq.gz,35.31353,-118.08886,test_F04,x
-test_F04,test_F04_S64_L001,genome1,test,test_F04_S64_L001_R1_001.fastq.gz,test_F04_S64_L001_R2_001.fastq.gz,35.31353,-118.08886,test_F04_S64_L001,x
-test_G04,test_G04_S76_L001,genome1,test,test_G04_S76_L001_R1_001.fastq.gz,test_G04_S76_L001_R2_001.fastq.gz,35.31353,-118.08886,test_G04_S76_L001,x
-test_G04,test_G04,genome1,test,test_G04_R1.fastq.gz,test_G04_R2.fastq.gz,35.31353,-118.08886,test_G04,x
-test_H04,test_H04_S88_L001,genome1,test,test_H04_S88_L001_R1_001.fastq.gz,test_H04_S88_L001_R2_001.fastq.gz,35.3572,-120.30679,test_H04_S88_L001,x
-test_H04,test_H04,genome1,test,test_H04_R1.fastq.gz,test_H04_R2.fastq.gz,35.3572,-120.30679,test_H04,x
-test_A05,test_A05_S5_L001,genome1,test,test_A05_S5_L001_R1_001.fastq.gz,test_A05_S5_L001_R2_001.fastq.gz,35.3572,-120.30679,test_A05_S5_L001,x
-test_A05,test_A05,genome1,test,test_A05_R1.fastq.gz,test_A05_R2.fastq.gz,35.3572,-120.30679,test_A05,x
-test_B05,test_B05_S17_L001,genome1,test,test_B05_S17_L001_R1_001.fastq.gz,test_B05_S17_L001_R2_001.fastq.gz,35.3572,-120.30679,test_B05_S17_L001,x
-test_B05,test_B05,genome1,test,test_B05_R1.fastq.gz,test_B05_R2.fastq.gz,35.3572,-120.30679,test_B05,x
-test_C05,test_C05,genome1,test,test_C05_R1.fastq.gz,test_C05_R2.fastq.gz,35.3572,-120.30679,test_C05,x
-test_C05,test_C05_S29_L001,genome1,test,test_C05_S29_L001_R1_001.fastq.gz,test_C05_S29_L001_R2_001.fastq.gz,35.3572,-120.30679,test_C05_S29_L001,x
-test_D05,test_D05,genome1,test,test_D05_R1.fastq.gz,test_D05_R2.fastq.gz,35.3572,-120.30679,test_D05,x
-test_D05,test_D05_S41_L001,genome1,test,test_D05_S41_L001_R1_001.fastq.gz,test_D05_S41_L001_R2_001.fastq.gz,35.3572,-120.30679,test_D05_S41_L001,x
-test_E05,test_E05_S53_L001,genome1,test,test_E05_S53_L001_R1_001.fastq.gz,test_E05_S53_L001_R2_001.fastq.gz,35.2680047,-120.4482756,test_E05_S53_L001,x
-test_E05,test_E05,genome1,test,test_E05_R1.fastq.gz,test_E05_R2.fastq.gz,35.2680047,-120.4482756,test_E05,x
-test_F05,test_F05,genome1,test,test_F05_R1.fastq.gz,test_F05_R2.fastq.gz,35.2680047,-120.4482756,test_F05,x
-test_F05,test_F05_S65_L001,genome1,test,test_F05_S65_L001_R1_001.fastq.gz,test_F05_S65_L001_R2_001.fastq.gz,35.2680047,-120.4482756,test_F05_S65_L001,x
-test_G05,test_G05,genome1,test,test_G05_R1.fastq.gz,test_G05_R2.fastq.gz,35.2680047,-120.4482756,test_G05,x
-test_G05,test_G05_S77_L001,genome1,test,test_G05_S77_L001_R1_001.fastq.gz,test_G05_S77_L001_R2_001.fastq.gz,35.2680047,-120.4482756,test_G05_S77_L001,x
-test_H05,test_H05,genome1,test,test_H05_R1.fastq.gz,test_H05_R2.fastq.gz,35.2680047,-120.4482756,test_H05,x
-test_H05,test_H05_S89_L001,genome1,test,test_H05_S89_L001_R1_001.fastq.gz,test_H05_S89_L001_R2_001.fastq.gz,35.2680047,-120.4482756,test_H05_S89_L001,x
-test_A06,test_A06,genome1,test,test_A06_R1.fastq.gz,test_A06_R2.fastq.gz,36.04454,-117.91501,test_A06,x
-test_A06,test_A06_S6_L001,genome1,test,test_A06_S6_L001_R1_001.fastq.gz,test_A06_S6_L001_R2_001.fastq.gz,36.04454,-117.91501,test_A06_S6_L001,x
-test_B06,test_B06_S18_L001,genome1,test,test_B06_S18_L001_R1_001.fastq.gz,test_B06_S18_L001_R2_001.fastq.gz,36.04454,-117.91501,test_B06_S18_L001,x
-test_B06,test_B06,genome1,test,test_B06_R1.fastq.gz,test_B06_R2.fastq.gz,36.04454,-117.91501,test_B06,x
-test_C06,test_C06,genome1,test,test_C06_R1.fastq.gz,test_C06_R2.fastq.gz,36.04454,-117.91501,test_C06,x
-test_C06,test_C06_S30_L001,genome1,test,test_C06_S30_L001_R1_001.fastq.gz,test_C06_S30_L001_R2_001.fastq.gz,36.04454,-117.91501,test_C06_S30_L001,x
-test_D06,test_D06_S42_L001,genome1,test,test_D06_S42_L001_R1_001.fastq.gz,test_D06_S42_L001_R2_001.fastq.gz,37.69792,-120.26773,test_D06_S42_L001,x
-test_D06,test_D06,genome1,test,test_D06_R1.fastq.gz,test_D06_R2.fastq.gz,37.69792,-120.26773,test_D06,x
-test_E06,test_E06,genome1,test,test_E06_R1.fastq.gz,test_E06_R2.fastq.gz,37.69792,-120.26773,test_E06,x
-test_E06,test_E06_S54_L001,genome1,test,test_E06_S54_L001_R1_001.fastq.gz,test_E06_S54_L001_R2_001.fastq.gz,37.69792,-120.26773,test_E06_S54_L001,x
-test_F06,test_F06_S66_L001,genome1,test,test_F06_S66_L001_R1_001.fastq.gz,test_F06_S66_L001_R2_001.fastq.gz,37.69792,-120.26773,test_F06_S66_L001,x
-test_F06,test_F06,genome1,test,test_F06_R1.fastq.gz,test_F06_R2.fastq.gz,37.69792,-120.26773,test_F06,x
-test_G06,test_G06_S78_L001,genome1,test,test_G06_S78_L001_R1_001.fastq.gz,test_G06_S78_L001_R2_001.fastq.gz,37.69792,-120.26773,test_G06_S78_L001,x
-test_G06,test_G06,genome1,test,test_G06_R1.fastq.gz,test_G06_R2.fastq.gz,37.69792,-120.26773,test_G06,x
-test_H06,test_H06,genome1,test,test_H06_R1.fastq.gz,test_H06_R2.fastq.gz,36.16023,-121.48945,test_H06,x
-test_H06,test_H06_S90_L001,genome1,test,test_H06_S90_L001_R1_001.fastq.gz,test_H06_S90_L001_R2_001.fastq.gz,36.16023,-121.48945,test_H06_S90_L001,x
-test_A07,test_A07_S7_L001,genome1,test,test_A07_S7_L001_R1_001.fastq.gz,test_A07_S7_L001_R2_001.fastq.gz,36.16023,-121.48945,test_A07_S7_L001,x
-test_A07,test_A07,genome1,test,test_A07_R1.fastq.gz,test_A07_R2.fastq.gz,36.16023,-121.48945,test_A07,x
-test_B07,test_B07,genome1,test,test_B07_R1.fastq.gz,test_B07_R2.fastq.gz,36.16023,-121.48945,test_B07,x
-test_B07,test_B07_S19_L001,genome1,test,test_B07_S19_L001_R1_001.fastq.gz,test_B07_S19_L001_R2_001.fastq.gz,36.16023,-121.48945,test_B07_S19_L001,x
-test_C07,test_C07_S31_L001,genome1,test,test_C07_S31_L001_R1_001.fastq.gz,test_C07_S31_L001_R2_001.fastq.gz,37.4343008,-118.3681869,test_C07_S31_L001,x
-test_C07,test_C07,genome1,test,test_C07_R1.fastq.gz,test_C07_R2.fastq.gz,37.4343008,-118.3681869,test_C07,x
-test_D07,test_D07_S43_L001,genome1,test,test_D07_S43_L001_R1_001.fastq.gz,test_D07_S43_L001_R2_001.fastq.gz,37.4343008,-118.3681869,test_D07_S43_L001,x
-test_D07,test_D07,genome1,test,test_D07_R1.fastq.gz,test_D07_R2.fastq.gz,37.4343008,-118.3681869,test_D07,x
-test_E07,test_E07,genome1,test,test_E07_R1.fastq.gz,test_E07_R2.fastq.gz,37.4343008,-118.3681869,test_E07,x
-test_E07,test_E07_S55_L001,genome1,test,test_E07_S55_L001_R1_001.fastq.gz,test_E07_S55_L001_R2_001.fastq.gz,37.4343008,-118.3681869,test_E07_S55_L001,x
-test_F07,test_F07,genome1,test,test_F07_R1.fastq.gz,test_F07_R2.fastq.gz,37.4306649,-118.4182199,test_F07,x
-test_F07,test_F07_S67_L001,genome1,test,test_F07_S67_L001_R1_001.fastq.gz,test_F07_S67_L001_R2_001.fastq.gz,37.4306649,-118.4182199,test_F07_S67_L001,x
-test_G07,test_G07,genome1,test,test_G07_R1.fastq.gz,test_G07_R2.fastq.gz,37.4012819,-118.4272425,test_G07,x
-test_G07,test_G07_S79_L001,genome1,test,test_G07_S79_L001_R1_001.fastq.gz,test_G07_S79_L001_R2_001.fastq.gz,37.4012819,-118.4272425,test_G07_S79_L001,x
-test_H07,test_H07,genome1,test,test_H07_R1.fastq.gz,test_H07_R2.fastq.gz,37.3521172,-118.3503229,test_H07,x
-test_H07,test_H07_S91_L001,genome1,test,test_H07_S91_L001_R1_001.fastq.gz,test_H07_S91_L001_R2_001.fastq.gz,37.3521172,-118.3503229,test_H07_S91_L001,x
-test_A08,test_A08_S8_L001,genome1,test,test_A08_S8_L001_R1_001.fastq.gz,test_A08_S8_L001_R2_001.fastq.gz,37.3521172,-118.3503229,test_A08_S8_L001,x
-test_A08,test_A08,genome1,test,test_A08_R1.fastq.gz,test_A08_R2.fastq.gz,37.3521172,-118.3503229,test_A08,x
-test_B08,test_B08,genome1,test,test_B08_R1.fastq.gz,test_B08_R2.fastq.gz,37.3365136,-118.3493757,test_B08,x
-test_B08,test_B08_S20_L001,genome1,test,test_B08_S20_L001_R1_001.fastq.gz,test_B08_S20_L001_R2_001.fastq.gz,37.3365136,-118.3493757,test_B08_S20_L001,x
-test_C08,test_C08,genome1,test,test_C08_R1.fastq.gz,test_C08_R2.fastq.gz,37.3294551,-118.4641188,test_C08,x
-test_C08,test_C08_S32_L001,genome1,test,test_C08_S32_L001_R1_001.fastq.gz,test_C08_S32_L001_R2_001.fastq.gz,37.3294551,-118.4641188,test_C08_S32_L001,x
-test_D08,test_D08,genome1,test,test_D08_R1.fastq.gz,test_D08_R2.fastq.gz,37.3222784,-118.4276573,test_D08,x
-test_D08,test_D08_S44_L001,genome1,test,test_D08_S44_L001_R1_001.fastq.gz,test_D08_S44_L001_R2_001.fastq.gz,37.3222784,-118.4276573,test_D08_S44_L001,x
-test_E08,test_E08,genome1,test,test_E08_R1.fastq.gz,test_E08_R2.fastq.gz,37.3222784,-118.4276573,test_E08,x
-test_E08,test_E08_S56_L001,genome1,test,test_E08_S56_L001_R1_001.fastq.gz,test_E08_S56_L001_R2_001.fastq.gz,37.3222784,-118.4276573,test_E08_S56_L001,x
-test_F08,test_F08,genome1,test,test_F08_R1.fastq.gz,test_F08_R2.fastq.gz,37.2928647,-118.3130339,test_F08,x
-test_F08,test_F08_S68_L001,genome1,test,test_F08_S68_L001_R1_001.fastq.gz,test_F08_S68_L001_R2_001.fastq.gz,37.2928647,-118.3130339,test_F08_S68_L001,x
-test_G08,test_G08_S80_L001,genome1,test,test_G08_S80_L001_R1_001.fastq.gz,test_G08_S80_L001_R2_001.fastq.gz,37.2928647,-118.3130339,test_G08_S80_L001,x
-test_G08,test_G08,genome1,test,test_G08_R1.fastq.gz,test_G08_R2.fastq.gz,37.2928647,-118.3130339,test_G08,x
-test_H08,test_H08,genome1,test,test_H08_R1.fastq.gz,test_H08_R2.fastq.gz,37.2928647,-118.3130339,test_H08,x
-test_H08,test_H08_S92_L001,genome1,test,test_H08_S92_L001_R1_001.fastq.gz,test_H08_S92_L001_R2_001.fastq.gz,37.2928647,-118.3130339,test_H08_S92_L001,x
-test_A09,test_A09,genome1,test,test_A09_R1.fastq.gz,test_A09_R2.fastq.gz,37.2928647,-118.3130339,test_A09,x
-test_A09,test_A09_S9_L001,genome1,test,test_A09_S9_L001_R1_001.fastq.gz,test_A09_S9_L001_R2_001.fastq.gz,37.2928647,-118.3130339,test_A09_S9_L001,x
-test_B09,test_B09,genome1,test,test_B09_R1.fastq.gz,test_B09_R2.fastq.gz,37.13932,-118.2707,test_B09,x
-test_B09,test_B09_S21_L001,genome1,test,test_B09_S21_L001_R1_001.fastq.gz,test_B09_S21_L001_R2_001.fastq.gz,37.13932,-118.2707,test_B09_S21_L001,x
-test_C09,test_C09,genome1,test,test_C09_R1.fastq.gz,test_C09_R2.fastq.gz,37.13932,-118.2707,test_C09,x
-test_C09,test_C09_S33_L001,genome1,test,test_C09_S33_L001_R1_001.fastq.gz,test_C09_S33_L001_R2_001.fastq.gz,37.13932,-118.2707,test_C09_S33_L001,x
-test_D09,test_D09_S45_L001,genome1,test,test_D09_S45_L001_R1_001.fastq.gz,test_D09_S45_L001_R2_001.fastq.gz,36.80243,-118.15371,test_D09_S45_L001,x
-test_D09,test_D09,genome1,test,test_D09_R1.fastq.gz,test_D09_R2.fastq.gz,36.80243,-118.15371,test_D09,x
-test_E09,test_E09,genome1,test,test_E09_R1.fastq.gz,test_E09_R2.fastq.gz,36.80243,-118.15371,test_E09,x
-test_E09,test_E09_S57_L001,genome1,test,test_E09_S57_L001_R1_001.fastq.gz,test_E09_S57_L001_R2_001.fastq.gz,36.80243,-118.15371,test_E09_S57_L001,x
-test_F09,test_F09,genome1,test,test_F09_R1.fastq.gz,test_F09_R2.fastq.gz,36.80243,-118.15371,test_F09,x
-test_F09,test_F09_S69_L001,genome1,test,test_F09_S69_L001_R1_001.fastq.gz,test_F09_S69_L001_R2_001.fastq.gz,36.80243,-118.15371,test_F09_S69_L001,x
-test_G09,test_G09,genome1,test,test_G09_R1.fastq.gz,test_G09_R2.fastq.gz,37.32226,-118.42213,test_G09,x
-test_G09,test_G09_S81_L001,genome1,test,test_G09_S81_L001_R1_001.fastq.gz,test_G09_S81_L001_R2_001.fastq.gz,37.32226,-118.42213,test_G09_S81_L001,x
-test_H09,test_H09_S93_L001,genome1,test,test_H09_S93_L001_R1_001.fastq.gz,test_H09_S93_L001_R2_001.fastq.gz,37.32226,-118.42213,test_H09_S93_L001,x
-test_H09,test_H09,genome1,test,test_H09_R1.fastq.gz,test_H09_R2.fastq.gz,37.32226,-118.42213,test_H09,x
-test_A10,test_A10,genome1,test,test_A10_R1.fastq.gz,test_A10_R2.fastq.gz,37.32226,-118.42213,test_A10,x
-test_A10,test_A10_S10_L001,genome1,test,test_A10_S10_L001_R1_001.fastq.gz,test_A10_S10_L001_R2_001.fastq.gz,37.32226,-118.42213,test_A10_S10_L001,x
-test_B10,test_B10,genome1,test,test_B10_R1.fastq.gz,test_B10_R2.fastq.gz,37.3072308,-118.3900358,test_B10,x
-test_B10,test_B10_S22_L001,genome1,test,test_B10_S22_L001_R1_001.fastq.gz,test_B10_S22_L001_R2_001.fastq.gz,37.3072308,-118.3900358,test_B10_S22_L001,x
-test_C10,test_C10,genome1,test,test_C10_R1.fastq.gz,test_C10_R2.fastq.gz,37.3072308,-118.3900358,test_C10,x
-test_C10,test_C10_S34_L001,genome1,test,test_C10_S34_L001_R1_001.fastq.gz,test_C10_S34_L001_R2_001.fastq.gz,37.3072308,-118.3900358,test_C10_S34_L001,x
-test_D10,test_D10,genome1,test,test_D10_R1.fastq.gz,test_D10_R2.fastq.gz,37.3071268,-118.3721058,test_D10,x
-test_D10,test_D10_S46_L001,genome1,test,test_D10_S46_L001_R1_001.fastq.gz,test_D10_S46_L001_R2_001.fastq.gz,37.3071268,-118.3721058,test_D10_S46_L001,x
-test_E10,test_E10,genome1,test,test_E10_R1.fastq.gz,test_E10_R2.fastq.gz,37.3071268,-118.3721058,test_E10,x
-test_E10,test_E10_S58_L001,genome1,test,test_E10_S58_L001_R1_001.fastq.gz,test_E10_S58_L001_R2_001.fastq.gz,37.3071268,-118.3721058,test_E10_S58_L001,x
-test_F10,test_F10,genome1,test,test_F10_R1.fastq.gz,test_F10_R2.fastq.gz,37.3143988,-118.3719325,test_F10,x
-test_F10,test_F10_S70_L001,genome1,test,test_F10_S70_L001_R1_001.fastq.gz,test_F10_S70_L001_R2_001.fastq.gz,37.3143988,-118.3719325,test_F10_S70_L001,x
-test_G10,test_G10_S82_L001,genome1,test,test_G10_S82_L001_R1_001.fastq.gz,test_G10_S82_L001_R2_001.fastq.gz,37.3107449,-118.3719762,test_G10_S82_L001,x
-test_G10,test_G10,genome1,test,test_G10_R1.fastq.gz,test_G10_R2.fastq.gz,37.3107449,-118.3719762,test_G10,x
-test_H10,test_H10_S94_L001,genome1,test,test_H10_S94_L001_R1_001.fastq.gz,test_H10_S94_L001_R2_001.fastq.gz,37.143437,-118.24242,test_H10_S94_L001,x
-test_H10,test_H10,genome1,test,test_H10_R1.fastq.gz,test_H10_R2.fastq.gz,37.143437,-118.24242,test_H10,x
-test_A11,test_A11,genome1,test,test_A11_R1.fastq.gz,test_A11_R2.fastq.gz,37.143437,-118.24242,test_A11,x
-test_A11,test_A11_S11_L001,genome1,test,test_A11_S11_L001_R1_001.fastq.gz,test_A11_S11_L001_R2_001.fastq.gz,37.143437,-118.24242,test_A11_S11_L001,x
-test_B11,test_B11,genome1,test,test_B11_R1.fastq.gz,test_B11_R2.fastq.gz,37.143437,-118.24242,test_B11,x
-test_B11,test_B11_S23_L001,genome1,test,test_B11_S23_L001_R1_001.fastq.gz,test_B11_S23_L001_R2_001.fastq.gz,37.143437,-118.24242,test_B11_S23_L001,x
-test_C11,test_C11,genome1,test,test_C11_R1.fastq.gz,test_C11_R2.fastq.gz,36.80786,-118.16565,test_C11,x
-test_C11,test_C11_S35_L001,genome1,test,test_C11_S35_L001_R1_001.fastq.gz,test_C11_S35_L001_R2_001.fastq.gz,36.80786,-118.16565,test_C11_S35_L001,x
-test_D11,test_D11_S47_L001,genome1,test,test_D11_S47_L001_R1_001.fastq.gz,test_D11_S47_L001_R2_001.fastq.gz,36.80786,-118.16565,test_D11_S47_L001,x
-test_D11,test_D11,genome1,test,test_D11_R1.fastq.gz,test_D11_R2.fastq.gz,36.80786,-118.16565,test_D11,x
-test_E11,test_E11,genome1,test,test_E11_R1.fastq.gz,test_E11_R2.fastq.gz,36.80786,-118.16565,test_E11,x
-test_E11,test_E11_S59_L001,genome1,test,test_E11_S59_L001_R1_001.fastq.gz,test_E11_S59_L001_R2_001.fastq.gz,36.80786,-118.16565,test_E11_S59_L001,x
-test_F11,test_F11_S71_L001,genome1,test,test_F11_S71_L001_R1_001.fastq.gz,test_F11_S71_L001_R2_001.fastq.gz,37.4552086,-118.4039482,test_F11_S71_L001,x
-test_F11,test_F11,genome1,test,test_F11_R1.fastq.gz,test_F11_R2.fastq.gz,37.4552086,-118.4039482,test_F11,x
-test_G11,test_G11,genome1,test,test_G11_R1.fastq.gz,test_G11_R2.fastq.gz,37.4087202,-118.3819057,test_G11,x
-test_G11,test_G11_S83_L001,genome1,test,test_G11_S83_L001_R1_001.fastq.gz,test_G11_S83_L001_R2_001.fastq.gz,37.4087202,-118.3819057,test_G11_S83_L001,x
-test_H11,test_H11,genome1,test,test_H11_R1.fastq.gz,test_H11_R2.fastq.gz,37.4087202,-118.3819057,test_H11,x
-test_H11,test_H11_S95_L001,genome1,test,test_H11_S95_L001_R1_001.fastq.gz,test_H11_S95_L001_R2_001.fastq.gz,37.4087202,-118.3819057,test_H11_S95_L001,x
-test_A12,test_A12,genome1,test,test_A12_R1.fastq.gz,test_A12_R2.fastq.gz,37.4087202,-118.3819057,test_A12,x
-test_A12,test_A12_S12_L001,genome1,test,test_A12_S12_L001_R1_001.fastq.gz,test_A12_S12_L001_R2_001.fastq.gz,37.4087202,-118.3819057,test_A12_S12_L001,x
-test_B12,test_B12,genome1,test,test_B12_R1.fastq.gz,test_B12_R2.fastq.gz,37.4087202,-118.3819057,test_B12,x
-test_B12,test_B12_S24_L001,genome1,test,test_B12_S24_L001_R1_001.fastq.gz,test_B12_S24_L001_R2_001.fastq.gz,37.4087202,-118.3819057,test_B12_S24_L001,x
-test_C12,test_C12,genome1,test,test_C12_R1.fastq.gz,test_C12_R2.fastq.gz,37.3440208,-118.3455811,test_C12,x
-test_C12,test_C12_S36_L001,genome1,test,test_C12_S36_L001_R1_001.fastq.gz,test_C12_S36_L001_R2_001.fastq.gz,37.3440208,-118.3455811,test_C12_S36_L001,x
-test_D12,test_D12_S48_L001,genome1,test,test_D12_S48_L001_R1_001.fastq.gz,test_D12_S48_L001_R2_001.fastq.gz,37.3440208,-118.3455811,test_D12_S48_L001,x
-test_D12,test_D12,genome1,test,test_D12_R1.fastq.gz,test_D12_R2.fastq.gz,37.3440208,-118.3455811,test_D12,x
-test_E12,test_E12,genome1,test,test_E12_R1.fastq.gz,test_E12_R2.fastq.gz,37.3440208,-118.3455811,test_E12,x
-test_E12,test_E12_S60_L001,genome1,test,test_E12_S60_L001_R1_001.fastq.gz,test_E12_S60_L001_R2_001.fastq.gz,37.3440208,-118.3455811,test_E12_S60_L001,x
-test_F12,test_F12,genome1,test,test_F12_R1.fastq.gz,test_F12_R2.fastq.gz,37.2859424,-118.313121,test_F12,x
-test_F12,test_F12_S72_L001,genome1,test,test_F12_S72_L001_R1_001.fastq.gz,test_F12_S72_L001_R2_001.fastq.gz,37.2859424,-118.313121,test_F12_S72_L001,x
-test_G12,test_G12,genome1,test,test_G12_R1.fastq.gz,test_G12_R2.fastq.gz,37.2859424,-118.313121,test_G12,x
-test_G12,test_G12_S84_L001,genome1,test,test_G12_S84_L001_R1_001.fastq.gz,test_G12_S84_L001_R2_001.fastq.gz,37.2859424,-118.313121,test_G12_S84_L001,x
-test_H12,test_H12,genome1,test,test_H12_R1.fastq.gz,test_H12_R2.fastq.gz,37.2859424,-118.313121,test_H12,x
-test_H12,test_H12_S96_L001,genome1,test,test_H12_S96_L001_R1_001.fastq.gz,test_H12_S96_L001_R2_001.fastq.gz,37.2859424,-118.313121,test_H12_S96_L001,x
-testI_A01,testI_A01,genome1,test,testI_A01_R1.fastq.gz,testI_A01_R2.fastq.gz,37.2859424,-118.313121,testI_A01,x
-testI_B01,testI_B01,genome1,test,testI_B01_R1.fastq.gz,testI_B01_R2.fastq.gz,37.3178459,-118.3719749,testI_B01,x
-testI_C01,testI_C01,genome1,test,testI_C01_R1.fastq.gz,testI_C01_R2.fastq.gz,37.3178459,-118.3719749,testI_C01,x
-testI_D01,testI_D01,genome1,test,testI_D01_R1.fastq.gz,testI_D01_R2.fastq.gz,37.3178459,-118.3719749,testI_D01,x
-testI_E01,testI_E01,genome1,test,testI_E01_R1.fastq.gz,testI_E01_R2.fastq.gz,37.3178459,-118.3719749,testI_E01,x
-testI_F01,testI_F01,genome1,test,testI_F01_R1.fastq.gz,testI_F01_R2.fastq.gz,37.3071268,-118.3721058,testI_F01,x
-testI_G01,testI_G01,genome1,test,testI_G01_R1.fastq.gz,testI_G01_R2.fastq.gz,37.4319511,-118.3819852,testI_G01,x
-testI_H01,testI_H01,genome1,test,testI_H01_R1.fastq.gz,testI_H01_R2.fastq.gz,37.4319511,-118.3819852,testI_H01,x
-testI_A02,testI_A02,genome1,test,testI_A02_R1.fastq.gz,testI_A02_R2.fastq.gz,37.4319511,-118.3819852,testI_A02,x
-testI_B02,testI_B02,genome1,test,testI_B02_R1.fastq.gz,testI_B02_R2.fastq.gz,37.398646,-118.355055,testI_B02,x
-testI_C02,testI_C02,genome1,test,testI_C02_R1.fastq.gz,testI_C02_R2.fastq.gz,37.4481,-118.3945,testI_C02,x
-testI_D02,testI_D02,genome1,test,testI_D02_R1.fastq.gz,testI_D02_R2.fastq.gz,37.4481,-118.3945,testI_D02,x
-testI_E02,testI_E02,genome1,test,testI_E02_R1.fastq.gz,testI_E02_R2.fastq.gz,37.3179155,-118.4037828,testI_E02,x
-testI_F02,testI_F02,genome1,test,testI_F02_R1.fastq.gz,testI_F02_R2.fastq.gz,37.3179155,-118.4037828,testI_F02,x
-testI_G02,testI_G02,genome1,test,testI_G02_R1.fastq.gz,testI_G02_R2.fastq.gz,37.3293336,-118.4584599,testI_G02,x
-testI_H02,testI_H02,genome1,test,testI_H02_R1.fastq.gz,testI_H02_R2.fastq.gz,37.3293336,-118.4584599,testI_H02,x
-testI_A03,testI_A03,genome1,test,testI_A03_R1.fastq.gz,testI_A03_R2.fastq.gz,37.3293336,-118.4584599,testI_A03,x
-testI_B03,testI_B03,genome1,test,testI_B03_R1.fastq.gz,testI_B03_R2.fastq.gz,37.3293336,-118.4584599,testI_B03,x
-testI_C03,testI_C03,genome1,test,testI_C03_R1.fastq.gz,testI_C03_R2.fastq.gz,37.4019387,-118.3352434,testI_C03,x
-testI_D03,testI_D03,genome1,test,testI_D03_R1.fastq.gz,testI_D03_R2.fastq.gz,37.4019387,-118.3352434,testI_D03,x
-testI_E03,testI_E03,genome1,test,testI_E03_R1.fastq.gz,testI_E03_R2.fastq.gz,37.4019387,-118.3352434,testI_E03,x
-testI_F03,testI_F03,genome1,test,testI_F03_R1.fastq.gz,testI_F03_R2.fastq.gz,37.3873652,-118.3319068,testI_F03,x
-testI_G03,testI_G03,genome1,test,testI_G03_R1.fastq.gz,testI_G03_R2.fastq.gz,37.3945214,-118.3363653,testI_G03,x
-testI_H03,testI_H03,genome1,test,testI_H03_R1.fastq.gz,testI_H03_R2.fastq.gz,37.4089418,-118.3318361,testI_H03,x
-testI_A04,testI_A04,genome1,test,testI_A04_R1.fastq.gz,testI_A04_R2.fastq.gz,37.4089418,-118.3318361,testI_A04,x
-testI_B04,testI_B04,genome1,test,testI_B04_R1.fastq.gz,testI_B04_R2.fastq.gz,37.3945214,-118.3363653,testI_B04,x
-testI_C04,testI_C04,genome1,test,testI_C04_R1.fastq.gz,testI_C04_R2.fastq.gz,37.3944865,-118.3319682,testI_C04,x
-testI_D04,testI_D04,genome1,test,testI_D04_R1.fastq.gz,testI_D04_R2.fastq.gz,37.3944865,-118.3319682,testI_D04,x
-testI_E04,testI_E04,genome1,test,testI_E04_R1.fastq.gz,testI_E04_R2.fastq.gz,37.3615988,-118.5831825,testI_E04,x
-testI_F04,testI_F04,genome1,test,testI_F04_R1.fastq.gz,testI_F04_R2.fastq.gz,37.3615988,-118.5831825,testI_F04,x
-testI_G04,testI_G04,genome1,test,testI_G04_R1.fastq.gz,testI_G04_R2.fastq.gz,37.3615988,-118.5831825,testI_G04,x
-testI_H04,testI_H04,genome1,test,testI_H04_R1.fastq.gz,testI_H04_R2.fastq.gz,37.38029,-118.50387,testI_H04,x
-testI_A05,testI_A05,genome1,test,testI_A05_R1.fastq.gz,testI_A05_R2.fastq.gz,37.38029,-118.50387,testI_A05,x
-testI_B05,testI_B05,genome1,test,testI_B05_R1.fastq.gz,testI_B05_R2.fastq.gz,37.38029,-118.50387,testI_B05,x
-testI_C05,testI_C05,genome1,test,testI_C05_R1.fastq.gz,testI_C05_R2.fastq.gz,37.38029,-118.50387,testI_C05,x
-testI_D05,testI_D05,genome1,test,testI_D05_R1.fastq.gz,testI_D05_R2.fastq.gz,37.4089068,-118.3273518,testI_D05,x
-testI_E05,testI_E05,genome1,test,testI_E05_R1.fastq.gz,testI_E05_R2.fastq.gz,37.4089068,-118.3273518,testI_E05,x
-testI_F05,testI_F05,genome1,test,testI_F05_R1.fastq.gz,testI_F05_R2.fastq.gz,37.40084,-118.38065,testI_F05,x
-testI_G05,testI_G05,genome1,test,testI_G05_R1.fastq.gz,testI_G05_R2.fastq.gz,37.40084,-118.38065,testI_G05,x
-testI_H05,testI_H05,genome1,test,testI_H05_R1.fastq.gz,testI_H05_R2.fastq.gz,37.394555,-118.3273985,testI_H05,x
-testI_A06,testI_A06,genome1,test,testI_A06_R1.fastq.gz,testI_A06_R2.fastq.gz,37.394555,-118.3273985,testI_A06,x
-testI_B06,testI_B06,genome1,test,testI_B06_R1.fastq.gz,testI_B06_R2.fastq.gz,37.82036667,-118.4671333,testI_B06,x
-testI_C06,testI_C06,genome1,test,testI_C06_R1.fastq.gz,testI_C06_R2.fastq.gz,37.82036667,-118.4671333,testI_C06,x
-testI_D06,testI_D06,genome1,test,testI_D06_R1.fastq.gz,testI_D06_R2.fastq.gz,37.82036667,-118.4671333,testI_D06,x
-testI_E06,testI_E06,genome1,test,testI_E06_R1.fastq.gz,testI_E06_R2.fastq.gz,37.51016,-118.36735,testI_E06,x
-testI_F06,testI_F06,genome1,test,testI_F06_R1.fastq.gz,testI_F06_R2.fastq.gz,37.51016,-118.36735,testI_F06,x
-testI_G06,testI_G06,genome1,test,testI_G06_R1.fastq.gz,testI_G06_R2.fastq.gz,37.51016,-118.36735,testI_G06,x
-testI_H06,testI_H06,genome1,test,testI_H06_R1.fastq.gz,testI_H06_R2.fastq.gz,34.42766,-118.16865,testI_H06,x
-testI_A07,testI_A07,genome1,test,testI_A07_R1.fastq.gz,testI_A07_R2.fastq.gz,34.43413,-118.15831,testI_A07,x
-testI_B07,testI_B07,genome1,test,testI_B07_R1.fastq.gz,testI_B07_R2.fastq.gz,34.43413,-118.15831,testI_B07,x
-testI_C07,testI_C07,genome1,test,testI_C07_R1.fastq.gz,testI_C07_R2.fastq.gz,32.909121,-118.496472,testI_C07,x
-testI_D07,testI_D07,genome1,test,testI_D07_R1.fastq.gz,testI_D07_R2.fastq.gz,32.909121,-118.496472,testI_D07,x
-testI_E07,testI_E07,genome1,test,testI_E07_R1.fastq.gz,testI_E07_R2.fastq.gz,32.909121,-118.496472,testI_E07,x
-testI_F07,testI_F07,genome1,test,testI_F07_R1.fastq.gz,testI_F07_R2.fastq.gz,32.9078387,-118.5314941,testI_F07,x
diff --git a/config/test_qc_gls_config.yaml b/config/test_qc_gls_config.yaml
deleted file mode 100644
index d77b7cba..00000000
--- a/config/test_qc_gls_config.yaml
+++ /dev/null
@@ -1,75 +0,0 @@
-##############################
-# Variables you need to change
-##############################
-
-samples: "config/test_coords.csv" # name of the sample metadata CSV
-tmp_dir: "tmp/" # directory path for a temp dir
-split_by_n: True #set to False to split by chromosome/scaffold; set to True to split on runs of Ns within chromosomes/scaffolds.
-sentieon: False #set to True if you want to use sentieon, False if you want GATK
-sentieon_lic: "" #set to path of sentieon license
-remote_reads: False # set if you want reads to be on google cloud storage remote
-remote_reads_prefix: "" # set to google bucket name where reads live
-##############################
-# Variables you *might* need to change
-##############################
-
-# if using the BAM -> VCF workflows
-bamsForGatk: "fastq2bam/01_mappedReads/" # directory containing BAMs for GATK. If you used the fastq -> bam workflow above, simply keep the default value; must be followed by a "/"
-bamsForFB: "fastq2bam/01_mappedReads/" # directory containing BAMs for Freebayes. If you used the fastq -> bam workflow above, simply keep the default value; must be followed by a "/"
-bam_suffix: "_final.bam" # the suffix for your BAM files that follow all the sample names. If you use the fastq->BAM workflow above, simply keep the default value
-
-# These parameters control how the genome gets split into intervals according to Nmers in the reference genome
-maxNumIntervals: 10 # the maximum number of intervals when split_by_n is False. If your reference genome has hundreds of scaffolds, it can be helpful to set this to less than number of scaffolds.
-minNmer: 500 # the minimum Nmer used to split up the genome; e.g. a value of 200 means only Nmers 200 or greater are used to define the boundaries of intervals. The minimum is 50.
-maxIntervalLen: 15000000 # the desired maximum size of an interval for calling variants; more than 2Mb is a good starting point
-maxBpPerList: 15000000 # the desired maximum number of bp per list file for GATK4; list files potentially contain many small intervals, and we cap the fraction of the genome contained in each list file here
-maxIntervalsPerList: 200 # the desired maximum number of intervals per list file for GATK4; this prevents list files from containing thousands of small intervals, which can slow parts of GATK4. Default is good.
-maxDP_fb: 200 # this is the maximum depth parameter used for freebayes, site with more than this number are ignored, change according to expected depth
-
-## Coverage options ##
-## default pipeline is optimized for low coverage data - if using high coverage data (> 10x), uncomment high coverage options and comment out low coverage options
-
-# low coverage options (< 10x)
-minP: 1
-minD: 1
-
-# high coverage options (> 10x)
-#minP: 2
-#minD: 4
-
-## QC options ##
-nClusters: 3
-
-##############################
-# Variables you DO NOT need to change
-# if you do, ensure all path/directory names are followed by a "/"
-# These variables control the output of the fastq2bam workflow. You don't need to change these, but if you do please have a forward slash follow name!
-##############################
-
-output: "results/"
-fastqDir: "data/fastq/" # this is where raw fastqs are downloaded
-refGenomeDir: "data/genome/" # where reference genomes go
-fastq2bamDir: "fastq2bam/"
-fastqFilterDir: "00_fastqFiltered/" # new directory created for filtered fastq reads
-bamDir: "01_mappedReads/" # new directory created for BAM files
-sumstatDir: "02_bamSumstats/" # new directory created for BAM summary statistics
-
-# These variables control the output of the GATK4 workflow (please have forward slash follow name!)
-gatkDir: "gatk/"
-gvcfDir: "03_gvcfs/"
-dbDir: "04_genomicsDB/"
-vcfDir_gatk: "05_vcfs/"
-qcDir: "06_QC/"
-intDir: "intervalFiles/"
-
-# These variables control the output of the FreeBayes workflow (please have forward slash follow name!)
-fbDir: "freebayes/"
-intervalDir: "00_intervals/"
-vcfDir_fb: "01_vcfs_byinterval/"
-
-#information for read groups
-flowcell: "FLOWCELL"
-platform: "ILLUMINA"
-
-cluster_config: "profiles/slurm/cluster_config.yml"
-test_qc: True
diff --git a/docs/CCGPhorizontalblack.png b/docs/CCGPhorizontalblack.png
deleted file mode 100644
index 2142c945..00000000
Binary files a/docs/CCGPhorizontalblack.png and /dev/null differ
diff --git a/docs/CCGPhorizontalblue.jpeg b/docs/CCGPhorizontalblue.jpeg
deleted file mode 100644
index d4cb8aff..00000000
Binary files a/docs/CCGPhorizontalblue.jpeg and /dev/null differ
diff --git a/docs/ccgp_workflow.png b/docs/ccgp_workflow.png
deleted file mode 100644
index 481ea2f1..00000000
Binary files a/docs/ccgp_workflow.png and /dev/null differ
diff --git a/docs/conf.py b/docs/conf.py
new file mode 100644
index 00000000..5c894051
--- /dev/null
+++ b/docs/conf.py
@@ -0,0 +1,38 @@
+# Configuration file for the Sphinx documentation builder.
+
+# -- Project information
+
+project = 'snpArcher'
+copyright = '2023, Cade Mirchandani'
+author = 'Cade Mirchandani'
+
+release = '0.1'
+version = '0.1.0'
+
+# -- General configuration
+
+extensions = [
+ 'sphinx.ext.duration',
+ 'sphinx.ext.doctest',
+ 'sphinx.ext.autodoc',
+ 'sphinx.ext.autosummary',
+ 'sphinx.ext.intersphinx',
+ 'myst_parser'
+]
+
+intersphinx_mapping = {
+ 'python': ('https://docs.python.org/3/', None),
+ 'sphinx': ('https://www.sphinx-doc.org/en/master/', None),
+}
+intersphinx_disabled_domains = ['std']
+
+templates_path = ['_templates']
+myst_enable_extensions = [
+ "html_image"
+]
+# -- Options for HTML output
+
+html_theme = 'sphinx_rtd_theme'
+
+# -- Options for EPUB output
+epub_show_urls = 'footnote'
diff --git a/docs/datasets.md b/docs/datasets.md
new file mode 100644
index 00000000..ea09e7c0
--- /dev/null
+++ b/docs/datasets.md
@@ -0,0 +1,4 @@
+# Datasets Produced by snpArcher
+A number of resequencing datasets have been run with snpArcher generating consistent variant calls, available via [Globus](https://www.globus.org/) in the [Comparative Population Genomics Data collection](https://app.globus.org/file-manager?origin_id=a6580c44-09fd-11ee-be16-195c41bc0be4&origin_path=%2F). Details of data processing are described [here](https://www.biorxiv.org/content/10.1101/2023.06.22.546168v1). If you use any of these datasets in your projects, please cite both the snpArcher paper and the original data producers.
+
+If you would like to contribute datasets you have created using snpArcher, please get in touch!
\ No newline at end of file
diff --git a/docs/examples.md b/docs/examples.md
new file mode 100644
index 00000000..d16affcd
--- /dev/null
+++ b/docs/examples.md
@@ -0,0 +1,140 @@
+# Examples
+On this page you will find an example project scenario and how to setup and run it using snpArcher.
+
+In this example, we have 10 resequenced individuals we would like to generate variant calls for. We will cover creating the sample sheet, selecting config options, and running the workflow.
+## Directory structure
+First, let's setup our directories as suggested in our [executing](./executing.md#optional-directory-setup) instructions. Let's assume we are working in a directory called `workdir/`, and the snpArcher repository has already been cloned there. We have also already created the `snparcher` conda env as instructed in the [setup docs](./setup.md#environment-setup).
+
+1. Let's create a directory to organize this project and future ones, call it `projects`. Then, create a new directory for this project, we'll call it `secretarybird_reseq`.
+```
+.
+├── projects
+│ └── secretarybird_reseq
+└── snpArcher
+ └── ...
+```
+```{note}
+Not all files and directories are shown, only relevant ones.
+```
+2. Copy the snpArcher config directory `snpArcher/config` to `projects/secretarybird_reseq`:
+```
+.
+├── projects
+│ └── secretarybird_reseq
+│ └── config
+│ └── config.yaml
+└── snpArcher
+ └── ...
+```
+
+3. Assume we already have all our sequence data and reference genome on our system, stored in a different location `/storage/data`. We do not need to move the raw data to our project directory.
+```{note}
+We'll cover the cases using SRA data and refSeq genomes later on in this example.
+```
+## Sample sheet setup
+Now we need to setup our sample sheet to inform snpArcher of our samples and their metadata. You can use any editor to create the sheet, as long as it is a CSV file. We will save the sample sheet in our project's config directory: `projects/secretarybird_reseq/samples.csv`. Below is the final sample sheet that we will use going forward, with explanations of each column following.
+
+For a more comprehensive explanation of the sample sheet, please refer to [here](./setup.md#creating-a-sample-sheet) for more details.
+
+
+### Final sample sheet
+```
+BioSample,LibraryName,Run,fq1,fq2,lat,long
+bird_1,bird_1_lib,1,/storage/data/bird_1_R1.fq.gz,/storage/data/bird_1_R2.fq.gz,-8.758119,-36.280061
+bird_2,bird_2_lib,2,/storage/data/bird_2_R1.fq.gz,/storage/data/bird_2_R2.fq.gz,-72.336165,35.751903
+bird_3,bird_3_lib,3,/storage/data/bird_3_R1.fq.gz,/storage/data/bird_3_R2.fq.gz,-11.874137,-5.382251
+bird_4,bird_4_lib,4,/storage/data/bird_4_R1.fq.gz,/storage/data/bird_4_R2.fq.gz,-73.235723,-145.261219
+bird_5,bird_5_lib,5,/storage/data/bird_5_R1.fq.gz,/storage/data/bird_5_R2.fq.gz,88.08701,-52.658705
+bird_6,bird_6_lib,6,/storage/data/bird_6_R1.fq.gz,/storage/data/bird_6_R2.fq.gz,69.640536,-12.971862
+bird_7,bird_7_lib,7,/storage/data/bird_7_R1.fq.gz,/storage/data/bird_7_R2.fq.gz,18.608941,-100.485774
+bird_8,bird_8_lib,8,/storage/data/bird_8_R1.fq.gz,/storage/data/bird_8_R2.fq.gz,-36.570632,-102.38721
+bird_9,bird_9_lib,9,/storage/data/bird_9_R1.fq.gz,/storage/data/bird_9_R2.fq.gz,-88.592265,157.406505
+bird_10,bird_10_lib,10,/storage/data/bird_10_R1.fq.gz,/storage/data/bird_10_R2.fq.gz,40.106437,-58.649016
+```
+### Description of Columns
+1. **BioSample**: This is the name for the sample.
+2. **LibraryName**: Identifier for the sample's sequencing library. This is especially important if you have samples that were sequenced multiple times across multiple lanes, which is not the case in this example. See [here](./setup.md#handling-samples-with-more-than-one-pair-of-reads) for more details.
+3. **Run**: If we were using reads from the SRA, this is where the sample's SRR accession would go. However, since we have local data, this just has to be a unique value.
+4. **fq1**: Path to the first read pair. Absolute paths are recommended. If we were using SRA data, this column should be omitted.
+5. **fq1**: Path to the second read pair. Same note as fq1.
+6. **lat**: Decimal latitude for the sample, used to generate map in QC module output.
+6. **long**: Decimal longitude for the sample, used to generate map in QC module output.
+
+```{note}
+If your project has multiple genomes, you can add the refPath and refGenome columns.
+```
+
+## Config file setup
+Now that we've created our sample sheet, we need to edit the config file we copied earlier: `projects/secretarybird_reseq/config.yaml`. This file controls the main options for controlling snpArcher's outputs. Refer to the [setup section](./setup.md#configuring-snparcher) for more details.
+
+In our example we are using all of the default options. This will configure snpArcher to perform variant calling using GATK with the scatter-by-intervals approach. Also, we have set our reference genome name and path since we want to use the same genome for all samples in our sample sheet.
+
+```{yaml}
+samples: "config/samples.csv" # path to the sample metadata CSV
+final_prefix: "" # prefix for final output files
+intervals: True #Set to True if you want to perform variant calling using interval approach.
+sentieon: False #set to True if you want to use sentieon, False if you want GATK
+sentieon_lic: "" #set to path of sentieon license
+remote_reads: False # Set True if reads are in a location seperate from --default-remote-prefix.
+bigtmp: "" #Set to a path with lots of free space to use for commands that require large amounts of temp space; defaults to system tmpdir if empty
+cov_filter: True #set to True if you want to include coverage thresholds in the callable sites bed file (default uses mappability only)
+generate_trackhub: True #Set to true if you want to generate a Genome Browser Trackhub. Dependent on postprocessing module.
+trackhub_email: "hi@email.com"
+##############################
+# Variables you *might* need to change
+##############################
+
+# Set reference genome here if you would like to you use the same reference genome for all samples in sample sheet. See docs for more info.
+refGenome: "bird_genome" # Name for reference genome
+refPath: "/storage/data/bird.fa.gz"
+```
+
+## Profile setup
+Snakemake uses profile YAML files to specify commonly used command line arguments, so you don't have to remember all of the arguments you need. Read more about profiles [here](https://snakemake.readthedocs.io/en/stable/executing/cli.html#profiles). To specify a profile, you can use the `--workflow-profile` option when running Snakemake. snpArcher comes with two profiles, `default` and `slurm`, found in the `profiles` directory of the repository. The default profile will be used automatically when running Snakemake, even without the flag `--workflow-profile`.
+
+The profile also enables you to specify the compute resources any of snpArcher's rules can use. This is done via the YAML keys `default-resources`, `set-resources`, and `set-threads`. `default-resources` will apply to all rules, and `set-resources` can be applied to indiviudal rules, overriding what the default was set to. There is no way to set a default thread value.
+
+First, we will specify how many threads each rule can use. This is the same using the default or SLURM profile. Both profiles come with reasonable default thread values, but you may need to adjust based on your system or cluster.
+
+Let's say we wanted the alignment step (bwa mem) to use more threads:
+```{yaml}
+# ...
+set-threads:
+ bwa_map: 16 # Changed from 8 to 16.
+# ...
+```
+Next, we will specify memory and other resources. This step only applies if you are running on a SLURM cluster.
+
+In our example cluster, we have two compute partitions, "short" and "long". So we want to put long running jobs on the "long" partition, and the rest on "short". Additionally, the "short" partition has a timelimit of 1 hour and "long" 10 hours, so we will specify that.
+
+First, lets specify the default resources:
+```{yaml}
+default-resources:
+ mem_mb: attempt * 2000
+ mem_mb_reduced: (attempt * 2000) * 0.9 # Mem allocated to java for GATK rules (tries to prevent OOM errors)
+ slurm_partition: "short" # This line was changed
+ slurm_account: # Same as sbatch -A. Not all clusters use this.
+ runtime: 60 # In minutes
+```
+Then, lets modify the specific resources for the GATK HaplotypeCaller step:
+```{yaml}
+set-resources:
+# ... other rules
+ bam2gvcf: # HaplotypeCaller <--- This line was uncommented
+# mem_mb: attempt * 2000
+# mem_mb_reduced: (attempt * 2000) * 0.9 # Mem allocated to java (tries to prevent OOM errors)
+ slurm_partition: "long" # This line was changed
+ runtime: 600 # This line was changed
+```
+
+## Running the workflow
+We are now ready to run the workflow! From our working directory we can run the command:
+```
+snakemake -s snpArcher/workflow/Snakefile -d projects/secretarybird_reseq
+```
+This instructs Snakemake to use snpArcher's workflow file, and to run in the project directory we setup using the config and sample sheet we setup there. This will also use the default profile found in `snpArcher/profiles/default`.
+
+If we were on a SLURM cluster, we would specify the slurm profile:
+```
+snakemake -s snpArcher/workflow/Snakefile -d projects/secretarybird_reseq --workflow-profile snpArcher/profiles/slurm
+```
\ No newline at end of file
diff --git a/docs/executing.md b/docs/executing.md
new file mode 100644
index 00000000..836905f7
--- /dev/null
+++ b/docs/executing.md
@@ -0,0 +1,121 @@
+# Running snpArcher
+## Setup
+Please refer to our [setup instructions](./setup.md) to prepare the snpArcher environment and requisite files.
+## Test datasets
+To test that your environment is properly setup, you can run a quick test with the following command:
+```
+snakemake -d .test/ecoli --cores 1 --use-conda
+```
+If this runs without errors, you are ready to go!
+## Using the Dry-run option
+Snakemake offers the `--dry-run (-n)` CLI option to perform a dry-run of the workflow to show what jobs would be run. We recommend doing this before executing snpArcher to ensure that the sample sheet was setup correctly, and Snakemake has correctly built the workflow DAG.
+## Local Execution
+Once you have setup the requisite configuration files and sample sheet, executing snpArcher on your local machine is as simple as running the Snakemake command with the number of cores you would like to use. For example, to use 8 cores you would run:
+```
+snakemake --cores 8 --use-conda
+```
+
+### Optional directory setup
+To maintain organization across many different projects, you may consider creating a new directory for each project you run using snpArcher. This way, each of your project directories will contain the configuration files used for that run. Below is an example directory structure:
+
+```
+.
+├── snpArcher
+├── project_1/
+│ ├── config/
+│ │ ├── config.yaml
+│ │ └── samples.csv
+│ ├── data
+│ └── results
+└── project_2/
+ ├── config/
+ │ ├── config.yaml
+ │ └── samples.csv
+ └── data
+```
+
+When creating a new directory for an analysis, ensure that you copy the `config` directory from the snpArcher directory to your new directory.
+
+Then, to run snpArcher on `project_2` from our example, we would execute the command:
+```
+snakemake -s ./snpArcher/workflow/Snakefile -d ./project_2 
+
+Have resequencing data as fastq files and a reference genome? Want a VCF file of genotypes? Use snpRarcher as your one-stop shop to quickly and efficiently produce an analysis-ready dataset. No need to create a hand-tailored workflow cobbled together by tape and error-ridden chatGPT code — use snpArcher for all your variant calling needs on your laptop, on your server, or up in the clouds.
+
+snpArcher is a reproducible workflow optimized for nonmodel organisms and comparisons across datasets, built on the [Snakemake](https://snakemake.readthedocs.io/en/stable/index.html#) workflow management system. It provides a streamlined approach to dataset acquisition, variant calling, quality control, and downstream analysis.
+
+Snakemake makes it easy to bundle together the many steps involved in running a bioinformatics pipeline. The workflow involves mapping reads to a reference genome, calling SNPs using GATK's haplotypecaller, and calling variants at the population level using GATK's combineGVCFs. Each of these steps can be slow and tiresome to run on their own and the workflow has been carefully designed and tested to maximize efficiency. We use intervals to break up jobs into smaller chunks so that time and memory-hungry steps like haplotypecaller run quickly. If you have access to a Sentieon license for accelerated variant calling, we include options for this.
+
+Finally, the pipeline makes it easy to evaluate how the data looks. Review the HTML file in the QC folder at the end of a run to see how your samples relate to each other and also a number of metrics for evaluating variant-calling quality.
+
+Remember to examine the config.yaml file to edit options for each step. We have carefully chosen default options that should work for most users, but these can be tweaked here.
+
+## Requirements
+- Illumina paired-end fastq files for one or more individuals
+- A reference genome
+- A sample sheet with sample names matched to the read names
+- Snakemake and Mamba installed on your system
+- If using Google Cloud, you will need to have set up an account on the GCP console
+
+## Using snpArcher
+- To get started quickly, check out the quick start tutorial!
+- Otherwise start [here](./setup.md).
+
+## Citing
+- Please cite our preprint [here](https://www.biorxiv.org/content/10.1101/2023.06.22.546168v1)
+- Also, make sure to cite the tools you used within snpArcher.
+
+## Contributing to snpArcher
+- If you encounter a bug or want to request a feature, please open a issue on our [github page](https://github.com/harvardinformatics/snpArcher).
+- If you'd like to contribute a module, check out our [module contribution guidelines](./modules.md#module-contribution-guidelines).
+
+```{toctree}
+:hidden: True
+./setup.md
+./executing.md
+./examples.md
+./modules.md
+./datasets.md
+```
diff --git a/docs/modules.md b/docs/modules.md
new file mode 100644
index 00000000..e21a49bb
--- /dev/null
+++ b/docs/modules.md
@@ -0,0 +1,47 @@
+# Modules
+A key goal in the design of the snpArcher pipeline is to allow seamless extensibility with downstream processing. We implement this using Snakemake modules, which allow additional rules to easily extend the main pipeline. We present several modular extensions of snpArcher here, but we hope also that user-developed modules will grow the set of tools linked to snpArcher in order to facilitate diverse analysis.
+## Module Contribution Guidelines
+We developed a set of criteria for including additional user-developed modules into snpArcher. This project is designed to be modular and easily extensible as we and workflow users develop additional features and downstream analysis tools. To ensure that contributed modules are reproducible and easily implemented, we propose the following evaluation criteria:
+
+1. Each module must include Snakemake workflow that defines necessary environments using Conda.
+2. The module must be freely distributed via Github with sufficient documentation that users can adapt it to their needs
+3. The module must provide a unit test based on either existing test datasets available from the main snpArcher workflow or via a module-specific minimal test dataset
+4. Each module should be registered within the main project page to enhance discoverability and ensure the above criteria are met.
+
+If you are interested in developing a module please reach out via email or Github, we'd love to know and chat about it.
+## Quality Control
+The quality control module aggregates various statistics from the workflow and produces preliminary analyses and plots in an interactive HTML file, offering visualizations of summary statistics related to population structure, batch effects, sequencing depth, genetic relatedness, geography, and admixture. Most summaries are based on a random sample of 100,000 SNPs, while others provide high-level summaries of the full variant dataset. These visualizations help identify outliers, potential biases, and sequencing artifacts, such as cryptic genetic variation, batch effects, and reference bias. Additionally, an interactive heatmap aids in quickly identifying close relatives within the dataset, and spatial maps provide a visualization of PCA clusters in space.
+### Config Options
+| Option | Description | Type |
+| ---- | -------------| ------ |
+|`nClusters`| Number of clusters for PCA| `int`|
+|`GoogleAPIKey`| Google Maps API key (optional).| `str`|
+|`min_depth`| Samples with average depth below this will be excluded for QC analysis| `int`|
+
+```{note}
+To generate the QC dashboard, you must have at least 2 samples specified in your sample sheet.
+```
+```{note}
+The output of the QC module should not be considered a final analysis and is solely intended to direct quality control of the dataset.
+```
+## Postprocessing
+The postprocessing module is designed to be run after snpArcher has intially been run and you have determined if there are samples that you would like to exclude from downstream analyses. In order to trigger this module, you must add the `SampleType` column to your sample sheet, and mark samples for inclusion with the value `include` and exclusion with the value `exclude`.
+
+This module produces a filtered VCF by filtering excluded samples as well as sites not passing default and user defined filters.
+### Config Options
+| Option | Description | Type |
+| ---- | -------------| ------ |
+|`contig_size`| SNPs on contigs this size or smaller will be excluded from 'clean' VCF | `int`|
+|`maf`| SNPs with MAF below this will be excluded from clean VCF| `float`|
+|`missingness`| SNPs with missingness below this will be excluded from clean VCF| `float`|
+|`scaffolds_to_exclude` | Comma separated, no spaces list of scaffolds/contigs to exclude from clean VCF|
+
+```{hint}
+If you'd like to run the postprocessing module by default, you can add the `SampleType` column in your sample sheet, and mark all samples as `include`.
+```
+## Trackhubs
+The trackhub module generates UCSC Genome Browser track files to explore population variation data from the VCF produced by snpArcher. This module computes and generates genome browser tracks for traditional population genomic summary statistics such as windowed estimates of Tajima’s D, SNP density, Pi, Minor Allele Frequency, SNP depth. To trigger this module, you must set the [config](./setup.md#core-configuration) option to `True` and supply a email (a requirement for tracks displayed on the UCSC Genome Browser).
+
+```{warning}
+The Trackhubs module is dependent on the postprocessing module.
+```
diff --git a/docs/requirements.txt b/docs/requirements.txt
new file mode 100644
index 00000000..82621bec
--- /dev/null
+++ b/docs/requirements.txt
@@ -0,0 +1,5 @@
+# Defining the exact version will make sure things don't break
+sphinx==5.3.0
+sphinx_rtd_theme==1.1.1
+readthedocs-sphinx-search==0.1.1
+myst-parser==1.0.0
diff --git a/docs/setup.md b/docs/setup.md
new file mode 100644
index 00000000..65c7c4e4
--- /dev/null
+++ b/docs/setup.md
@@ -0,0 +1,195 @@
+# Setting up snpArcher
+## Environment Setup
+First, you will need to have [Mamba](https://mamba.readthedocs.io/en/latest/mamba-installation.html#mamba-install) installed. Follow the link and use the "Fresh Install (recommended)" directions.
+
+Mamba is a faster version of conda. Conda is a package manager that makes it easy to create local environments with pre-configured versioning for your favorite packages.
+
+Once Mamba is installed, create a conda environment with snakemake. These are the only two dependencies you need for the pipeline to work, the workflow will create mamba environments for each rule, and there is no need to install each package separately.
+
+```
+mamba create -c conda-forge -c bioconda -n snparcher "snakemake>=8" "python==3.11.4"
+mamba activate snparcher
+```
+If you encounter issues, please see the [Snakemake docs](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html) for detailed installation instructions.
+
+Next, clone the [snpArcher github repo](https://github.com/harvardinformatics/snpArcher) to your machine:
+```
+git clone https://github.com/harvardinformatics/snpArcher.git
+cd snpArcher
+```
+
+## Creating a sample sheet
+In order to determine what outputs to create, snpArcher requires sample sheet file. This comma-separated file contains the required sample metadata about your samples in order to run the workflow. At a minimum, the snpArcher pipeline requires that each sample have a unique sample name, a reference genome accession or a path to a fasta file, and an SRA accession, or path to two paired end fastq files.
+
+Below are all of the accepted fields for a sample sheet:
+| Field | Description |
+| ---- | -------------|
+| BioSample | The name of the sample. This will be the sample name in the final VCF |
+| LibraryName | LibraryID for sample, this can be the same or different than BioSample |
+| Run | The SRR for the sample, if applicable. If not, must be some **unique** value. It is often the lane number if samples are sequenced on multiple lanes. |
+| refGenome | Reference genome accession, if applicable. *See note* |
+| refPath | Path to local reference genome, if applicable. *See note* |
+| BioProject | If applicable. Otherwise any value is acceptable. |
+| fq1 | Optional if no SRR value in Run. Path to read 1 for sample |
+| fq2 | Optional if no SRR value in Run. Path to read 2 for sample |
+| SampleType | Optional. Triggers postproccesing module. Accepted values are 'include' or 'exclude' |
+| lat | Optional. Decimal latitude for sample, required to generate map in QC dashboard. |
+| long | Optional. Decimal longitude for sample, required to generate map in QC dashboard. |
+
+```{note}
+refGenome is always required. refPath specifying the path to a reference fasta file is optional, but when specified, a name for the assembly (in refGenome) must also be included.
+
+If you are using the same reference genome for all samples in your sample sheet, you can omit the refGenome and/or refPath column from the sample sheet and specify these fields in the config file. See [config setup below](#configuring-snparcher) for more details.
+```
+
+It is important to note that samples are proccessed together based on their `refGenome` metadata, so **all BioSamples that share a reference genome will ultimately end up in the same final vcf file.** If you are mapping multiple populations / species to a single reference genome, and want separate VCF files for each population / species, you will need to split your final vcf after the pipeline completes, or run multiple indpendent sample sheets in different results directories.
+
+If your reads (and, optionally, your local reference genome) are stored in somewhere seperate of the workflow (e.g.: a scratch disk) then you can specify the path to your reads using the `fq1` and `fq2` fields, and the location of your reference genome fasta in the `refPath` field.
+
+### Using data from NCBI SRA
+If you'd like to reanalyze an existing NCBI SRA BioProject, please follow these instructions to quickly create a sample sheet.
+
+1. Go to the BioProject overview web page on the SRA.
+2. In the subheading `Project Data` there is a table with the columns `Resource Name` and `Number of Links`. Click the link in the `Number of Links` column in the `SRA Experiments` row. You will be redirected to a search results page.
+3. Near the top of the search results page, click the link `Send results to Run Selector`
+4. On the Run Selector page, you can select/deselect samples you'd like to include/exclude in your sample sheet by using the checkboxes.
+5. Once you are done selecting samples, click the `Metadata` button in the `Download` column in the table near the middle of the page.
+6. This will download a a comma separated file called `SraRunTable.txt`.
+7. Open the file in the editor of your choice, and add a column named `refGenome`. In this column, enter the reference genome accession you want to use for every row in the sheet.
+8. Save the sample sheet, it is now ready to use.
+
+### Using local data
+A python script `workflow/snparcher_utils/write_samples.py` is included to help write the sample sheet for you. In order to use this script, you must have organized all of your fastq files in to one directory. The script requies you provide a file with one sample per name that maps uniquely to a pair of fastq files in the afformentioned directory. The script also requires either a reference genome accession or path to reference fasta.
+
+```{note}
+This script cannot currently handle multiple sequencing runs per sample. Please see below for how to handle this case.
+```
+
+Usage details:
+
+|Argument| Description|
+| ------ | ---------- |
+| `-s / --sample_list` | Path to a sample list. One sample per line |
+| `-f / --fastq_dir` | Path to directory containing ALL fastq files. It is assumed that each fastq file will contain the sample name uniquely. |
+| `-r / --ref` | Path to reference fasta. Mutually exclusive with -a|
+| `-a / --acc` | NCBI accession of reference. Mutually exclusive with -r|
+
+#### Handling samples with more than one pair of reads
+
+In order to specify samples that were sequenced multiple times in your sample sheet, you must:
+1. Create a duplicate row for each unit of sequencing
+2. Ensure the `BioSample` value is the same across all rows for the sample.
+3. Give each row a unique `Run` value. This allows snpArcher to collect all read pairs for a `BioSample`. All runs for a sample will be mapped separately to the genome and subsequently merged.
+4. Give each row a unique `LibraryName` value, if applicable. Used for marking duplicates, `LibraryName` should be the same in cases where the same library was sequenced multiple times. If a sample had multiple libraries prepared for it, then `LibraryName` should be unique for each library. See [here](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) for more info.
+
+For example, consider we have 2 samples: `A` and `B`. `Sample A` was sequenced 3 times, 2 of which were derived from the same library prep, and the other was a unique library. `Sample B `was only sequenced once. Below is how the sample sheet would look in order to define these relationships. Note, only the relevant fields have been included.
+
+| BioSample | LibraryName | Run |
+| --------- | ----------- | --- |
+| sample_A | lib_A_1 | 1 |
+| sample_A | lib_A_1 | 2 |
+| sample_A | lib_A_2 | 3 |
+| sample_B | lib_B_1 | 4 |
+
+
+## Configuring snpArcher
+
+Workflow variables such as output file prefix, tool settings, and other options are set in `config/config.yaml`. Resource settings such as threads and memory are controlled per tool in the `profiles/default/config.yaml`.
+
+### Core configuration
+The following options in `config/config.yaml` must be set before running snpArcher:
+
+| Option | Description | Type | Required | Default |
+| ---- | -------------| ------ |------ | ------ |
+| `samples` | Path to CSV sample sheet.| `str` | `True` | `None` |
+| `final_prefix` | Prefix to name final output files with (e.g. VCF) | `str` | `True` | `None` |
+| `intervals` | Use SplitByN interval approach for GATK variant calling | `bool` | `True` | `True` |
+| `sentieon` | Use Sentieon tools instead of GATK for variant calling | `bool` | `True` | `False` |
+| `sentieon_lic` | If using Sentieon tools, provide license here | `str` | `True` if `sentieon==True`| `None` |
+| `remote_reads` | Use remote storage provider for reads. | `bool` | `False` | `False`|
+| `bigtmp` | Set a directory for TMP. Default is $TMPDIR env var | `str` | `False` | `None` |
+| `cov_filter` | Use coverage thresholds for filtering callable sites | `bool`| `True` | `True` |
+| `generate_trackhub` | Generate population genomics stats trackhub | `bool`| `True` | `True` |
+| `trackhub_email` | Trackhubs require an email address | `str` | `True` if `generate_trackhub==True` | `None` |
+| `refGenome` | Reference genome name or accession | `str` | `True` if not provided in sample sheet | `None` |
+| `refPath` | Path to reference genome if not using NCBI genome accession | `str` | `False` | `None` |
+
+
+### Other options
+The following options can be adjusted based on your needs and your dataset.
+
+#### Variant Calling Options
+| Option | Description | Type |
+| ---- | -------------| ------ |
+|`minNmer`| The minimum span of Ns to split reference genome at for interval generation | `int`|
+|`num_gvcf_intervals` | The maximum number of GVCF intervals to create. Actual number of intervals may be less if reference genome is highly contiguous. | `int`|
+|`db_scatter_factor` | Used to calculate number of DB intervals to create. `num_db_intervals = (scatter_factor * num_samples * num_gvcf_intervals)`. Recommend <1 | `float`|
+| `minP` | Controls `--min-pruning` in GATK HaplotypeCaller. Recommend 1 for low coverage (<10x), 2 for high coverage (>10x) | `int` |
+| `minD` | Controls `--min-dangling-branch-length` in GATK HaplotypeCaller. Recommend 1 for low coverage (<10x), 4 for high coverage (>10x) | `int` |
+| `ploidy` | Ploidy for variant calling step. | `int` |
+
+#### Callable Sites Options
+| Option | Description | Type |
+| ---- | -------------| ------ |
+|`mappability_min`| Genomic regions with mappability score less than this will be removed from callable sites. | `int`|
+|`mappability_k`| Kmer size to compute mappability. | `int`|
+|`mappability_merge`| Merge passing mappability regions separated by this or fewer basepairs into a signle region | `int`|
+|`cov_merge`| Merge passing coverage regions separated by this or fewer basepairs into a signle region | `int`|
+
+#### Coverage Filtering Options
+If `cov_filter` in the core options is set to `True`, then the following options can be adjusted to the user's needs. Coverage filtering can be handled 3 ways:
+
+1. Hard upper and lower thresholds: regions with a mean coverage that falls within these thresholds are considered callable.
+
+| Option | Description | Type |
+| ---- | -------------| ------ |
+|`cov_threshold_lower`| Lower coverage threshold| `int`|
+|`cov_threshold_upper`| Upper coverage threshold| `int`|
+
+2. Standard deviations: regions with a mean coverage that is within N standard deviations (assumes Poisson distribution) are considered callable.
+
+| Option | Description | Type |
+| ---- | -------------| ------ |
+|`cov_threshold_stdev`| Number of standard deviations is considered callable | `int`|
+
+3. Absolute scaling: Thresholds set by factor N. Lower bowund is (global mean coverage / N), upper bound (global mean coverage * N). A region is callable if its mean coverage is within these bounds.
+
+| Option | Description | Type |
+| ---- | -------------| ------ |
+|`cov_threshold_rel`| Scaling factor for coverage threshold| `int`|
+
+```{note}
+In order to use one of the above coverage filtering approaches, you must set the options of the desired approach, and leave the others blank.
+```
+
+#### Module Options
+Please refer to the [modules page](./modules.md) for each module's options.
+
+### Resources
+Compute resources (threads, memory, etc.) as well as Snakemake arguments are set by the workflow profile located in `profiles/default/config.yaml`. This profile is used by default when running Snakemake. To specify a different profile, use the `--workflow-profile` option in your Snakemake command.
+
+In the profile you can set resources to be applied to all rules via the `default-resources` key. You can override this default per-rule by uncommenting that rule under the `set-threads` and/or `set-resources` key.
+
+#### Threads
+The profile controls how many `threads` (or CPU cores) a rule can use via the `set-thread` key. We have provided reasonable defaults, though you may need to adjust depending on the resources available on your system/cluster.
+```{note}
+Many rules can only use 1 thread, and providing more threads **will not** decrease runtime or improve performance. Please refer to the `profiles/default/config.yaml` for details.
+```
+
+#### Memory and other resources
+The profile controls how much memory and what other resources a rule can use via the `set-resources` key. When executing snpArcher on a cluster or the cloud, specifying memory is important as these environments will typically kill jobs that use more memory than they requested.
+
+Other resources, such as `slurm_partition`, `runtime`, etc. can also be set here if they are required by your cluster. We have provided a SLURM profile `profiles/slurm` that has the most common SLURM resources.
+
+```{note}
+Snakemake allows you to dynamically assign resources. We use the `attempt` keyword to specify memory. For example. `attempt * 2000` will provide 2GB on the first attempt of the rule, if the rule fails (out of memory) then on the second attempt it will be provided 4GB. This behavior requires the `-T/--retries` Snakemake option.
+```
+
+
+
+
+
+
+
+
+
diff --git a/profiles/default/config.yaml b/profiles/default/config.yaml
new file mode 100644
index 00000000..7ef1eebd
--- /dev/null
+++ b/profiles/default/config.yaml
@@ -0,0 +1,235 @@
+use-conda: True
+
+# These resources will be applied to all rules. Can be overriden on a per-rule basis below.
+default-resources:
+ mem_mb: attempt * 2000
+ mem_mb_reduced: (attempt * 2000) * 0.9 # Mem allocated to java for GATK rules (tries to prevent OOM errors)
+
+# Control number of threads each rule will use.
+set-threads:
+ # Reference Genome Processing. Does NOT use more than 1 thread.
+ download_reference: 1
+ index_reference: 1
+ # Interval Generation. Does NOT use more than 1 thread.
+ format_interval_list: 1
+ create_gvcf_intervals: 1
+ create_db_intervals: 1
+ picard_intervals: 1
+ # Mappability
+ genmap: 1 # Can use more than 1 thread
+ mappability_bed: 1 # Does NOT use more than 1 thread
+ # Fastq Processing. Can use more than 1 thread.
+ get_fastq_pe: 1
+ fastp: 1
+ # Alignment. Can use more than 1 thread, except merge_bams.
+ bwa_map: 1
+ dedup: 1
+ merge_bams: 1 # Does NOT use more than 1 thread.
+ # GVCF
+ bam2gvcf: 1 # Should be run with no more than 2 threads.
+ concat_gvcfs: 1 # Does NOT use more than 1 thread.
+ create_db_mapfile: 1 # Does NOT use more than 1 thread.
+ gvcf2DB: 1 # Should be run with no more than 2 threads.
+ # VCF
+ DB2vcf: 1 # Should be run with no more than 2 threads.
+ filterVcfs: 1 # Should be run with no more than 2 threads.
+ sort_gatherVcfs: 1 # Should be run with no more than 2 threads.
+ # Callable Bed
+ compute_d4: 1 # Can use more than 1 thread
+ create_cov_bed: 1 # Does NOT use more than 1 thread.
+ merge_d4: 1 # Does NOT use more than 1 thread.
+ # Summary Stats Does NOT use more than 1 thread.
+ bam_sumstats: 1
+ collect_covstats: 1
+ collect_fastp_stats: 1
+ collect_sumstats: 1
+ # QC Module Does NOT use more than 1 thread.
+ qc_admixture: 1
+ qc_check_fai: 1
+ qc_generate_coords_file: 1
+ qc_plink: 1
+ qc_qc_plots: 1
+ qc_setup_admixture: 1
+ qc_subsample_snps: 1
+ qc_vcftools_individuals: 1
+ # MK Module Does NOT use more than 1 thread.
+ mk_degenotate: 1
+ mk_prep_genome: 1
+ mk_split_samples: 1
+ # Postprocess Module Does NOT use more than 1 thread.
+ postprocess_strict_filter: 1
+ postprocess_basic_filter: 1
+ postprocess_filter_individuals: 1
+ postprocess_subset_indels: 1
+ postprocess_subset_snps: 1
+ postprocess_update_bed: 1
+ # Trackhub Module Does NOT use more than 1 thread.
+ trackhub_bcftools_depth: 1
+ trackhub_bedgraph_to_bigwig: 1
+ trackhub_calc_pi: 1
+ trackhub_calc_snpden: 1
+ trackhub_calc_tajima: 1
+ trackhub_chrom_sizes: 1
+ trackhub_convert_to_bedgraph: 1
+ trackhub_strip_vcf: 1
+ trackhub_vcftools_freq: 1
+ trackhub_write_hub_files: 1
+ # Sentieon Tools. Can use more than 1 thread, except sentieon_bam_stats.
+ sentieon_map: 1
+ sentieon_dedup: 1
+ sentieon_haplotyper: 1
+ sentieon_combine_gvcf: 1
+ sentieon_bam_stats: 1 # Does NOT use more than 1 thread.
+
+# # Control other resources used by each rule.
+# set-resources:
+# # Reference Genome Processing
+# copy_reference:
+# mem_mb: attempt * 200
+# download_reference:
+# mem_mb: attempt * 200
+# index_reference:
+# mem_mb: attempt * 200
+
+# # Interval Generation
+# format_interval_list:
+# mem_mb: attempt * 200
+# create_gvcf_intervals:
+# mem_mb: attempt * 200
+# create_db_intervals:
+# mem_mb: attempt * 200
+# picard_intervals:
+# mem_mb: attempt * 200
+
+# # Mappability
+# genmap:
+# mem_mb: attempt * 200
+# mappability_bed:
+# mem_mb: attempt * 200
+
+# # Fastq Processing
+# get_fastq_pe:
+# mem_mb: attempt * 200
+# fastp:
+# mem_mb: attempt * 200
+
+# # Alignment
+# bwa_map:
+# mem_mb: attempt * 200
+# dedup:
+# mem_mb: attempt * 200
+# merge_bams:
+# mem_mb: attempt * 200
+
+# # GVCF
+# bam2gvcf: # HaplotypeCaller
+# mem_mb: attempt * 200 # Mem allocated to the snakemake job
+# mem_mb_reduced: (attempt * 200) * 0.9 # Mem allocated to java (tries to prevent OOM errors)
+# concat_gvcfs:
+# mem_mb: attempt * 200
+# create_db_mapfile:
+# mem_mb: attempt * 200
+# gvcf2DB: # GenomicsDBImport
+# mem_mb: attempt * 200
+# mem_mb_reduced: (attempt * 200) * 0.9 # Mem allocated to java (tries to prevent OOM errors)
+
+# # VCF
+# DB2vcf: # GenotypeGVCFs
+# mem_mb: attempt * 200
+# mem_mb_reduced: (attempt * 200) * 0.9 # Mem allocated to java (tries to prevent OOM errors)
+# filterVcfs:
+# mem_mb: attempt * 200
+# sort_gatherVcfs:
+# mem_mb: attempt * 200
+
+# # Callable Bed
+# compute_d4:
+# mem_mb: attempt * 200
+# create_cov_bed:
+# mem_mb: attempt * 200
+# merge_d4:
+# mem_mb: attempt * 200
+
+# # Summary Stats
+# bam_sumstats:
+# mem_mb: attempt * 200
+# collect_covstats:
+# mem_mb: attempt * 200
+# collect_fastp_stats:
+# mem_mb: attempt * 200
+# collect_sumstats:
+# mem_mb: attempt * 200
+
+# # QC Module
+# qc_admixture:
+# mem_mb: attempt * 200
+# qc_check_fai:
+# mem_mb: attempt * 200
+# qc_generate_coords_file:
+# mem_mb: attempt * 200
+# qc_plink:
+# mem_mb: attempt * 200
+# qc_qc_plots:
+# mem_mb: attempt * 200
+# qc_setup_admixture:
+# mem_mb: attempt * 200
+# qc_subsample_snps:
+# mem_mb: attempt * 200
+# qc_vcftools_individuals:
+# mem_mb: attempt * 200
+
+# # MK Module
+# mk_degenotate:
+# mem_mb: attempt * 200
+# mk_prep_genome:
+# mem_mb: attempt * 200
+# mk_split_samples:
+# mem_mb: attempt * 200
+
+# # Postprocess Module
+# postprocess_strict_filter:
+# mem_mb: attempt * 200
+# postprocess_basic_filter:
+# mem_mb: attempt * 200
+# postprocess_filter_individuals:
+# mem_mb: attempt * 200
+# postprocess_subset_indels:
+# mem_mb: attempt * 200
+# postprocess_subset_snps:
+# mem_mb: attempt * 200
+# postprocess_update_bed:
+# mem_mb: attempt * 200
+
+# # Trackhub Module
+# trackhub_bcftools_depth:
+# mem_mb: attempt * 200
+# trackhub_bedgraph_to_bigwig:
+# mem_mb: attempt * 200
+# trackhub_calc_pi:
+# mem_mb: attempt * 200
+# trackhub_calc_snpden:
+# mem_mb: attempt * 200
+# trackhub_calc_tajima:
+# mem_mb: attempt * 200
+# trackhub_chrom_sizes:
+# mem_mb: attempt * 200
+# trackhub_convert_to_bedgraph:
+# mem_mb: attempt * 200
+# trackhub_strip_vcf:
+# mem_mb: attempt * 200
+# trackhub_vcftools_freq:
+# mem_mb: attempt * 200
+# trackhub_write_hub_files:
+# mem_mb: attempt * 200
+
+# # Sentieon Tools
+# sentieon_map:
+# mem_mb: attempt * 200
+# sentieon_dedup:
+# mem_mb: attempt * 200
+# sentieon_haplotyper:
+# mem_mb: attempt * 200
+# sentieon_combine_gvcf:
+# mem_mb: attempt * 200
+# sentieon_bam_stats:
+# mem_mb: attempt * 200
diff --git a/profiles/gls-gatk/config.yaml b/profiles/gls-gatk/config.yaml
deleted file mode 100644
index eb8d958e..00000000
--- a/profiles/gls-gatk/config.yaml
+++ /dev/null
@@ -1,8 +0,0 @@
-google-lifesciences: True
-google-lifesciences-region : 'us-central1'
-use-conda: True
-jobs: 150
-preemption-default: 5
-preemptible-rules: ["fastp=5", "bwa_map=5","merge_bams=5","dedup=5","bam2gvcf=10","create_db_mapfile=5","gvcf2DB=5","DB2vcf=5","filterVcfs=5","sort_gatherVcfs=5","qc_check_fai=5","qc_vcftools_individuals=5","qc_subsample_snps=5","qc_plink=5","qc_admixture=5","qc_generate_coords_file=5","qc_qc_plots=5", "collect_sumstats=2"]
-default-resources: ["machine_type=n2-standard", "disk_mb=512000"]
-
diff --git a/profiles/gls-sentieon/config.yaml b/profiles/gls-sentieon/config.yaml
deleted file mode 100644
index c4263a5a..00000000
--- a/profiles/gls-sentieon/config.yaml
+++ /dev/null
@@ -1,8 +0,0 @@
-google-lifesciences: True
-google-lifesciences-region : 'us-central1'
-use-conda: True
-jobs: 150
-preemption-default: 5
-preemptible-rules: ["fastp=5", "sentieon_map=5","merge_bams=5","sentieon_dedup=5","sentieon_haplotyper=10","sentieon_combine_gvcf=5","qc_check_fai=5","qc_vcftools_individuals=5","qc_subsample_snps=5","qc_plink=5","qc_admixture=5","qc_generate_coords_file=5","qc_qc_plots=5", "collect_sumstats=2"]
-default-resources: ["machine_type=n2-standard", "disk_mb=512000"]
-
diff --git a/profiles/slurm/cluster_config.yml b/profiles/slurm/cluster_config.yml
deleted file mode 100644
index 4b5d9c20..00000000
--- a/profiles/slurm/cluster_config.yml
+++ /dev/null
@@ -1,6 +0,0 @@
-__default__:
- partition: "holy-cow,holy-smokes"
- time: 9000
- nodes: 1
- mem: 1000
- output: "logs/slurm/slurm-%j.out"
diff --git a/profiles/slurm/config.yaml b/profiles/slurm/config.yaml
index 0fcf4d2e..baeee2c1 100644
--- a/profiles/slurm/config.yaml
+++ b/profiles/slurm/config.yaml
@@ -1,14 +1,424 @@
-jobs: 1000
-restart-times: 3
-max-jobs-per-second: 1
-max-status-checks-per-second: 3
-cores: 10000
-#local-cores: 1
-latency-wait: 60
+executor: slurm
use-conda: True
-keep-going: True
-rerun-incomplete: True
-jobscript: "slurm-jobscript.sh"
-cluster: "slurm-submit.py"
-cluster-status: "slurm-status.py"
+jobs: 100 # Have up to N jobs submitted at any given time
+latency-wait: 100 # Wait N seconds for output files due to latency
+retries: 0 # Retry jobs N times.
+# These resources will be applied to all rules. Can be overriden on a per-rule basis below.
+default-resources:
+ mem_mb: attempt * 2000
+ mem_mb_reduced: (attempt * 2000) * 0.9 # Mem allocated to java for GATK rules (tries to prevent OOM errors)
+ slurm_partition: ""
+ slurm_account: # Same as sbatch -A. Not all clusters use this.
+ runtime: 30 # In minutes
+
+
+# Control number of threads each rule will use.
+set-threads:
+ # Reference Genome Processing. Does NOT use more than 1 thread.
+ download_reference: 1
+ index_reference: 1
+ # Interval Generation. Does NOT use more than 1 thread.
+ format_interval_list: 1
+ create_gvcf_intervals: 1
+ create_db_intervals: 1
+ picard_intervals: 1
+ # Mappability
+ genmap: 1 # Can use more than 1 thread
+ mappability_bed: 1 # Does NOT use more than 1 thread
+ # Fastq Processing. Can use more than 1 thread.
+ get_fastq_pe: 1
+ fastp: 1
+ # Alignment. Can use more than 1 thread, except merge_bams.
+ bwa_map: 1
+ dedup: 1
+ merge_bams: 1 # Does NOT use more than 1 thread.
+ # GVCF
+ bam2gvcf: 1 # Should be run with no more than 2 threads.
+ concat_gvcfs: 1 # Does NOT use more than 1 thread.
+ create_db_mapfile: 1 # Does NOT use more than 1 thread.
+ gvcf2DB: 1 # Should be run with no more than 2 threads.
+ # VCF
+ DB2vcf: 1 # Should be run with no more than 2 threads.
+ filterVcfs: 1 # Should be run with no more than 2 threads.
+ sort_gatherVcfs: 1 # Should be run with no more than 2 threads.
+ # Callable Bed
+ compute_d4: 1 # Can use more than 1 thread
+ create_cov_bed: 1 # Does NOT use more than 1 thread.
+ merge_d4: 1 # Does NOT use more than 1 thread.
+ # Summary Stats Does NOT use more than 1 thread.
+ bam_sumstats: 1
+ collect_covstats: 1
+ collect_fastp_stats: 1
+ collect_sumstats: 1
+ # QC Module Does NOT use more than 1 thread.
+ qc_admixture: 1
+ qc_check_fai: 1
+ qc_generate_coords_file: 1
+ qc_plink: 1
+ qc_qc_plots: 1
+ qc_setup_admixture: 1
+ qc_subsample_snps: 1
+ qc_vcftools_individuals: 1
+ # MK Module Does NOT use more than 1 thread.
+ mk_degenotate: 1
+ mk_prep_genome: 1
+ mk_split_samples: 1
+ # Postprocess Module Does NOT use more than 1 thread.
+ postprocess_strict_filter: 1
+ postprocess_basic_filter: 1
+ postprocess_filter_individuals: 1
+ postprocess_subset_indels: 1
+ postprocess_subset_snps: 1
+ postprocess_update_bed: 1
+ # Trackhub Module Does NOT use more than 1 thread.
+ trackhub_bcftools_depth: 1
+ trackhub_bedgraph_to_bigwig: 1
+ trackhub_calc_pi: 1
+ trackhub_calc_snpden: 1
+ trackhub_calc_tajima: 1
+ trackhub_chrom_sizes: 1
+ trackhub_convert_to_bedgraph: 1
+ trackhub_strip_vcf: 1
+ trackhub_vcftools_freq: 1
+ trackhub_write_hub_files: 1
+ # Sentieon Tools. Can use more than 1 thread, except sentieon_bam_stats.
+ sentieon_map: 1
+ sentieon_dedup: 1
+ sentieon_haplotyper: 1
+ sentieon_combine_gvcf: 1
+ sentieon_bam_stats: 1 # Does NOT use more than 1 thread.
+
+# Control other resources used by each rule.
+# set-resources:
+# # Reference Genome Processing
+# copy_reference:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# download_reference:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# index_reference:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Interval Generation
+# format_interval_list:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# create_gvcf_intervals:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# create_db_intervals:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# picard_intervals:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Mappability
+# genmap:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# mappability_bed:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Fastq Processing
+# get_fastq_pe:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# fastp:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Alignment
+# bwa_map:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# dedup:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# merge_bams:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # GVCF
+# bam2gvcf: # HaplotypeCaller
+# mem_mb: attempt * 2000
+# mem_mb_reduced: (attempt * 2000) * 0.9 # Mem allocated to java (tries to prevent OOM errors)
+# slurm_partition:
+# runtime:
+# cpus_per_task: # Mem allocated to the snakemake job
+# concat_gvcfs:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# create_db_mapfile:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# gvcf2DB: # GenomicsDBImport
+# mem_mb: attempt * 2000
+# mem_mb_reduced: (attempt * 2000) * 0.9 # Mem allocated to java (tries to prevent OOM errors)
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # VCF
+# DB2vcf: # GenotypeGVCFs
+# mem_mb: attempt * 2000
+# mem_mb_reduced: (attempt * 2000) * 0.9 # Mem allocated to java (tries to prevent OOM errors)
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# filterVcfs:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# sort_gatherVcfs:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Callable Bed
+# compute_d4:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# create_cov_bed:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# merge_d4:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Summary Stats
+# bam_sumstats:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# collect_covstats:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# collect_fastp_stats:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# collect_sumstats:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # QC Module
+# qc_admixture:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# qc_check_fai:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# qc_generate_coords_file:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# qc_plink:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# qc_qc_plots:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# qc_setup_admixture:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# qc_subsample_snps:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# qc_vcftools_individuals:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # MK Module
+# mk_degenotate:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# mk_prep_genome:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# mk_split_samples:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Postprocess Module
+# postprocess_strict_filter:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# postprocess_basic_filter:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# postprocess_filter_individuals:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# postprocess_subset_indels:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# postprocess_subset_snps:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# postprocess_update_bed:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Trackhub Module
+# trackhub_bcftools_depth:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_bedgraph_to_bigwig:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_calc_pi:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_calc_snpden:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_calc_tajima:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_chrom_sizes:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_convert_to_bedgraph:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_strip_vcf:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_vcftools_freq:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# trackhub_write_hub_files:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+
+# # Sentieon Tools
+# sentieon_map:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# sentieon_dedup:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# sentieon_haplotyper:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# sentieon_combine_gvcf:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
+# sentieon_bam_stats:
+# mem_mb: attempt * 2000
+# slurm_partition:
+# runtime:
+# cpus_per_task:
diff --git a/profiles/slurm/slurm-jobscript.sh b/profiles/slurm/slurm-jobscript.sh
deleted file mode 100755
index 391741ef..00000000
--- a/profiles/slurm/slurm-jobscript.sh
+++ /dev/null
@@ -1,3 +0,0 @@
-#!/bin/bash
-# properties = {properties}
-{exec_job}
diff --git a/profiles/slurm/slurm-status.py b/profiles/slurm/slurm-status.py
deleted file mode 100755
index 1929b785..00000000
--- a/profiles/slurm/slurm-status.py
+++ /dev/null
@@ -1,69 +0,0 @@
-#!/usr/bin/env python3
-import re
-import subprocess as sp
-import shlex
-import sys
-import time
-import logging
-
-logger = logging.getLogger("__name__")
-
-STATUS_ATTEMPTS = 20
-
-jobid = sys.argv[1]
-
-for i in range(STATUS_ATTEMPTS):
- try:
- sacct_res = sp.check_output(shlex.split("sacct -P -b -j {} -n".format(jobid)))
- res = {
- x.split("|")[0]: x.split("|")[1]
- for x in sacct_res.decode().strip().split("\n")
- }
- break
- except sp.CalledProcessError as e:
- logger.error("sacct process error")
- logger.error(e)
- except IndexError as e:
- pass
- # Try getting job with scontrol instead in case sacct is misconfigured
- try:
- sctrl_res = sp.check_output(
- shlex.split("scontrol -o show job {}".format(jobid))
- )
- m = re.search("JobState=(\w+)", sctrl_res.decode())
- res = {jobid: m.group(1)}
- break
- except sp.CalledProcessError as e:
- logger.error("scontrol process error")
- logger.error(e)
- if i >= STATUS_ATTEMPTS - 1:
- print("failed")
- exit(0)
- else:
- time.sleep(1)
-
-status = res[jobid]
-
-if status == "BOOT_FAIL":
- print("failed")
-elif status == "OUT_OF_MEMORY":
- print("failed")
-elif status.startswith("CANCELLED"):
- print("failed")
-elif status == "COMPLETED":
- print("success")
-elif status == "DEADLINE":
- print("failed")
-elif status == "FAILED":
- print("failed")
-elif status == "NODE_FAIL":
- print("failed")
-elif status == "PREEMPTED":
- print("failed")
-elif status == "TIMEOUT":
- print("failed")
-# Unclear whether SUSPENDED should be treated as running or failed
-elif status == "SUSPENDED":
- print("failed")
-else:
- print("running")
diff --git a/profiles/slurm/slurm-submit.py b/profiles/slurm/slurm-submit.py
deleted file mode 100755
index 2614f324..00000000
--- a/profiles/slurm/slurm-submit.py
+++ /dev/null
@@ -1,59 +0,0 @@
-#!/usr/bin/env python
-"""
-Snakemake SLURM submit script.
-"""
-import warnings # use warnings.warn() rather than print() to output info in this script
-
-from snakemake.utils import read_job_properties
-
-import slurm_utils
-
-# cookiecutter arguments
-SBATCH_DEFAULTS = """ """
-CLUSTER_CONFIG = "./cluster_config.yml"
-ADVANCED_ARGUMENT_CONVERSION = {"yes": True, "no": False}["no"]
-
-RESOURCE_MAPPING = {
- "time": ("time", "runtime", "walltime"),
- "mem": ("mem", "mem_mb", "ram", "memory"),
- "mem-per-cpu": ("mem-per-cpu", "mem_per_cpu", "mem_per_thread"),
- "nodes": ("nodes", "nnodes")
-}
-
-# parse job
-jobscript = slurm_utils.parse_jobscript()
-job_properties = read_job_properties(jobscript)
-
-sbatch_options = {}
-cluster_config = slurm_utils.load_cluster_config(CLUSTER_CONFIG)
-
-# 1) sbatch default arguments
-sbatch_options.update(slurm_utils.parse_sbatch_defaults(SBATCH_DEFAULTS))
-
-# 2) cluster_config defaults
-sbatch_options.update(cluster_config["__default__"])
-
-# 3) Convert resources (no unit conversion!) and threads
-sbatch_options.update(
- slurm_utils.convert_job_properties(job_properties, RESOURCE_MAPPING)
-)
-
-# 4) cluster_config for particular rule
-sbatch_options.update(cluster_config.get(job_properties.get("rule"), {}))
-
-# 5) cluster_config options
-sbatch_options.update(job_properties.get("cluster", {}))
-
-# 6) Advanced conversion of parameters
-if ADVANCED_ARGUMENT_CONVERSION:
- sbatch_options = slurm_utils.advanced_argument_conversion(sbatch_options)
-
-#7) Format pattern in snakemake style
-sbatch_options = slurm_utils.format_values(sbatch_options, job_properties)
-
-# ensure sbatch output dirs exist
-for o in ("output", "error"):
- slurm_utils.ensure_dirs_exist(sbatch_options[o]) if o in sbatch_options else None
-
-# submit job and echo id back to Snakemake (must be the only stdout)
-print(slurm_utils.submit_job(jobscript, **sbatch_options))
diff --git a/profiles/slurm/slurm_utils.py b/profiles/slurm/slurm_utils.py
deleted file mode 100644
index a420d1c3..00000000
--- a/profiles/slurm/slurm_utils.py
+++ /dev/null
@@ -1,260 +0,0 @@
-#!/usr/bin/env python3
-import os
-import re
-import math
-import argparse
-import subprocess
-
-from snakemake import io
-from snakemake.io import Wildcards
-from snakemake.utils import SequenceFormatter, AlwaysQuotedFormatter, QuotedFormatter
-from snakemake.exceptions import WorkflowError
-
-def parse_jobscript():
- """Minimal CLI to require/only accept single positional argument."""
- p = argparse.ArgumentParser(description="SLURM snakemake submit script")
- p.add_argument("jobscript", help="Snakemake jobscript with job properties.")
- return p.parse_args().jobscript
-
-
-def parse_sbatch_defaults(parsed):
- """Unpack SBATCH_DEFAULTS."""
- d = parsed.split() if type(parsed) == str else parsed
- args = {k.strip().strip("-"): v.strip() for k, v in [a.split("=") for a in d]}
- return args
-
-
-def load_cluster_config(path):
- """Load config to dict either from absolute path or relative to profile dir."""
- if path:
- path = os.path.join(os.path.dirname(__file__), os.path.expandvars(path))
- dcc = io.load_configfile(path)
- else:
- dcc = {}
- if "__default__" not in dcc:
- dcc["__default__"] = {}
- return dcc
-
-# adapted from format function in snakemake.utils
-def format(_pattern, _quote_all=False, **kwargs):
- """Format a pattern in Snakemake style.
- This means that keywords embedded in braces are replaced by any variable
- values that are available in the current namespace.
- """
- fmt = SequenceFormatter(separator=" ")
- if _quote_all:
- fmt.element_formatter = AlwaysQuotedFormatter()
- else:
- fmt.element_formatter = QuotedFormatter()
- try:
- return fmt.format(_pattern, **kwargs)
- except KeyError as ex:
- raise NameError(
- "The name {} is unknown in this context. Please "
- "make sure that you defined that variable. "
- "Also note that braces not used for variable access "
- "have to be escaped by repeating them "
- )
-
-# adapted from Job.format_wildcards in snakemake.jobs
-def format_wildcards(string, job_properties):
- """ Format a string with variables from the job. """
-
- class Job(object):
- def __init__(self, job_properties):
- for key in job_properties:
- setattr(self, key, job_properties[key])
- job = Job(job_properties)
- job._format_params = Wildcards(fromdict=job_properties['params'])
- job._format_wildcards = Wildcards(fromdict=job_properties['wildcards'])
- _variables = dict()
- _variables.update(
- dict(
- params=job._format_params,
- wildcards=job._format_wildcards,
- rule=job.rule
- )
- )
- try:
- return format(string, **_variables)
- except NameError as ex:
- raise WorkflowError(
- "NameError with group job {}: {}".format(job.jobid, str(ex))
- )
- except IndexError as ex:
- raise WorkflowError(
- "IndexError with group job {}: {}".format(job.jobid, str(ex))
- )
-
-# adapted from ClusterExecutor.cluster_params function in snakemake.executor
-def format_values(dictionary, job_properties):
- formatted = dictionary.copy()
- for key, value in list(formatted.items()):
- if isinstance(value, str):
- try:
- formatted[key] = format_wildcards(value, job_properties)
- except NameError as e:
- msg = (
- "Failed to format cluster config "
- "entry for job {}.".format(job_properties['rule'])
- )
- raise WorkflowError(msg, e)
- return formatted
-
-def convert_job_properties(job_properties, resource_mapping={}):
- options = {}
- resources = job_properties.get("resources", {})
- for k, v in resource_mapping.items():
- options.update({k: resources[i] for i in v if i in resources})
-
- if "threads" in job_properties:
- options["cpus-per-task"] = job_properties["threads"]
- return options
-
-
-def ensure_dirs_exist(path):
- """Ensure output folder for Slurm log files exist."""
- di = os.path.dirname(path)
- if di == "":
- return
- if not os.path.exists(di):
- os.makedirs(di, exist_ok=True)
- return
-
-
-def submit_job(jobscript, **sbatch_options):
- """Submit jobscript and return jobid."""
- optsbatch_options = [f"--{k}={v}" for k, v in sbatch_options.items()]
- try:
- res = subprocess.check_output(["sbatch"] + optsbatch_options + [jobscript])
- except subprocess.CalledProcessError as e:
- raise e
- # Get jobid
- res = res.decode()
- try:
- jobid = re.search(r"Submitted batch job (\d+)", res).group(1)
- except Exception as e:
- raise e
- return jobid
-
-
-def advanced_argument_conversion(arg_dict):
- """Experimental adjustment of sbatch arguments to the given or default partition.
- """
- adjusted_args = {}
-
- partition = arg_dict.get("partition", None) or _get_default_partition()
- constraint = arg_dict.get("constraint", None)
- ncpus = int(arg_dict.get("cpus-per-task", 1))
- nodes = int(arg_dict.get("nodes", 1))
- mem = arg_dict.get("mem", None)
- # Determine partition with features. If no constraints have been set,
- # select the partition with lowest memory
- try:
- config = _get_cluster_configuration(partition)
- mem_feat = _get_features_and_memory(partition)
- MEMORY_PER_PARTITION = _get_available_memory(mem_feat, constraint)
- MEMORY_PER_CPU = MEMORY_PER_PARTITION / int(config["cpus"])
- except Exception as e:
- print(e)
- raise e
-
- # Adjust memory in the single-node case only; getting the
- # functionality right for multi-node multi-cpu jobs requires more
- # development
- if "nodes" not in arg_dict or nodes == 1:
- if mem:
- adjusted_args["mem"] = min(int(mem), MEMORY_PER_PARTITION)
- AVAILABLE_MEM = ncpus * MEMORY_PER_CPU
- if adjusted_args["mem"] > AVAILABLE_MEM:
- adjusted_args["cpus-per-task"] = int(
- math.ceil(int(mem) / MEMORY_PER_CPU)
- )
- adjusted_args["cpus-per-task"] = min(int(config["cpus"]), ncpus)
- else:
- if nodes == 1:
- # Allocate at least as many tasks as requested nodes
- adjusted_args["cpus-per-task"] = nodes
- # Update time. If requested time is larger than maximum allowed time, reset
- try:
- if "time" in arg_dict:
- adjusted_args["time"] = min(int(config["time"]), int(arg_dict["time"]))
- except Exception as e:
- print(e)
- raise e
- # update and return
- arg_dict.update(adjusted_args)
- return arg_dict
-
-
-def _get_default_partition():
- """Retrieve default partition for cluster"""
- res = subprocess.check_output(["sinfo", "-O", "partition"])
- m = re.search(r"(?PsnpArcher Track Hub Description
+ +Introduction
+To facilitate downstream data exploration and as an example of the module development components of this work, we + developed a module to generate UCSC Genome Browser track files to explore population variation data (see preprint for details).
+ +Track Descriptions
+ +Tajima’s D
+This track provides windowed estimates of Tajima’s D, a population genetic statistic that measures the departure from neutral evolution in a DNA sequence.
+ +SNP Density
+This track displays the density of single nucleotide polymorphisms (SNPs) across the genome, showing regions with high or low levels of genetic variation.
+ +Pi
+The Pi track represents the average number of nucleotide differences per site between any two sequences in a population, providing an estimate of genetic diversity.
+ +Minor Allele Frequency
+This track shows the frequency of the less common allele at a SNP locus, providing insights into the genetic variation within a population.
+ +SNP Depth
+The SNP Depth track displays the number of reads or sequencing depth at each SNP position, indicating the coverage and quality of the variant calls.
+ +Non Callable Sites
+The Non Callable Sites track highlights regions in the genome that are considered non-callable, meaning that they have low sequencing coverage or other technical limitations that make it difficult to accurately determine genetic variation in those regions.
+ + + \ No newline at end of file diff --git a/workflow/modules/trackhub/scripts/write_hub_files.py b/workflow/modules/trackhub/scripts/write_hub_files.py index 308c33a1..bf0a3ce6 100644 --- a/workflow/modules/trackhub/scripts/write_hub_files.py +++ b/workflow/modules/trackhub/scripts/write_hub_files.py @@ -1,55 +1,64 @@ from os.path import basename +import shutil +import os # https://genome.ucsc.edu/goldenPath/help/hgTracksHelp.html#UseOneFile hub_text = """hub {genome} shortLabel {genome} snpArcher Track Hub longLabel {genome} snpArcher Track Hub useOneFile on -email {email}\n""" +descriptionUrl index.html +email {email} +genome {genome}\n""" vcf_track_txt = """track VCF bigDataUrl {vcf_file} shortLabel VCF longLabel VCF visibility squish +html index.html type vcfTabix\n""" window_parent_txt = """track {track_type} compositeTrack on shortLabel {track_type} longLabel {track_type} -color 200,100,0 +color {color} altColor 0,102,204 autoScale on type bigWig allButtonPair on +html index.html visibility full\n""" window_track_txt = """track {track_name} parent {parent} on bigDataUrl {data_url} type bigWig +visibility {vis} shortLabel {label} longLabel {label}\n""" allele_freq_txt = """track MinorAlleleFrequency bigDataUrl {data_url} type bigWig -color 200,100,0 +color 88,85,120 altColor 0,102,204 autoScale on visibility full shortLabel Minor Allele Frequency +html index.html longLabel Minor Allele Frequency\n""" snp_depth_txt = """track SNPDepth bigDataUrl {data_url} type bigWig -color 200,100,0 +color 120,172,145 altColor 0,102,204 autoScale on visibility full shortLabel SNP Depth +html index.html longLabel SNP Depth\n""" coverage_track_txt = """track NonCallableSites @@ -57,9 +66,82 @@ shortLabel Non Callable Sites type bigBed longLabel Non Callable Sites -color 0,0,255 +color 0,0,0 +html index.html visibility dense\n""" +COLORS = { + "Tajima": "(70,130,180)", + "SNP-Density": "(186,85,211)", + "Pi": "(248,174,51)", +} + +html = """ + + +snpArcher Track Hub Description
+ +Introduction
+To facilitate downstream data exploration and as an example of the module development components of this work, we + developed a module to generate UCSC Genome Browser track files to explore population variation data (see preprint for details).
+ +Track Descriptions
+ +Tajima’s D
+This track provides windowed estimates of Tajima’s D, a population genetic statistic that measures the departure from neutral evolution in a DNA sequence.
+ +SNP Density
+This track displays the density of single nucleotide polymorphisms (SNPs) across the genome, showing regions with high or low levels of genetic variation.
+ +Pi
+The Pi track represents the average number of nucleotide differences per site between any two sequences in a population, providing an estimate of genetic diversity.
+ +Minor Allele Frequency
+This track shows the frequency of the less common allele at a SNP locus, providing insights into the genetic variation within a population.
+ +SNP Depth
+The SNP Depth track displays the number of reads or sequencing depth at each SNP position, indicating the coverage and quality of the variant calls.
+ +Non Callable Sites
+The Non Callable Sites track highlights regions in the genome that are considered non-callable, meaning that they have low sequencing coverage or other technical limitations that make it difficult to accurately determine genetic variation in those regions.
+ + + + + + + +""" + def human_format(num): num = float("{:.3g}".format(num)) @@ -78,29 +160,45 @@ def main(): trackhub_windows = snakemake.params["windows"] # noqa: F821 vcf_file = basename(snakemake.input["vcf"][0]) # noqa: F821 cov_file = basename(snakemake.input["callable_sites"][0]) # noqa: F821 - freq_file = basename(snakemake.input["allele_freq"][0]) # noqa: F821 + freq_file = basename(snakemake.input["allele_freq"][0]) # noqa: F821 depth_file = basename(snakemake.input["depth"][0]) # noqa: F821 genome = snakemake.params["refGenome"] # noqa: F821 trackhub_file = snakemake.output["trackhub_file"] # noqa: F821 + html_file = snakemake.output["html"] # noqa: F821 + + with open(html_file, "w") as f: + f.write(html) with open(trackhub_file, "w") as out: print(hub_text.format(genome=genome, email=email), file=out) - print(f"genome {genome}\n", file=out) print(vcf_track_txt.format(vcf_file=vcf_file), file=out) print(coverage_track_txt.format(cov_file=cov_file), file=out) print(allele_freq_txt.format(data_url=freq_file), file=out) print(snp_depth_txt.format(data_url=depth_file), file=out) for file in file_types: - print(window_parent_txt.format(track_type=file), file=out) + print( + window_parent_txt.format(track_type=file, color=COLORS[file]), file=out + ) for window in trackhub_windows: track_name = f"{file}_{human_format(window)}_bp_bins" - label = f"{file} {human_format(window)} bp bins" + label = f"{file}_{human_format(window)}_bp bins" url = f"{file}_{window}.bw" + if window == 1000: + vis = "True" + else: + vis = "False" print( - window_track_txt.format(track_name=track_name, label=label, parent=file, data_url=url), + window_track_txt.format( + track_name=track_name, + label=label, + parent=file, + data_url=url, + vis=vis, + ), file=out, ) + if __name__ == "__main__": main() diff --git a/workflow/rules/bam2vcf_gatk.smk b/workflow/rules/bam2vcf_gatk.smk index 6e361a6a..03712b3b 100644 --- a/workflow/rules/bam2vcf_gatk.smk +++ b/workflow/rules/bam2vcf_gatk.smk @@ -14,11 +14,6 @@ rule bam2gvcf: output: gvcf = "results/{refGenome}/gvcfs/{sample}.g.vcf.gz", tbi = "results/{refGenome}/gvcfs/{sample}.g.vcf.gz.tbi" - resources: - #!The -Xmx value the tool is run with should be less than the total amount of physical memory available by at least a few GB - # subtract that memory here - mem_mb = lambda wildcards, attempt: attempt * resources['bam2gvcf']['mem'], # this is the overall memory requested - reduced = lambda wildcards, attempt: attempt * (resources['bam2gvcf']['mem'] - 3000) # this is the maximum amount given to java log: "logs/{refGenome}/gatk_hc/{sample}.txt" benchmark: @@ -26,14 +21,16 @@ rule bam2gvcf: params: minPrun = config['minP'], minDang = config['minD'], + ploidy = config['ploidy'] conda: "../envs/bam2vcf.yml" shell: "gatk HaplotypeCaller " - "--java-options \"-Xmx{resources.reduced}m\" " + "--java-options \"-Xmx{resources.mem_mb_reduced}m\" " "-R {input.ref} " "-I {input.bam} " "-O {output.gvcf} " + "-ploidy {params.ploidy} " "--emit-ref-confidence GVCF --min-pruning {params.minPrun} --min-dangling-branch-length {params.minDang} &> {log}" rule create_db_mapfile: @@ -75,9 +72,6 @@ rule gvcf2DB: output: db = temp(directory("results/{refGenome}/genomics_db_import/DB")), tar = temp("results/{refGenome}/genomics_db_import/DB.tar"), - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['gvcf2DB']['mem'], # this is the overall memory requested - reduced = lambda wildcards, attempt: int(attempt * resources['gvcf2DB']['mem'] * 0.80) # this is the maximum amount given to java log: "logs/{refGenome}/gatk_db_import.txt" benchmark: @@ -90,7 +84,7 @@ rule gvcf2DB: """ export TILEDB_DISABLE_FILE_LOCKING=1 gatk GenomicsDBImport \ - --java-options '-Xmx{resources.reduced}m -Xms{resources.reduced}m' \ + --java-options '-Xmx{resources.mem_mb_reduced}m -Xms{resources.mem_mb_reduced}m' \ --genomicsdb-shared-posixfs-optimizations true \ --batch-size 25 \ --genomicsdb-workspace-path {output.db} \ @@ -116,9 +110,6 @@ rule DB2vcf: params: het = config['het_prior'], db = lambda wc, input: input.db[:-4] - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['DB2vcf']['mem'], # this is the overall memory requested - reduced = lambda wildcards, attempt: attempt * (resources['DB2vcf']['mem'] - 3000) # this is the maximum amount given to java log: "logs/{refGenome}/gatk_genotype_gvcfs.txt" benchmark: @@ -129,7 +120,7 @@ rule DB2vcf: """ tar -xf {input.db} gatk GenotypeGVCFs \ - --java-options '-Xmx{resources.reduced}m -Xms{resources.reduced}m' \ + --java-options '-Xmx{resources.mem_mb_reduced}m -Xms{resources.mem_mb_reduced}m' \ -R {input.ref} \ --heterozygosity {params.het} \ --genomicsdb-shared-posixfs-optimizations true \ @@ -151,8 +142,6 @@ rule filterVcfs: vcfidx = temp("results/{refGenome}/vcfs/filtered.vcf.gz.tbi") conda: "../envs/bam2vcf.yml" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['filterVcfs']['mem'] # this is the overall memory requested log: "logs/{refGenome}/gatk_filter.txt" benchmark: @@ -186,9 +175,6 @@ rule sort_gatherVcfs: "logs/{refGenome}/sort_gather_vcfs/{prefix}_log.txt" benchmark: "benchmarks/{refGenome}/sort_gather_vcfs/{prefix}_benchmark.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['gatherVcfs']['mem'], # this is the overall memory requested - reduced = lambda wildcards, attempt: attempt * (resources['gatherVcfs']['mem'] - 2000) # this is the maximum amount given to java shell: """ bcftools sort -Oz -o {output.vcfFinal} {input.vcf} 2>> {log} diff --git a/workflow/rules/bam2vcf_gatk_intervals.smk b/workflow/rules/bam2vcf_gatk_intervals.smk index 79e68c8c..79431d3b 100644 --- a/workflow/rules/bam2vcf_gatk_intervals.smk +++ b/workflow/rules/bam2vcf_gatk_intervals.smk @@ -14,11 +14,6 @@ rule bam2gvcf: output: gvcf = "results/{refGenome}/interval_gvcfs/{sample}/{l}.raw.g.vcf.gz", gvcf_idx = "results/{refGenome}/interval_gvcfs/{sample}/{l}.raw.g.vcf.gz.tbi" - resources: - #!The -Xmx value the tool is run with should be less than the total amount of physical memory available by at least a few GB - # subtract that memory here - mem_mb = lambda wildcards, attempt: attempt * resources['bam2gvcf']['mem'], # this is the overall memory requested - reduced = lambda wildcards, attempt: attempt * (resources['bam2gvcf']['mem'] - 3000) # this is the maximum amount given to java log: "logs/{refGenome}/gatk_hc/{sample}/{l}.txt" benchmark: @@ -26,20 +21,25 @@ rule bam2gvcf: params: minPrun = config['minP'], minDang = config['minD'], + ploidy = config['ploidy'], conda: "../envs/bam2vcf.yml" shell: - "gatk HaplotypeCaller " - "--java-options \"-Xmx{resources.reduced}m\" " - "-R {input.ref} " - "-I {input.bam} " - "-O {output.gvcf} " - "-L {input.l} " - "--emit-ref-confidence GVCF --min-pruning {params.minPrun} --min-dangling-branch-length {params.minDang} &> {log}" + """ + gatk HaplotypeCaller \ + --java-options -Xmx{resources.mem_mb_reduced}m \ + -R {input.ref} \ + -I {input.bam} \ + -O {output.gvcf} \ + -L {input.l} \ + -ploidy {params.ploidy} \ + --emit-ref-confidence GVCF --min-pruning {params.minPrun} --min-dangling-branch-length {params.minDang} &> {log} + """ rule concat_gvcfs: input: - unpack(get_interval_gvcfs) + gvcfs = get_interval_gvcfs, + tbis = get_interval_gvcfs_idx output: gvcf = "results/{refGenome}/gvcfs/{sample}.g.vcf.gz", tbi = "results/{refGenome}/gvcfs/{sample}.g.vcf.gz.tbi" @@ -48,7 +48,6 @@ rule concat_gvcfs: benchmark: "benchmarks/{refGenome}/concat_gvcfs/{sample}.txt" resources: - mem_mb = lambda wildcards, attempt: attempt * resources['gatherVcfs']['mem'], # this is the overall memory requested tmpdir = get_big_temp conda: "../envs/bcftools.yml" @@ -74,7 +73,7 @@ rule create_db_mapfile: rule gvcf2DB: """ - todo + Create GenomicsDB. """ input: unpack(get_gvcfs_db), @@ -83,15 +82,12 @@ rule gvcf2DB: output: db = temp(directory("results/{refGenome}/genomics_db_import/DB_L{l}")), tar = temp("results/{refGenome}/genomics_db_import/DB_L{l}.tar"), - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['gvcf2DB']['mem'], # this is the overall memory requested - reduced = lambda wildcards, attempt: int(attempt * resources['gvcf2DB']['mem'] * 0.80) # this is the maximum amount given to java log: "logs/{refGenome}/gatk_db_import/{l}.txt" - resources: - tmpdir = get_big_temp benchmark: "benchmarks/{refGenome}/gatk_db_import/{l}.txt" + resources: + tmpdir = get_big_temp conda: "../envs/bam2vcf.yml" shell: @@ -100,7 +96,7 @@ rule gvcf2DB: """ export TILEDB_DISABLE_FILE_LOCKING=1 gatk GenomicsDBImport \ - --java-options '-Xmx{resources.reduced}m -Xms{resources.reduced}m' \ + --java-options '-Xmx{resources.mem_mb_reduced}m -Xms{resources.mem_mb_reduced}m' \ --genomicsdb-shared-posixfs-optimizations true \ --batch-size 25 \ --genomicsdb-workspace-path {output.db} \ @@ -109,7 +105,7 @@ rule gvcf2DB: --tmp-dir {resources.tmpdir} \ --sample-name-map {input.db_mapfile} &> {log} - tar --overwrite -cf {output.tar} {output.db} + tar -cf {output.tar} {output.db} """ rule DB2vcf: @@ -120,6 +116,8 @@ rule DB2vcf: input: db = "results/{refGenome}/genomics_db_import/DB_L{l}.tar", ref = "results/{refGenome}/data/genome/{refGenome}.fna", + fai = "results/{refGenome}/data/genome/{refGenome}.fna.fai", + dictf = "results/{refGenome}/data/genome/{refGenome}.dict", output: vcf = temp("results/{refGenome}/vcfs/intervals/L{l}.vcf.gz"), vcfidx = temp("results/{refGenome}/vcfs/intervals/L{l}.vcf.gz.tbi"), @@ -127,8 +125,6 @@ rule DB2vcf: het = config['het_prior'], db = lambda wc, input: input.db[:-4] resources: - mem_mb = lambda wildcards, attempt: attempt * resources['DB2vcf']['mem'], # this is the overall memory requested - reduced = lambda wildcards, attempt: attempt * (resources['DB2vcf']['mem'] - 3000), # this is the maximum amount given to java tmpdir = get_big_temp log: "logs/{refGenome}/gatk_genotype_gvcfs/{l}.txt" @@ -140,7 +136,7 @@ rule DB2vcf: """ tar -xf {input.db} gatk GenotypeGVCFs \ - --java-options '-Xmx{resources.reduced}m -Xms{resources.reduced}m' \ + --java-options '-Xmx{resources.mem_mb_reduced}m -Xms{resources.mem_mb_reduced}m' \ -R {input.ref} \ --heterozygosity {params.het} \ --genomicsdb-shared-posixfs-optimizations true \ @@ -156,14 +152,14 @@ rule filterVcfs: input: vcf = "results/{refGenome}/vcfs/intervals/L{l}.vcf.gz", vcfidx = "results/{refGenome}/vcfs/intervals/L{l}.vcf.gz.tbi", - ref = "results/{refGenome}/data/genome/{refGenome}.fna" + ref = "results/{refGenome}/data/genome/{refGenome}.fna", + fai = "results/{refGenome}/data/genome/{refGenome}.fna.fai", + dictf = "results/{refGenome}/data/genome/{refGenome}.dict", output: vcf = temp("results/{refGenome}/vcfs/intervals/filtered_L{l}.vcf.gz"), vcfidx = temp("results/{refGenome}/vcfs/intervals/filtered_L{l}.vcf.gz.tbi") conda: "../envs/bam2vcf.yml" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['filterVcfs']['mem'] # this is the overall memory requested log: "logs/{refGenome}/gatk_filter/{l}.txt" benchmark: @@ -187,7 +183,7 @@ rule filterVcfs: rule sort_gatherVcfs: input: vcfs = get_interval_vcfs, - tbis = get_interval_vcf_tbis + tbis = get_interval_vcfs_idx output: vcfFinal = "results/{refGenome}/{prefix}_raw.vcf.gz", vcfFinalidx = "results/{refGenome}/{prefix}_raw.vcf.gz.tbi" @@ -198,7 +194,6 @@ rule sort_gatherVcfs: benchmark: "benchmarks/{refGenome}/sort_gather_vcfs/{prefix}_benchmark.txt" resources: - mem_mb = lambda wildcards, attempt: attempt * resources['gatherVcfs']['mem'], # this is the overall memory requested tmpdir = get_big_temp shell: """ diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 94d0b0f5..5f54a822 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -1,103 +1,191 @@ -import glob -import re -import os import sys +import os import tempfile import random import string import statistics +from pathlib import Path from collections import defaultdict -from urllib.request import urlopen + import pandas as pd +import snparcher_utils from yaml import safe_load -from collections import defaultdict, deque -from snakemake.exceptions import WorkflowError -samples = pd.read_table(config["samples"], sep=",", dtype=str).replace(' ', '_', regex=True) + +DEFAULT_STORAGE_PREFIX = workflow.default_remote_prefix +GS_READS_PREFIX = config["remote_reads_prefix"] +samples = snparcher_utils.parse_sample_sheet(config) + with open(config["resource_config"], "r") as f: resources = safe_load(f) def get_output(): + if config["final_prefix"] == "": - raise(WorkflowError("'final_prefix' is not set in config.")) + raise (WorkflowError("'final_prefix' is not set in config.")) out = [] - genomes = samples['refGenome'].unique().tolist() - sample_counts = samples.drop_duplicates(subset = ["BioSample", "refGenome"]).value_counts(subset=['refGenome']) #get BioSample for each refGenome + genomes = samples["refGenome"].unique().tolist() + sample_counts = samples.drop_duplicates( + subset=["BioSample", "refGenome"] + ).value_counts( + subset=["refGenome"] + ) # get BioSample for each refGenome out.extend + if config["final_prefix"] == "": + raise (WorkflowError("'final_prefix' is not set in config.")) + out = [] + + sample_counts = samples.drop_duplicates( + subset=["BioSample", "refGenome"] + ).value_counts( + subset=["refGenome"] + ) # get BioSample for each refGenome + for ref in genomes: - out.extend(expand("results/{refGenome}/{prefix}_raw.vcf.gz", refGenome=ref, prefix=config['final_prefix'])) - out.extend(expand("results/{refGenome}/summary_stats/{prefix}_bam_sumstats.txt", refGenome=ref, prefix=config['final_prefix'])) - out.extend(expand("results/{refGenome}/{prefix}_callable_sites.bed", refGenome=ref, prefix=config['final_prefix'])) + # Workaround for Snakemake issue 2762. There is problem with running nested checkpoints in snakemake8. Adding mapfile in rule all forces gvcf interval checkpoint to run. + # This is actually kind of a good thing to do since it makes dryrun more clear (shows all bam>gvcf jobs now). + out.extend(expand("results/{refGenome}/genomics_db_import/DB_mapfile.txt", refGenome=ref)) + out.extend( + expand( "results/{refGenome}/{prefix}_raw.vcf.gz",refGenome=ref, prefix=config["final_prefix"])) + out.extend( + expand( "results/{refGenome}/summary_stats/{prefix}_bam_sumstats.txt", refGenome=ref, prefix=config["final_prefix"])) + out.extend( + expand("results/{refGenome}/{prefix}_callable_sites.bed", refGenome=ref, prefix=config["final_prefix"])) if sample_counts[ref] > 2: out.append(rules.qc_all.input) if "SampleType" in samples.columns: out.append(rules.postprocess_all.input) - if all(i in samples['SampleType'].tolist() for i in ["ingroup", "outgroup"]): + if all( + i in samples["SampleType"].tolist() for i in ["ingroup", "outgroup"] + ): out.append(rules.mk_all.input) if config["generate_trackhub"]: if not config["trackhub_email"]: - raise(WorkflowError("If generating trackhub, you must provide an email in the config file.")) + raise ( + WorkflowError( + "If generating trackhub, you must provide an email in the config file." + ) + ) out.append(rules.trackhub_all.input) return out + def merge_bams_input(wc): - return expand("results/{{refGenome}}/bams/preMerge/{{sample}}/{run}.bam", run=samples.loc[samples['BioSample'] == wc.sample]['Run'].tolist()) + return expand( + "results/{{refGenome}}/bams/preMerge/{{sample}}/{run}.bam", + run=samples.loc[samples["BioSample"] == wc.sample]["Run"].tolist(), + ) -def get_ref(wildcards): - if 'refPath' in samples.columns: - _refs = samples.loc[(samples['refGenome'] == wildcards.refGenome)]['refPath'].dropna().unique().tolist() - for ref in _refs: - if workflow.default_remote_prefix == "": - if not os.path.exists(ref): - raise WorkflowError(f"Reference genome {ref} does not exist") - elif ref.rsplit(".", 1)[1] == '.gz': - raise WorkflowError(f"Reference genome {ref} must be unzipped first.") - return _refs +def setup_curlrc(): + curlrc_path = Path("~/.curlrc").expanduser() + marker = "# Added by snpArcher" + entry = f"-L {marker}\n" + if curlrc_path.exists(): + with curlrc_path.open("r+") as f: + if "-L" not in f.read(): + f.write(f"\n{entry}\n") + logger.info(f"Added -L to {curlrc_path} for pyd4") + else: - return [] - + with curlrc_path.open("a+") as f: + f.write(f"{entry}\n") + +def cleanup_curlrc(): + curlrc_path = Path("~/.curlrc").expanduser() + marker = "# Added by snpArcher" + entry = f"-L {marker}\n" + logger.info(f"Removing -L we added from {curlrc_path}...") + if curlrc_path.exists(): + with curlrc_path.open("r") as f: + lines = f.readlines() + # remove entry if its there + new_lines = [line for line in lines if line.strip() != entry.strip()] + if len(new_lines) == 0: + # our entry was only thing there, we can delete .curlrc + curlrc_path.unlink() + else: + # write back any options that were there. + with curlrc_path.open("w") as f: + f.writelines(new_lines) + +def get_ref(wildcards): + + if "refPath" in samples.columns: + _refs = ( + samples.loc[(samples["refGenome"] == wildcards.refGenome)]["refPath"] + .dropna() + .unique() + .tolist() + ) + if _refs: + return _refs + # if not user-specified refpath, force MissingInputError in copy_ref with dummyfile, which allows download_ref to run b/c of ruleorder. + logger.info(f"refPath specified in sample sheet header, but no path provided for refGenome '{wildcards.refGenome}'\n" + + f"Will try to download '{wildcards.refGenome}' from NCBI. If this is a genome accession, you can ignore this warning.") + return [] + def sentieon_combine_gvcf_cmd_line(wc): - gvcfs = sentieon_combine_gvcf_input(wc)['gvcfs'] + gvcfs = sentieon_combine_gvcf_input(wc)["gvcfs"] return " ".join(["-v " + gvcf for gvcf in gvcfs]) + def get_interval_gvcfs(wc): checkpoint_output = checkpoints.create_gvcf_intervals.get(**wc).output[0] - list_files = [os.path.basename(x) for x in glob.glob(os.path.join(checkpoint_output, "*.interval_list"))] + with checkpoint_output.open() as f: + lines = [l.strip() for l in f.readlines()] + list_files = [os.path.basename(x) for x in lines] list_numbers = [f.replace("-scattered.interval_list", "") for f in list_files] - gvcfs = expand("results/{{refGenome}}/interval_gvcfs/{{sample}}/{l}.raw.g.vcf.gz", l=list_numbers) - tbis = expand("results/{{refGenome}}/interval_gvcfs/{{sample}}/{l}.raw.g.vcf.gz.tbi", l=list_numbers) - return {"gvcfs": gvcfs, "tbis": tbis} + gvcfs = expand( + "results/{{refGenome}}/interval_gvcfs/{{sample}}/{l}.raw.g.vcf.gz", + l=list_numbers, + ) + + return gvcfs + + +def get_interval_gvcfs_idx(wc): + tbis = [f + ".tbi" for f in get_interval_gvcfs(wc)] + return tbis + def get_db_interval_count(wc): - _samples = samples.loc[(samples['refGenome'] == wc.refGenome)]['BioSample'].unique().tolist() - out = max(int((config["db_scatter_factor"]) * len(_samples) * config["num_gvcf_intervals"]), 1) + _samples = ( + samples.loc[(samples["refGenome"] == wc.refGenome)]["BioSample"] + .unique() + .tolist() + ) + out = max( + int((config["db_scatter_factor"]) * len(_samples) * config["num_gvcf_intervals"]),1) return out + def get_interval_vcfs(wc): checkpoint_output = checkpoints.create_db_intervals.get(**wc).output[0] - list_files = [os.path.basename(x) for x in glob.glob(os.path.join(checkpoint_output, "*.interval_list"))] + with checkpoint_output.open() as f: + lines = [l.strip() for l in f.readlines()] + list_files = [os.path.basename(x) for x in lines] + list_numbers = [f.replace("-scattered.interval_list", "") for f in list_files] vcfs = expand("results/{{refGenome}}/vcfs/intervals/filtered_L{l}.vcf.gz", l=list_numbers) - tbis = expand("results/{{refGenome}}/vcfs/intervals/filtered_L{l}.vcf.gz.tbi", l=list_numbers) + return vcfs -def get_interval_vcf_tbis(wc): - checkpoint_output = checkpoints.create_db_intervals.get(**wc).output[0] - list_files = [os.path.basename(x) for x in glob.glob(os.path.join(checkpoint_output, "*.interval_list"))] - list_numbers = [f.replace("-scattered.interval_list", "") for f in list_files] - vcfs = expand("results/{{refGenome}}/vcfs/intervals/filtered_L{l}.vcf.gz", l=list_numbers) - tbis = expand("results/{{refGenome}}/vcfs/intervals/filtered_L{l}.vcf.gz.tbi", l=list_numbers) + +def get_interval_vcfs_idx(wc): + tbis = [f + ".tbi" for f in get_interval_vcfs(wc)] return tbis - + + def get_gvcfs_db(wc): - _samples = samples.loc[(samples['refGenome'] == wc.refGenome)]['BioSample'].unique().tolist() + _samples = samples.loc[(samples["refGenome"] == wc.refGenome)]["BioSample"].unique().tolist() gvcfs = expand("results/{{refGenome}}/gvcfs/{sample}.g.vcf.gz", sample=_samples) tbis = expand("results/{{refGenome}}/gvcfs/{sample}.g.vcf.gz.tbi", sample=_samples) return {"gvcfs": gvcfs, "tbis": tbis} + def dedup_input(wc): - runs = samples.loc[samples['BioSample'] == wc.sample]['Run'].tolist() - + runs = samples.loc[samples["BioSample"] == wc.sample]["Run"].tolist() + if len(runs) == 1: bam = expand("results/{{refGenome}}/bams/preMerge/{{sample}}/{run}.bam", run=runs) bai = expand("results/{{refGenome}}/bams/preMerge/{{sample}}/{run}.bam.bai", run=runs) @@ -106,59 +194,72 @@ def dedup_input(wc): bai = "results/{refGenome}/bams/postMerge/{sample}.bam.bai" return {"bam": bam, "bai": bai} + def sentieon_combine_gvcf_input(wc): - _samples = samples['BioSample'].unique().tolist() + _samples = samples["BioSample"].unique().tolist() gvcfs = expand("results/{{refGenome}}/gvcfs/{sample}.g.vcf.gz", sample=_samples) tbis = expand("results/{{refGenome}}/gvcfs/{sample}.g.vcf.gz.tbi", sample=_samples) return {"gvcfs": gvcfs, "tbis": tbis} + def get_reads(wc): """Returns local read files if present. Defaults to SRR if no local reads in sample sheet.""" - if config['remote_reads']: + if config["remote_reads"]: return get_remote_reads(wc) - else: - row = samples.loc[samples['Run'] == wc.run] - if 'fq1' in samples.columns and 'fq2' in samples.columns: + row = samples.loc[samples["Run"] == wc.run] + r1 = f"results/data/fastq/{wc.refGenome}/{wc.sample}/{wc.run}_1.fastq.gz" + r2 = f"results/data/fastq/{wc.refGenome}/{wc.sample}/{wc.run}_2.fastq.gz" + if "fq1" in samples.columns and "fq2" in samples.columns: + if row["fq1"].notnull().any() and row["fq2"].notnull().any(): + r1 = row.fq1.item() + r2 = row.fq2.item() + if os.path.exists(row.fq1.item()) and os.path.exists(row.fq2.item()): - r1 = row.fq1.item() - r2 = row.fq2.item() return {"r1": r1, "r2": r2} else: - raise WorkflowError(f"fq1 and fq2 specified for {wc.sample}, but files were not found.") + raise WorkflowError( + f"fq1 and fq2 specified for {wc.sample}, but files were not found." + ) else: - r1 = f"results/data/fastq/{wc.refGenome}/{wc.sample}/{wc.run}_1.fastq.gz", - r2 = f"results/data/fastq/{wc.refGenome}/{wc.sample}/{wc.run}_2.fastq.gz" + # this allows mixed srr and user-specified paths for reads return {"r1": r1, "r2": r2} + else: + return {"r1": r1, "r2": r2} def get_remote_reads(wildcards): """Use this for reads on a different remote bucket than the default.""" # print(wildcards) - row = samples.loc[samples['Run'] == wildcards.run] + row = samples.loc[samples["Run"] == wildcards.run] r1 = GS.remote(os.path.join(GS_READS_PREFIX, row.fq1.item())) r2 = GS.remote(os.path.join(GS_READS_PREFIX, row.fq2.item())) return {"r1": r1, "r2": r2} def collect_fastp_stats_input(wc): - return expand("results/{{refGenome}}/summary_stats/{{sample}}/{run}.fastp.out", run=samples.loc[samples['BioSample'] == wc.sample]['Run'].tolist()) + return expand( + "results/{{refGenome}}/summary_stats/{{sample}}/{run}.fastp.out", + run=samples.loc[samples["BioSample"] == wc.sample]["Run"].tolist(), + ) + def get_read_group(wc): """Denote sample name and library_id in read group.""" - return r"'@RG\tID:{lib}\tSM:{sample}\tPL:ILLUMINA'".format( - sample=wc.sample, - lib=wc.run + libname = samples.loc[samples["Run"] == wc.run]["LibraryName"].tolist()[0] + return r"'@RG\tID:{lib}\tSM:{sample}\tLB:{lib}\tPL:ILLUMINA'".format( + sample=wc.sample, lib=libname ) + def get_input_sumstats(wildcards): - _samples = samples.loc[(samples['refGenome'] == wildcards.refGenome)]['BioSample'].unique().tolist() + _samples = samples.loc[(samples["refGenome"] == wildcards.refGenome)]["BioSample"].unique().tolist() aln = expand("results/{{refGenome}}/summary_stats/{sample}_AlnSumMets.txt", sample=_samples) cov = expand("results/{{refGenome}}/summary_stats/{sample}_coverage.txt", sample=_samples) fastp = expand("results/{{refGenome}}/summary_stats/{sample}_fastp.out", sample=_samples) - insert = expand("results/{{refGenome}}/summary_stats/{sample}_insert_metrics.txt", sample=_samples) + insert = expand("results/{{refGenome}}/summary_stats/{sample}_insert_metrics.txt",sample=_samples) qd = expand("results/{{refGenome}}/summary_stats/{sample}_qd_metrics.txt", sample=_samples) mq = expand("results/{{refGenome}}/summary_stats/{sample}_mq_metrics.txt", sample=_samples) gc = expand("results/{{refGenome}}/summary_stats/{sample}_gc_metrics.txt", sample=_samples) gc_summary = expand("results/{{refGenome}}/summary_stats/{sample}_gc_summary.txt", sample=_samples) - if config['sentieon']: + if config["sentieon"]: out = { "alnSumMetsFiles": aln, "fastpFiles": fastp, @@ -167,7 +268,7 @@ def get_input_sumstats(wildcards): "qc_files": qd, "mq_files": mq, "gc_files": gc, - "gc_summary": gc_summary + "gc_summary": gc_summary, } return out else: @@ -178,37 +279,68 @@ def get_input_sumstats(wildcards): } return out + def get_input_for_mapfile(wildcards): - sample_names = samples.loc[(samples['refGenome'] == wildcards.refGenome)]['BioSample'].unique().tolist() + sample_names = samples.loc[(samples["refGenome"] == wildcards.refGenome)]["BioSample"].unique().tolist() return expand("results/{{refGenome}}/gvcfs/{sample}.g.vcf.gz", sample=sample_names) + def get_input_for_coverage(wildcards): # Gets the correct sample given the organism and reference genome for the bedgraph merge step - _samples = samples.loc[(samples['refGenome'] == wildcards.refGenome)]['BioSample'].unique().tolist() + _samples = samples.loc[(samples["refGenome"] == wildcards.refGenome)]["BioSample"].unique().tolist() + d4files = expand("results/{{refGenome}}/callable_sites/{sample}.per-base.d4", sample=_samples) - return {'d4files': d4files} + return {"d4files": d4files} + def get_input_covstats(wildcards): # Gets the correct sample given the organism and reference genome for the bedgraph merge step - _samples = samples.loc[(samples['refGenome'] == wildcards.refGenome)]['BioSample'].unique().tolist() - covstats = expand("results/{{refGenome}}/callable_sites/{sample}.mosdepth.summary.txt", sample=_samples) - return {'covStatFiles': covstats} + _samples = samples.loc[(samples["refGenome"] == wildcards.refGenome)]["BioSample"].unique().tolist() + + covstats = expand("results/{{refGenome}}/callable_sites/{sample}.mosdepth.summary.txt",sample=_samples) + return {"covStatFiles": covstats} + def get_bedgraphs(wildcards): """Snakemake seems to struggle with unpack() and default_remote_prefix. So we have to do this one by one.""" - _samples = samples.loc[(samples['Organism'] == wildcards.Organism) & (samples['refGenome'] == wildcards.refGenome)]['BioSample'].unique().tolist() - bedgraphFiles = expand(config['output'] + "{{Organism}}/{{refGenome}}/" + config['bamDir'] + "preMerge/{sample}" + ".sorted.bg", sample=_samples) + _samples = ( + samples.loc[ + (samples["Organism"] == wildcards.Organism) + & (samples["refGenome"] == wildcards.refGenome) + ]["BioSample"] + .unique() + .tolist() + ) + bedgraphFiles = expand( + config["output"] + + "{{Organism}}/{{refGenome}}/" + + config["bamDir"] + + "preMerge/{sample}" + + ".sorted.bg", + sample=_samples, + ) return bedgraphFiles + def get_big_temp(wildcards): - """Sets a temp dir for rules that need more temp space that is typical on some cluster environments. Defaults to system temp dir.""" - if config['bigtmp']: - if config['bigtmp'].endswith("/"): - return config['bigtmp'] + "".join(random.choices(string.ascii_uppercase, k=12)) + "/" - else: - return config['bigtmp'] + "/" + "".join(random.choices(string.ascii_uppercase, k=12)) + "/" - else: - return tempfile.gettempdir() + """Sets a temp dir for rules that need more temp space that is typical on some cluster environments. Defaults to system temp dir.""" + if config["bigtmp"]: + if config["bigtmp"].endswith("/"): + return ( + config["bigtmp"] + + "".join(random.choices(string.ascii_uppercase, k=12)) + + "/" + ) + else: + return ( + config["bigtmp"] + + "/" + + "".join(random.choices(string.ascii_uppercase, k=12)) + + "/" + ) + else: + return tempfile.gettempdir() + def collectCovStats(covSumFiles): covStats = {} @@ -228,19 +360,19 @@ def collectCovStats(covSumFiles): sampleCov[chrom].append(cov) else: sampleCov[chrom] = [cov] - + for chr in sampleCov: mean_cov = statistics.mean(sampleCov[chr]) - try: + try: std_cov = statistics.stdev(sampleCov[chr]) except: std_cov = "NA" - covStats[chr] = {"mean" : mean_cov, "stdev" : std_cov} - - return(covStats) - -def collectFastpOutput(fastpFiles): + covStats[chr] = {"mean": mean_cov, "stdev": std_cov} + + return covStats + +def collectFastpOutput(fastpFiles): FractionReadsPassFilter = defaultdict(float) NumReadsPassFilter = defaultdict(int) @@ -252,13 +384,11 @@ def collectFastpOutput(fastpFiles): f = open(fn, "r") for line in f: if "before filtering" in line: - line = next(f) line = line.split() unfiltered += int(line[2]) if "Filtering result" in line: - line = next(f) line = line.split() pass_filter += int(line[3]) @@ -293,7 +423,6 @@ def collectAlnSumMets(alnSumMetsFiles): def collectCoverageMetrics(coverageFiles): - SeqDepths = defaultdict(float) CoveredBases = defaultdict(float) @@ -351,7 +480,6 @@ def printBamSumStats( med_insert_sizes=None, med_abs_insert_std=None, ): - samples = depths.keys() if med_insert_sizes is None: with open(out_file, "w") as f: @@ -403,4 +531,3 @@ def printBamSumStats( med_abs_insert_std[samp], file=f, ) - diff --git a/workflow/rules/cov_filter.smk b/workflow/rules/cov_filter.smk index cc5c9f0e..1f841549 100644 --- a/workflow/rules/cov_filter.smk +++ b/workflow/rules/cov_filter.smk @@ -12,12 +12,8 @@ rule compute_d4: "logs/{refGenome}/compute_d4/{sample}.txt" benchmark: "benchmarks/{refGenome}/compute_d4/{sample}.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['compute_d4']['mem'] - threads: - resources['compute_d4']['threads'] params: - prefix = os.path.join(workflow.default_remote_prefix, "results/{refGenome}/callable_sites/{sample}") + prefix = os.path.join(DEFAULT_STORAGE_PREFIX, "results/{refGenome}/callable_sites/{sample}") shell: "mosdepth --d4 -t {threads} {params.prefix} {input.bam} &> {log}" @@ -32,8 +28,6 @@ rule merge_d4: "logs/{refGenome}/merge_d4/log.txt" benchmark: "benchmarks/{refGenome}/merge_d4/benchmark.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['merge_d4']['mem'] shell: "d4tools merge {input.d4files} {output} &> {log}" @@ -78,8 +72,6 @@ rule callable_bed: "../envs/cov_filter.yml" benchmark: "benchmarks/{refGenome}/callable_bed/{prefix}_benchmark.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['callable_bed']['mem'] params: merge = config['cov_merge'] shell: diff --git a/workflow/rules/fastq.smk b/workflow/rules/fastq.smk index bfad20ba..c378e13f 100644 --- a/workflow/rules/fastq.smk +++ b/workflow/rules/fastq.smk @@ -3,15 +3,12 @@ rule get_fastq_pe: temp("results/data/fastq/{refGenome}/{sample}/{run}_1.fastq.gz"), temp("results/data/fastq/{refGenome}/{sample}/{run}_2.fastq.gz") params: - outdir = os.path.join(workflow.default_remote_prefix, "results/data/fastq/{refGenome}/{sample}/") + outdir = os.path.join(DEFAULT_STORAGE_PREFIX, "results/data/fastq/{refGenome}/{sample}/") conda: "../envs/fastq2bam.yml" - threads: - resources['get_fastq_pe']['threads'] benchmark: "benchmarks/{refGenome}/getfastq/{sample}_{run}.txt" resources: - mem_mb = lambda wildcards, attempt: attempt * resources['get_fastq_pe']['mem'], tmpdir = get_big_temp shell: """ @@ -40,10 +37,6 @@ rule fastp: summ = "results/{refGenome}/summary_stats/{sample}/{run}.fastp.out" conda: "../envs/fastq2bam.yml" - threads: - resources['fastp']['threads'] - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['fastp']['mem'], log: "logs/{refGenome}/fastp/{sample}/{run}.txt" benchmark: diff --git a/workflow/rules/fastq2bam.smk b/workflow/rules/fastq2bam.smk index 71d35fc5..4bd70e48 100644 --- a/workflow/rules/fastq2bam.smk +++ b/workflow/rules/fastq2bam.smk @@ -11,10 +11,6 @@ rule bwa_map: rg = get_read_group conda: "../envs/fastq2bam.yml" - threads: - resources['bwa_map']['threads'] - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['bwa_map']['mem'] log: "logs/{refGenome}/bwa_mem/{sample}/{run}.txt" benchmark: @@ -34,8 +30,6 @@ rule merge_bams: "logs/{refGenome}/merge_bams/{sample}.txt" benchmark: "benchmarks/{refGenome}/merge_bams/{sample}.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['merge_bams']['mem'] shell: "samtools merge {output.bam} {input} && samtools index {output.bam} > {log}" @@ -47,9 +41,6 @@ rule dedup: dedupBai = "results/{refGenome}/bams/{sample}_final.bam.bai", conda: "../envs/sambamba.yml" - resources: - threads = resources['dedup']['threads'], - mem_mb = lambda wildcards, attempt: attempt * resources['dedup']['mem'] log: "logs/{refGenome}/sambamba_dedup/{sample}.txt" benchmark: diff --git a/workflow/rules/intervals.smk b/workflow/rules/intervals.smk index bda9ccd0..23ec658d 100644 --- a/workflow/rules/intervals.smk +++ b/workflow/rules/intervals.smk @@ -13,8 +13,6 @@ rule picard_intervals: "logs/{refGenome}/picard_intervals/log.txt" benchmark: "benchmarks/{refGenome}/picard_intervals/benchmark.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['process_ref']['mem'] shell: "picard ScatterIntervalsByNs REFERENCE={input.ref} OUTPUT={output.intervals} MAX_TO_MERGE={params.minNmer} OUTPUT_TYPE=ACGT &> {log}\n" @@ -36,8 +34,11 @@ rule format_interval_list: checkpoint create_db_intervals: input: ref = "results/{refGenome}/data/genome/{refGenome}.fna", + fai = "results/{refGenome}/data/genome/{refGenome}.fna.fai", + dictf = "results/{refGenome}/data/genome/{refGenome}.dict", intervals = "results/{refGenome}/intervals/master_interval_list.list" output: + fof = "results/{refGenome}/intervals/db_intervals/intervals.txt", out_dir = directory("results/{refGenome}/intervals/db_intervals"), params: max_intervals = get_db_interval_count @@ -50,16 +51,20 @@ checkpoint create_db_intervals: shell: """ gatk SplitIntervals -L {input.intervals} \ - -O {output} -R {input.ref} -scatter {params} \ + -O {output.out_dir} -R {input.ref} -scatter {params} \ -mode INTERVAL_SUBDIVISION \ --interval-merging-rule OVERLAPPING_ONLY &> {log} + ls -l {output.out_dir}/*scattered.interval_list > {output.fof} """ checkpoint create_gvcf_intervals: input: ref = "results/{refGenome}/data/genome/{refGenome}.fna", + fai = "results/{refGenome}/data/genome/{refGenome}.fna.fai", + dictf = "results/{refGenome}/data/genome/{refGenome}.dict", intervals = "results/{refGenome}/intervals/master_interval_list.list" output: + fof = "results/{refGenome}/intervals/gvcf_intervals/intervals.txt", out_dir = directory("results/{refGenome}/intervals/gvcf_intervals"), params: max_intervals = config["num_gvcf_intervals"] @@ -72,7 +77,8 @@ checkpoint create_gvcf_intervals: shell: """ gatk SplitIntervals -L {input.intervals} \ - -O {output} -R {input.ref} -scatter {params} \ + -O {output.out_dir} -R {input.ref} -scatter {params} \ -mode BALANCING_WITHOUT_INTERVAL_SUBDIVISION \ --interval-merging-rule OVERLAPPING_ONLY &> {log} + ls -l {output.out_dir}/*scattered.interval_list > {output.fof} """ \ No newline at end of file diff --git a/workflow/rules/mappability.smk b/workflow/rules/mappability.smk index b4bbb943..f5bc1754 100644 --- a/workflow/rules/mappability.smk +++ b/workflow/rules/mappability.smk @@ -5,8 +5,8 @@ rule genmap: bg = temp("results/{refGenome}/genmap/{refGenome}.genmap.bedgraph"), sorted_bg = "results/{refGenome}/genmap/sorted_mappability.bg" params: - indir = os.path.join(workflow.default_remote_prefix, "results/{refGenome}/genmap_index"), - outdir = os.path.join(workflow.default_remote_prefix, "results/{refGenome}/genmap"), + indir = os.path.join(DEFAULT_STORAGE_PREFIX, "results/{refGenome}/genmap_index"), + outdir = os.path.join(DEFAULT_STORAGE_PREFIX, "results/{refGenome}/genmap"), kmer = config['mappability_k'] log: "logs/{refGenome}/genmap/log.txt" @@ -14,10 +14,6 @@ rule genmap: "benchmarks/{refGenome}/genmap/benchmark.txt" conda: "../envs/mappability.yml" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['genmap']['mem'] - threads: - resources['genmap']['threads'] shell: # snakemake creates the output directory before the shell command, but genmap doesnt like this. so we remove the directory first. """ @@ -36,8 +32,6 @@ rule mappability_bed: "../envs/mappability.yml" benchmark: "benchmarks/{refGenome}/mapbed/{prefix}_benchmark.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['callable_bed']['mem'] params: merge = config['mappability_merge'], mappability = config['mappability_min'] diff --git a/workflow/rules/reference.smk b/workflow/rules/reference.smk index 0ad6ba44..d1605c98 100644 --- a/workflow/rules/reference.smk +++ b/workflow/rules/reference.smk @@ -1,4 +1,20 @@ ruleorder: download_reference > index_reference +# localrules: copy_reference, download_reference + +# This does not work with SLURM as of 4/3/24. See here for more info:https://github.com/snakemake/snakemake-executor-plugin-slurm/issues/60 +# rule copy_reference: +# """Copies user-specified reference genome path to results dir to maintain refGenome wildcard""" +# input: +# ref = get_ref +# output: +# ref = "results/{refGenome}/data/genome/{refGenome}.fna" +# log: +# "logs/{refGenome}/copy_ref/log.txt" +# shell: +# #probably don't need to unzip but might as well. +# """ +# gunzip -c {input.ref} 2> {log} > {output.ref} || cp {input.ref} {output.ref} &> {log} +# """ rule download_reference: input: @@ -23,9 +39,10 @@ rule download_reference: && 7z x {params.dataset} -aoa -o{params.outdir} \ && cat {params.outdir}/ncbi_dataset/data/{wildcards.refGenome}/*.fna > {output.ref} else - cp {input.ref} {output.ref} + gunzip -c {input.ref} 2> {log} > {output.ref} || cp {input.ref} {output.ref} &> {log} fi """ + rule index_reference: input: ref = "results/{refGenome}/data/genome/{refGenome}.fna" @@ -35,8 +52,6 @@ rule index_reference: dictf = "results/{refGenome}/data/genome/{refGenome}.dict" conda: "../envs/fastq2bam.yml" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['index_ref']['mem'] log: "logs/{refGenome}/index_ref/log.txt" benchmark: diff --git a/workflow/rules/sentieon.smk b/workflow/rules/sentieon.smk index d47cc6da..fb0f9e37 100644 --- a/workflow/rules/sentieon.smk +++ b/workflow/rules/sentieon.smk @@ -12,14 +12,10 @@ rule sentieon_map: lic = config['sentieon_lic'] conda: "../envs/sentieon.yml" - threads: resources['sentieon_map']['threads'] log: "logs/{refGenome}/sentieon_map/{sample}/{run}.txt" benchmark: "benchmarks/{refGenome}/sentieon_map/{sample}/{run}.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['sentieon_map']['mem'], - machine_type = resources['sentieon_map']['machine_type'] shell: """ export MALLOC_CONF=lg_dirty_mult:-1 @@ -39,8 +35,6 @@ rule merge_bams: "logs/{refGenome}/merge_bams/{sample}.txt" benchmark: "benchmarks/{refGenome}/merge_bams/{sample}.txt" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['merge_bams']['mem'] shell: "samtools merge {output.bam} {input} && samtools index {output.bam}" @@ -60,11 +54,6 @@ rule sentieon_dedup: "logs/{refGenome}/sentieon_dedup/{sample}.txt" benchmark: "benchmarks/{refGenome}/sentieon_dedup/{sample}.txt" - threads: - resources['sentieon_dedup']['threads'] - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['sentieon_dedup']['mem'], - machine_type = resources['sentieon_dedup']['machine_type'] shell: """ export SENTIEON_LICENSE={params.lic} @@ -80,14 +69,11 @@ rule sentieon_haplotyper: bam = "results/{refGenome}/bams/{sample}_final.bam", bai = "results/{refGenome}/bams/{sample}_final.bam.bai" params: - lic = config['sentieon_lic'] + lic = config['sentieon_lic'], + ploidy = config['ploidy'] output: gvcf = "results/{refGenome}/gvcfs/{sample}.g.vcf.gz", gvcf_idx = "results/{refGenome}/gvcfs/{sample}.g.vcf.gz.tbi", - threads: resources['sentieon_haplotyper']['threads'] - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['sentieon_haplotyper']['mem'], - machine_type = resources['sentieon_haplotyper']['machine_type'] conda: "../envs/sentieon.yml" log: @@ -97,7 +83,7 @@ rule sentieon_haplotyper: shell: """ export SENTIEON_LICENSE={params.lic} - sentieon driver -r {input.ref} -t {threads} -i {input.bam} --algo Haplotyper --genotype_model multinomial --emit_mode gvcf --emit_conf 30 --call_conf 30 {output.gvcf} 2> {log} + sentieon driver -r {input.ref} -t {threads} -i {input.bam} --algo Haplotyper --genotype_model multinomial --emit_mode gvcf --emit_conf 30 --call_conf 30 {output.gvcf} --ploidy {params.ploidy} 2> {log} """ rule sentieon_combine_gvcf: @@ -113,10 +99,6 @@ rule sentieon_combine_gvcf: glist = lambda wc, input: " ".join(["-v " + gvcf for gvcf in input['gvcfs']]), lic = config['sentieon_lic'] threads: resources['sentieon_combine_gvcf']['threads'] - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['sentieon_combine_gvcf']['mem'], - machine_type = resources['sentieon_combine_gvcf']['machine_type'], - disk_mb = resources['sentieon_combine_gvcf']['disk_mb'] conda: "../envs/sentieon.yml" log: @@ -144,8 +126,6 @@ rule filter_vcf: tbi = "results/{refGenome}/{prefix}_raw.vcf.gz.tbi" conda: "../envs/bam2vcf.yml" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['filterVcfs']['mem'], # this is the overall memory requested log: "logs/{refGenome}/sentieon_combine_gvcf/{prefix}_log.txt" benchmark: diff --git a/workflow/rules/sumstats.smk b/workflow/rules/sumstats.smk index 246a0eac..2124053d 100644 --- a/workflow/rules/sumstats.smk +++ b/workflow/rules/sumstats.smk @@ -8,8 +8,6 @@ rule bam_sumstats: alnSum = "results/{refGenome}/summary_stats/{sample}_AlnSumMets.txt", conda: "../envs/fastq2bam.yml" - resources: - mem_mb = lambda wildcards, attempt: attempt * resources['bam_sumstats']['mem'] shell: """ samtools coverage --output {output.cov} {input.bam} diff --git a/workflow/snparcher_utils/__init__.py b/workflow/snparcher_utils/__init__.py new file mode 100644 index 00000000..79c3923a --- /dev/null +++ b/workflow/snparcher_utils/__init__.py @@ -0,0 +1,47 @@ +import pandas as pd +from snakemake.exceptions import WorkflowError +from pathlib import Path + +def parse_sample_sheet(config: dict) -> pd.DataFrame: + samples = pd.read_table(config["samples"], sep=",", dtype=str).replace(' ', '_', regex=True) + config_genomes = get_config_genomes(config, samples) + refGenome = 'refGenome' in samples.columns and samples['refGenome'].notna().any() + refPath = 'refPath' in samples.columns and samples['refPath'].notna().any() + if not any([config_genomes, refGenome, refPath]): + raise WorkflowError("No 'refGenome' or 'refPath' found in config or sample sheet.") + if config_genomes is not None: + config_refGenome, config_refPath = config_genomes + samples["refGenome"] = config_refGenome + samples["refPath"] = config_refPath + if 'refPath' in samples.columns and samples['refPath'].notna().any(): + check_ref_paths(samples) + return samples + +def get_config_genomes(config: dict, samples: pd.DataFrame): + refGenome = config.get("refGenome", False) + refPath = config.get("refPath", False) + + if refGenome and refPath: + if 'refGenome' in samples.columns and samples['refGenome'].notna().any(): + raise WorkflowError("'refGenome' is set in sample sheet AND in config. These are mutually exclusive.") + return refGenome, refPath + elif refGenome and not refPath: + raise WorkflowError("'refGenome' is set in config, but 'refPath' is not. Both are required to use these settings.") + elif refPath and not refGenome: + raise WorkflowError("'refPath' is set in config, but 'refGenome' is not. Both are required to use these settings.") + return None + +def check_ref_paths(samples: pd.DataFrame) -> None: + """ + Checks reference paths to make sure they exist, otherwise we might try to download them based on refGenome. + Also make sure only one refPath per refGenome. + """ + for refname in samples["refGenome"].dropna().tolist(): + refs = samples[samples["refGenome"] == refname]["refPath"].dropna().unique().tolist() + if len(refs) > 1: + raise WorkflowError(f"refGenome '{refname}' has more than one unique 'refPath' specified: {refs}") + for ref in refs: + if not Path(ref).exists: + raise WorkflowError(f"refPath: '{ref}' was specified in sample sheet, but could not be found.") + + diff --git a/workflow/write_samples.py b/workflow/snparcher_utils/write_samples.py similarity index 86% rename from workflow/write_samples.py rename to workflow/snparcher_utils/write_samples.py index 97f183fa..13a37ede 100644 --- a/workflow/write_samples.py +++ b/workflow/snparcher_utils/write_samples.py @@ -39,17 +39,17 @@ def copy_reference(ref: Path) -> str: shutil.copyfile(ref, Path('data', 'genome', ref_name + ".fna")) return ref_name -def write_sample_sheet(sample_dict: dict, ref_name: str, ref_path: str, organism: str, ncbi_ref: bool) -> None: +def write_sample_sheet(sample_dict: dict, ref_name: str, ref_path: str, ncbi_ref: bool) -> None: """Writes the sample sheet""" with open(Path("../config", "samples.csv"), "w") as out: if (ncbi_ref): - out.write("BioSample,LibraryName,refGenome,Run,Organism,BioProject,fq1,fq2\n") + out.write("BioSample,LibraryName,refGenome,Run,BioProject,fq1,fq2\n") for i, (k, v) in enumerate(sample_dict.items()): - out.write(f"{k},lib_{k},{ref_name},{i},{organism},NaN,{v[0]},{v[1]}\n") + out.write(f"{k},lib_{k},{ref_name},{i},NaN,{v[0]},{v[1]}\n") else: - out.write("BioSample,LibraryName,refGenome,refPath,Run,Organism,BioProject,fq1,fq2\n") + out.write("BioSample,LibraryName,refGenome,refPath,Run,BioProject,fq1,fq2\n") for i, (k, v) in enumerate(sample_dict.items()): - out.write(f"{k},lib_{k},{ref_name},{ref_path}{i},{organism},NaN,{v[0]},{v[1]}\n") + out.write(f"{k},lib_{k},{ref_name},{ref_path},{i},NaN,{v[0]},{v[1]}\n") def main() -> None: @@ -58,7 +58,7 @@ def main() -> None: parser.add_argument('-s', '--sample_list', dest='samp', required=True, help="Specify path to sample list") parser.add_argument('-f', '--fastq_dir', dest='fastq', required=True, help="Specify path to fastq dir") parser.add_argument('-c', '--copy', dest='copyref', required=False, default=False, help="Copy reference genome to data/genome dir and unzip.") - parser.add_argument('-o', '--org', dest='org', required=True, help="Specify organism name") + group = parser.add_mutually_exclusive_group(required=True) group.add_argument('-r', '--ref', dest='ref', help="Specify path to reference genome. Mutually exclusive with -a/--acc.") group.add_argument('-a', '--acc', dest='acc', help="Specify reference genome accession. Mutually exclusive with -r/--ref") @@ -66,7 +66,7 @@ def main() -> None: sample_list = args.samp fastq_dir = Path(args.fastq) - organism = args.org + with open(sample_list, "r") as f: samples = read_sample_list(f) @@ -87,10 +87,12 @@ def main() -> None: ref_name = args.acc ref_path = "" - write_sample_sheet(sample_dict, ref_name, ref_path, organism, ncbi_ref) + + write_sample_sheet(sample_dict, ref_name, ref_path, ncbi_ref) + if cant_find: print("Couldnt' find fastqs for these files:") for name in cant_find: print(name) - -main() +if __name__ == "__main__": + main()