standard pipeline working partial cases

mmeds/config.py

@@ -27,11 +27,12 @@
 # this flag is just to disable that behavior as necessary
 DB_INSTALLED = True
 
+ROOT = Path(mmeds.__file__).parent.resolve()
+HTML_DIR = Path(html.__file__).parent.resolve()
+STORAGE_DIR = Path(resources.__file__).parent.resolve()
+SNAKEMAKE_DIR = Path(snakemake.__file__).parent.resolve()
+
 if TESTING:
-    ROOT = Path(mmeds.__file__).parent.resolve()
-    HTML_DIR = Path(html.__file__).parent.resolve()
-    STORAGE_DIR = Path(resources.__file__).parent.resolve()
-    SNAKEMAKE_DIR = Path(snakemake.__file__).parent.resolve()
     # If apache is trying to access apache's home directory for `mmeds_server_data` rather than your user account's home
     # you may need to hardcode the path here
     # TODO: Add mmeds configuration options for DATABASE_DIR
@@ -57,12 +58,9 @@
         ROOT = Path('/sc/arion/projects/clemej05a/matt/mmeds-meta/')
     else:
         DATABASE_DIR = Path('/sc/arion/projects/MMEDS/mmeds_server_data')
-        ROOT = Path('/hpc/users/mmedsadmin/www/mmeds-meta/')
+        # ROOT = Path('/hpc/users/mmedsadmin/www/mmeds-meta/')
 
-    HTML_DIR = ROOT / 'mmeds/html'
     CSS_DIR = ROOT / 'mmeds/CSS'
-    STORAGE_DIR = ROOT / 'mmeds/resources'
-    SNAKEMAKE_DIR = ROOT / 'mmeds/snakemake'
     TOOLS_DIR = ROOT / 'mmeds/tools'
     WWW_ROOT = "https://mmedsadmin.dmz.hpc.mssm.edu/"
     SERVER_ROOT = WWW_ROOT + "mmeds_app/"

mmeds/snakemake/rules/common.smk

@@ -14,7 +14,7 @@ def lefse_splits(wildcards):
         subclasses = []
         if "subclasses" in config and config["subclasses"]:
             subclasses = deepcopy(config["subclasses"])
-        
+
         if lefse_class not in subclasses:
             subclasses += [lefse_class]
 
@@ -23,7 +23,7 @@ def lefse_splits(wildcards):
                              feature_table=config["tables"], lefse_class=lefse_class, subclass=subclasses)
             continue
 
-            
+
         splits += expand("results/lefse_plot_strict.{feature_table}.{lefse_class}.{subclass}.pdf",
                          feature_table=config["tables"], lefse_class=lefse_class, subclass=subclasses)
         for i in range(len(categories)-1):
@@ -32,6 +32,11 @@ def lefse_splits(wildcards):
                                      feature_table=config["tables"], lefse_class=lefse_class, cat1=categories[i], cat2=categories[j], subclass=subclasses)
     return splits
 
+def demux_single_option(wildcards):
+    components = wildcards.sequencing_run.split("_")
+    if "MSQ" in components and int(components[-1]) > 90:
+        return "--p-no-golay-error-correction"
+    return "--p-rev-comp-mapping-barcodes"
 
 def get_tool_dir():
     return TOOLS_DIR

mmeds/snakemake/rules/demux_denoise.smk

@@ -1,19 +1,21 @@
 rule demux_single_barcodes:
     input:
-        seqs = "section_{sequencing_run}/qiime_artifact.qza",
+        seqs = "section_{sequencing_run}/qiime_import_artifact.qza",
         barcodes = "section_{sequencing_run}/qiime_mapping_file_{sequencing_run}.tsv"
     output:
         error_correction = "section_{sequencing_run}/error_correction.qza",
         demux_file = "section_{sequencing_run}/demux_file.qza"
     conda:
         "qiime2-2020.8.0"
+    params:
+        option = demux_single_option
     shell:
         "qiime demux emp-paired "
         "--i-seqs {input.seqs} "
         "--m-barcodes-file {input.barcodes} "
         "--m-barcodes-column BarcodeSequence "
-        "--p-rev-comp-maping-barcodes "
-        "--o-error-correction {output.error_correction} "
+        "{params.option} "
+        "--o-error-correction-details {output.error_correction} "
         "--o-per-sample-sequences {output.demux_file}"
 
 rule demux_dual_barcodes_pheniqs:

mmeds/snakemake/rules/diversity_analysis.smk

@@ -4,7 +4,7 @@ ruleorder: alpha_rarefaction_phylogenetic > alpha_rarefaction
 rule diversity_core_metrics_phylogenetic:
     threads: 10
     input:
-        feature_table = "tables/feature_table.qza",
+        feature_table = "tables/asv_table.qza",
         rooted_tree = "tables/rooted_tree.qza",
         mapping_file = "tables/qiime_mapping_file.tsv"
     output:
@@ -23,7 +23,7 @@ rule diversity_core_metrics_phylogenetic:
 rule diversity_core_metrics:
     threads: 10
     input:
-        feature_table = "tables/feature_table.qza",
+        feature_table = "tables/asv_table.qza",
         mapping_file = "tables/qiime_mapping_file.tsv"
     output:
         "diversity/core_metrics_results"
@@ -39,7 +39,7 @@ rule diversity_core_metrics:
 
 rule alpha_rarefaction_phylogenetic:
     input:
-        feature_table = "tables/feature_table.qza",
+        feature_table = "tables/asv_table.qza",
         rooted_tree = "tables/rooted_tree.qza",
         mapping_file = "tables/qiime_mapping_file.tsv"
     output:
@@ -56,7 +56,7 @@ rule alpha_rarefaction_phylogenetic:
 
 rule alpha_rarefaction:
     input:
-        feature_table = "tables/feature_table.qza",
+        feature_table = "tables/asv_table.qza",
         mapping_file = "tables/qiime_mapping_file.tsv"
     output:
         "diversity/alpha_rarefaction.qzv"

mmeds/snakemake/rules/process_sequencing_data.smk

@@ -3,7 +3,7 @@ rule merge_sequencing_runs:
         feature_tables = expand("section_{sequencing_run}/table_dada2.qza", sequencing_run=config['sequencing_runs']),
         rep_seqs = expand("section_{sequencing_run}/rep_seqs_dada2.qza", sequencing_run=config['sequencing_runs'])
     output:
-        feature_table = "tables/feature_table_no_reads_threshold.qza",
+        feature_table = "tables/asv_table_no_reads_threshold.qza",
         rep_seqs_table = "tables/rep_seqs_table.qza"
     conda:
         "qiime2-2020.8.0"
@@ -13,35 +13,20 @@ rule merge_sequencing_runs:
         qiime feature-table merge-seqs --i-data {input.rep_seqs} --o-merged-data {output.rep_seqs_table}
         """
 
-rule filter_table_by_threshold:
-    input:
-        feature_table = "tables/feature_table_no_reads_threshold.qza",
-        mapping_file = "tables/qiime_mapping_file.tsv"
-    output:
-        "tables/feature_table.qza"
-    conda:
-        "qiime2-2020.8.0"
-    shell:
-        "qiime feature-table filter-samples "
-        "--i-table {input.feature_table} "
-        "--m-metadata-file {input.mapping_file} "
-        "--p-min-frequency {config[sampling_depth]} "
-        "--o-filtered-table {output}"
-
 rule import_single_barcodes:
     input:
         import_dir = "section_{sequencing_run}/import_dir",
         forward_reads = "section_{sequencing_run}/import_dir/forward.fastq.gz",
         reverse_reads = "section_{sequencing_run}/import_dir/reverse.fastq.gz",
         barcodes = "section_{sequencing_run}/import_dir/barcodes.fastq.gz"
     output:
-        "section_{sequencing_run}/qiime_artifact.qza"
+        "section_{sequencing_run}/qiime_import_artifact.qza"
     conda:
         "qiime2-2020.8.0"
     shell:
         "qiime tools import "
         "--type EMPPairedEndSequences "
-        "--input-path {input} "
+        "--input-path {input.import_dir} "
         "--output-path {output}"
 
 rule import_pheniqs_sample_data:
@@ -79,13 +64,13 @@ rule make_pheniqs_config:
         "--mapping-file {input.mapping_file} "
         "--o-directory {input.dir}/pheniqs_output "
         "--o-config {output}"
-        
+
 rule build_phylogenetic_tree:
     input:
-        feature_table = "tables/feature_table.qza",
+        feature_table = "tables/asv_table.qza",
         rep_seqs = "tables/rep_seqs_table.qza"
     output:
-        feature_table_viz = "tables/feature_table_viz.qzv",
+        feature_table_viz = "tables/asv_table_viz.qzv",
         rooted_tree = "tables/rooted_tree.qza"
     conda:
         "qiime2-2020.8.0"
@@ -97,4 +82,3 @@ rule build_phylogenetic_tree:
         qiime phylogeny fasttree --i-alignment temp_files/masked_alignment.qza --o-tree temp_files/unrooted_tree.qza
         qiime phylogeny midpoint-root --i-tree temp_files/unrooted_tree.qza --o-rooted-tree tables/rooted_tree.qza
         """
-

mmeds/snakemake/rules/table_filtering.smk

@@ -1,9 +1,9 @@
 rule filter_table_by_threshold:
     input:
-        feature_table = "tables/feature_table_no_reads_threshold.qza",
+        feature_table = "tables/asv_table_no_reads_threshold.qza",
         mapping_file = "tables/qiime_mapping_file.tsv"
     output:
-        "tables/feature_table.qza"
+        "tables/asv_table.qza"
     conda:
         "qiime2-2020.8.0"
     shell:

mmeds/snakemake/rules/taxonomy.smk

@@ -65,7 +65,7 @@ rule classify_taxonomy_test:
 
 rule taxonomy_collapse:
     input:
-        feature_table = "tables/feature_table.qza",
+        feature_table = "tables/asv_table.qza",
         taxonomy = "tables/taxonomy.qza"
     output:
         "tables/taxa_table_L{level}.qza"
@@ -80,7 +80,7 @@ rule taxonomy_collapse:
 
 rule taxonomic_barplot:
     input:
-        feature_table = "tables/feature_table.qza",
+        feature_table = "tables/asv_table.qza",
         taxonomy = "tables/taxonomy.qza",
         mapping_file = "tables/qiime_mapping_file.tsv"
     output:

mmeds/snakemake/workflows/standard_pipeline.Snakefile

@@ -1,11 +1,13 @@
 configfile: "config_file.yaml"
 report: "report.rst"
 
+include: "{snakemake_dir}/common.smk"
 include: "{snakemake_dir}/demux_denoise.smk"
 include: "{snakemake_dir}/process_sequencing_data.smk"
 include: "{snakemake_dir}/diversity_analysis.smk"
 include: "{snakemake_dir}/taxonomy.smk"
 include: "{snakemake_dir}/differential_abundance.smk"
+include: "{snakemake_dir}/table_filtering.smk"
 
 rule results:
     input:

mmeds/tools/analysis.py

@@ -367,7 +367,7 @@ def initial_setup(self):
     def setup_analysis(self, summary=False):
         """ Setup error logs and jobfile. """
         self.jobtext.append(f"cd {self.run_dir}")
-        self.jobtest.append("ml anaconda3/2024.06")
+        self.jobtext.append("ml anaconda3/2024.06")
         self.jobtext.append("conda activate mmeds_test")
         self.jobtext.append("snakemake --dag | dot -Tpdf >| snakemake_dag.pdf")
         self.jobtext.append("snakemake --rulegraph | dot -Tpdf >| snakemake_rulegraph.pdf")