mrering

.gitignore

@@ -53,3 +53,9 @@ scripts/__pycache__/*
 
 config_361.yml
 config_499.yml
+.DS_Store
+test/QIIME_output/
+test/cluster_logs/
+test/manifest.csv
+test/total_read_counts.tsv
+test/report.RData

Snakefile

@@ -17,7 +17,7 @@ DENOISE_DIR = (QIIME_OUTPUT_DIR + "/denoise" +
                 "_fwd_" + str(trim_left_f) + "-" + str(trunc_len_f) +
                 "_rev_" + str(trim_left_r) + "-" + str(trunc_len_r)
                     )
-CORE_METRIC_DIR = (DENOISE_DIR + "/core-metrics" + 
+CORE_METRIC_DIR = (DENOISE_DIR + "/core-metrics" +
                     "_sampling_depth_" + str(config["diversity"]["sampling_depth"]))
 
 include: "rules/targets/targets.rules"
@@ -30,6 +30,7 @@ include: "rules/tree/tree.rules"
 include: "rules/diversity/diversity.rules"
 include: "rules/unassign/unassign.rules"
 include: "rules/report/report.rules"
+include: "rules/dada2_species/dada2.rules"
 
 workdir: PROJECT_DIR
 

rules/dada2_species/dada2.rules

@@ -3,9 +3,10 @@ rule all_dada2_species:
        TARGET_DADA2SP
 
 rule dada2_assignSpecies:
-    input: DENOISE_DIR + "/representative_seq_fasta/dna-sequences.fasta"
-    params: config["dada2"]["training_set"],
-            config["dada2"]["species_training_set"]
-    output: DADA2_DIR + "/dada2_assigned_species.fasta"
+    input:
+        DENOISE_DIR + "/representative_seq_fasta/dna-sequences.fasta"
+    output:
+        dada2_output = DENOISE_DIR + "/representative_seq_fasta/dada2/dada2_species_assignments.tsv",
+        dada2_rdata = DENOISE_DIR + "/representative_seq_fasta/dada2/dada2.RData"
     script:
-        str(config["dada2"]["rscript"]
+        str(config["dada2"]["rscript"])

rules/report/report.rules

@@ -6,13 +6,13 @@ rule generate_report:
     input:
         mapping_fp = MAPPING_FP,
         demux_count_fp = QIIME_OUTPUT_DIR + "/demux_stat/per-sample-fastq-counts.csv",
-        feature_table_fp = DENOISE_DIR + "/feature_table/feature-table.tsv", 
+        feature_table_fp = DENOISE_DIR + "/feature_table/feature-table.tsv",
         taxo_assignment_fp = DENOISE_DIR + "/taxonomy/taxonomy.tsv",
         faith_fp = CORE_METRIC_DIR + "/faith_pd.tsv",
         uu_fp = CORE_METRIC_DIR + "/uu.tsv",
         wu_fp = CORE_METRIC_DIR + "/wu.tsv"
     output:
         CORE_METRIC_DIR + "/basic_bioinformatics_report/Basic_Bioinformatics_Report.html"
     script:
-        str(config["report"]["rmd"]) 
+        str(config["report"]["rmd"])
 

rules/targets/targets.rules

@@ -20,6 +20,9 @@ TARGET_UNASSIGN = [str(DENOISE_DIR + "/representative_seq_fasta/unassigned/unass
 
 TARGET_REPORT = [str(CORE_METRIC_DIR + "/basic_bioinformatics_report/Basic_Bioinformatics_Report.html")]
 
+TARGET_DADA2SP = [str(DENOISE_DIR + "/representative_seq_fasta/dada2/dada2_species_assignments.tsv"),
+                  str(DENOISE_DIR + "/representative_seq_fasta/dada2/dada2.RData")]
+
 TARGET_ALL = (
     TARGET_DNABC +
     TARGET_QIIME_IMPORT +
@@ -30,5 +33,6 @@ TARGET_ALL = (
     TARGET_TREE +
     TARGET_DIVERSITY +
     TARGET_UNASSIGN +
-    TARGET_REPORT
+    TARGET_REPORT +
+    TARGET_DADA2SP
 )

scripts/Basic_Bioinformatics_Report.Rmd

@@ -422,3 +422,10 @@ ReadCounts <- s1 %>% select(SampleID, demux_Read_Counts, denoise_Read_Counts, QC
   arrange(QC_Read_Counts)
 datatable(ReadCounts, fillContainer = FALSE, options = list(pageLength = 10))
 ```
+
+```{r save environment}
+
+save.image(file = snakemake@config[["report"]][["report_rdata"]])
+
+```
+

scripts/dada2.R

@@ -0,0 +1,35 @@
+#!/usr/bin/env Rscript
+
+library(tidyverse)
+library(Biostrings)
+library(dada2)
+
+#from config.yml
+
+rdp_species <- snakemake@config[["dada2"]][["species_training_set"]]
+
+#from dada2.rules
+
+my_seqs <- readDNAStringSet(snakemake@input[[1]])
+
+#main
+
+seqs_df <- tibble(query_id = names(my_seqs),
+                  sequence = as.character(my_seqs))
+
+set.seed(100) # Initialize random number generator for reproducibility
+
+genus.species <- assignSpecies(seqs_df$sequence, rdp_species)
+
+genspec_df <- genus.species %>%
+    as.data.frame() %>%
+    rownames_to_column(var = "sequence") %>%
+    inner_join(seqs_df, by = "sequence") %>%
+    mutate(species = paste0(Genus, " ", Species)) %>% #making same as unassigner output
+    select(query_id, species, sequence)
+
+#from dada2.rules and targets.rules
+
+save.image(file = snakemake@output[["dada2_rdata"]])
+
+write_tsv(x = genspec_df, path = snakemake@output[["dada2_output"]])