Add an option to use different correlation methods

Author: ahmig

.test/targets.yaml

@@ -33,6 +33,9 @@ DEMIX:
   COV_CUTOFF: 0
   MIN_ABUNDANCE: 0.0001
 USE_BIONJ: false
+COR:
+  METHOD: "pearson"
+  EXACT: null
 PLOT_GENOME_REGIONS:
   "config/nsp_annotation.csv"
 REPORT_QMD:

config/README.md

@@ -140,6 +140,9 @@ All of the following variables are pre-defined in [config.yaml](/config/config.y
   - `ACCESSION_COLUMN`: name of the column that contains GISAID EPI identifiers in the input target metadata.
 - `DIVERSITY_REPS`: number of random sample subsets of the context dataset for the nucleotide diversity comparison.
 - `USE_BIONJ`: use the BIONJ algorithm ([Gascuel, 1997](https://doi.org/10.1093/oxfordjournals.molbev.a025808)) instead of NJ (neighbor-joining; [Saitou & Nei, 1987](https://doi.org/10.1093/oxfordjournals.molbev.a040454)) to reconstruct phylogenetic trees from pairwise distances.
+- `COR`: configuration for correlation analyses of allele frequency data over time and between variants. This parameter controls how correlation tests are performed using R's `cor.test` and `cor` functions (see [R documentation](https://search.r-project.org/CRAN/refmans/correlation/html/cor_test.html)).
+  - `METHOD`: correlation method to use. Valid options are "pearson" (default), "kendall", or "spearman".
+  - `EXACT`: boolean flag indicating whether to compute an exact p-value when possible. This option applies only to certain methods and may be set to `null` (default) to let R decide automatically.
 - `LOG_PY_FMT`: logging format string for Python scripts.
 - `PLOTS`: path of the R script that sets the design and style of data visualizations.
 - `PLOT_GENOME_REGIONS`: path of a CSV file containing genome regions, e.g. SARS-CoV-2 non-structural protein (NSP) coordinates, for data visualization.

config/config.yaml

@@ -33,6 +33,9 @@ GISAID:
 DIVERSITY_REPS:
   1000
 USE_BIONJ: false
+COR:
+  METHOD: "pearson"
+  EXACT: null
 LOG_PY_FMT:
   "%(asctime)s - %(name)s - %(levelname)-8s - %(message)s"
 PLOTS:

template.qmd

@@ -23,6 +23,7 @@ params:
   min_ivar_freq: ""
   workflow_version: ""
   use_bionj: ""
+  cor_method: ""
   div: ""
   freyja: ""
   tree: ""
@@ -81,8 +82,17 @@ n.samples <- table %>% pull(Sample) %>% unique() %>% length()
 vcf <- read_csv(params$heat_tab) 
 row.names(vcf) <- vcf$`...1`
 vcf <- vcf[-1]
-cor.mat <- cor(vcf)
+cor.mat <- cor(vcf, method = params$cor_method)
 cor.mat[is.na(cor.mat)] <- 0
+if (params$cor_method == "pearson") {
+  cor.method.name <- "Pearson's"
+} else if (params$cor_method == "spearman") {
+  cor.method.name <- "Spearman's"
+} else if (params$cor_method == "kendall") {
+  cor.method.name <- "Kendall's"
+} else {
+  cor.method.name <- "Undetermined"
+}
 
 # Distance tree
 if (params$use_bionj) {
@@ -217,7 +227,7 @@ at the bottom, and the latest, at the top.](`r params$SNV_s`){#fig-SNV_s}
 The correlation of the allele frequency of each NV with the time since the
 initial sampling has been calculated (@fig-volcano).
 
-![Pearson’s correlation coefficients and adjusted p-values of
+![`r cor.method.name` correlation coefficients and adjusted p-values of
 allele frequencies with time. Red dashed line indicates adjusted $p = 0.05$.
 Labeled dots represent nucleotide variants correlated with time
 (adjusted $p < 0.05$).](`r params$volcano`){#fig-volcano}
@@ -253,7 +263,7 @@ zooming in on specific regions.
 
 ```{r heatmap, echo = F, message = F, fig.align = 'center'}
 #| label: fig-heatmap
-#| fig-cap: "Interactive hierarchically clustered heatmap of the pairwise Pearson’s correlation coefficients between the time series of allele frequencies in the case study."
+#| fig-cap: "Interactive hierarchically clustered heatmap of the pairwise correlation coefficients between the time series of allele frequencies in the case study."
 
 heatmaply_cor(
   cor.mat,

workflow/rules/report.smk

@@ -140,7 +140,9 @@ rule evo_plots:
 rule snp_plots:
     conda: "../envs/renv.yaml"
     params:
-        design = config["PLOTS"]
+        design = config["PLOTS"],
+        cor_method = config["COR"]["METHOD"],
+        cor_exact = config["COR"]["EXACT"]
     input:
         vcf =  OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv",
         metadata = config["METADATA"]
@@ -196,7 +198,8 @@ rule report:
         ufboot_reps = config["UFBOOT_REPS"],
         shalrt_reps = config["SHALRT_REPS"],
         name = config["OUTPUT_NAME"],
-        use_bionj = config["USE_BIONJ"]
+        use_bionj = config["USE_BIONJ"],
+        cor_method = config["COR"]["METHOD"]
     output:
         html = report(OUTDIR/f"{OUTPUT_NAME}.report.html")
     log:
@@ -211,6 +214,7 @@ rule report:
                                                        min_ivar_freq='{params.min_ivar_freq}',\
                                                        workflow_version='{params.workflow_version}',\
                                                        use_bionj='{params.use_bionj}',\
+                                                       cor_method='{params.cor_method}',\
                                                        div='{input.diversity}',\
                                                        freyja ='{input.freyja}',\
                                                        tree = '{input.tree}',\

workflow/scripts/report/snp_plots.R

@@ -93,7 +93,9 @@ cor.df.fill <- lapply(
 
     test <- cor.test(
       df$ALT_FREQ,
-      df$interval
+      df$interval,
+      method = snakemake@params$cor_method,
+      exact = snakemake@params$cor_exact
     )
 
     pvalue <- p.adjust(
@@ -174,7 +176,7 @@ subset <- c(sign, dup) %>%
 # Plot height depending on the number of SNPs assuming 4 columns in the plot
 plot.height <- ceiling(length(subset) / 4) * 42
 
-log_info("PLotting SNPs trends in time")
+log_info("Plotting SNPs trends in time")
 panel <- vcf %>%
   filter(variant %in% subset) %>%
   ggplot() +
@@ -221,7 +223,7 @@ log_info("Saving coorelation table")
 cor.df.fill %>%
   transmute(
     NV = snp,
-    PearsonCor = cor,
+    correlation_coef = cor,
     adj_pvalue = p.value
   ) %>%
   write.csv(