From 4559c59ea5ffd9b9098776a42e48bd07cbab9b46 Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 16 Jan 2026 20:46:08 -0800 Subject: [PATCH 01/14] initial commit of deseq2 function --- single_cell/dge/deseq2_function.R | 64 +++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 single_cell/dge/deseq2_function.R diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R new file mode 100644 index 0000000..eb4deb5 --- /dev/null +++ b/single_cell/dge/deseq2_function.R @@ -0,0 +1,64 @@ +library(edgeR) +library(DESeq2) +library(tidyverse) +##### DGE analysis using DESeq2 +#' @param mat: feature (row) x sample (column) raw count matrix +#' @param meta: sample metadata, with rownames containing column names of mat +#' @param comparison: name of meta column to use DGE comparion. Must have exactly two levels +#' @param numerator: level to be used as numerator in log2FC calculation +#' @param denominator: level to be used as denominator in log2FC calculation +#' @param fixed_effects: a vector of meta column names to use as fixed effects +run_deseq2 <- function(mat, meta, comparison, numerator, denominator, fixed_effects=NULL) { + if(! all( colnames(mat) %in% rownames(meta)) ) { + stop("One or more samples in mat are not present in meta.") + } + # Subset and order metadata according to the order of samples in mat + meta = meta[ match(colnames(mat), rownames(meta)) , ] + + if( length(unique(as.vector(meta[,comparison]))) != 2) { + stop("Comparison column of meta does not have exactly two unique values.") + } + if( ! numerator %in% meta[,comparison] ) { + stop("Numerator value does not exist in comparison column of meta.") + } + if( ! denominator %in% meta[,comparison] ) { + stop("Denominator value does not exist in comparison column of meta.") + } + if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(meta)) ) { + stop("One or more fixed_effects do not exist in meta.") + } + + # Keep features that are expressed with CPM > 1 in more than 50% of samples in at least one comparison group. + # feats = names(which(apply(cpm(mat) > 1, 1, function(x) { + # any(tapply(x, meta[,comparison], mean) > 0.5) + # } ) + # ) + # ) + + # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least 60% of min_gsize samples, regardless of the sample group identity. + min_gsize = min(table(meta[,comparison])) + feats = names(which(rowSums(cpm(mat) > 1) > (min_gsize * 0.6))) + print(paste0("Number of final features: ", length(feats))) + mat = mat[feats,] + + # Make formula string + formula_str = paste("~", comparison) + if(! is.null(fixed_effects)) { + for(fe in fixed_effects) { + formula_str = paste(formula_str, "+", fe) + } + } + + # Perform DGE analysis + dds = DESeqDataSetFromMatrix(mat, colData = meta, design = as.formula( formula_str )) + dds = DESeq(dds) + + # Extract and format DGE results + res = results(dds, contrast = c(comparison, numerator, denominator)) %>% as.data.frame() + res = arrange(res, padj) + + res = dplyr::rename(res, "logFC"="log2FoldChange", "AveExpr"="baseMean", "P.Value"="pvalue", "adj.P.Val"="padj") + return(res) +} + + From 8eb4c51f042f13931b1587ec56fc4bb470c7ceaa Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 27 Feb 2026 11:11:02 -0800 Subject: [PATCH 02/14] updated input parameter names --- single_cell/dge/deseq2_function.R | 54 +++++++++++++++---------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index eb4deb5..72778c1 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -2,47 +2,47 @@ library(edgeR) library(DESeq2) library(tidyverse) ##### DGE analysis using DESeq2 -#' @param mat: feature (row) x sample (column) raw count matrix -#' @param meta: sample metadata, with rownames containing column names of mat -#' @param comparison: name of meta column to use DGE comparion. Must have exactly two levels -#' @param numerator: level to be used as numerator in log2FC calculation -#' @param denominator: level to be used as denominator in log2FC calculation -#' @param fixed_effects: a vector of meta column names to use as fixed effects -run_deseq2 <- function(mat, meta, comparison, numerator, denominator, fixed_effects=NULL) { - if(! all( colnames(mat) %in% rownames(meta)) ) { - stop("One or more samples in mat are not present in meta.") +#' @param counts: feature (row) x sample (column) raw count matrix +#' @param metadata: sample metadata, with rownames containing column names of \code{counts} +#' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels +#' @param case_group: level to be used as numerator in log2FC calculation +#' @param reference_group: level to be used as denominator in log2FC calculation +#' @param fixed_effects: a vector of \code{metadata} column names to use as fixed effects +run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL) { + if(! all( colnames(counts) %in% rownames(metadata)) ) { + stop("One or more samples in counts are not present in metadata.") } - # Subset and order metadata according to the order of samples in mat - meta = meta[ match(colnames(mat), rownames(meta)) , ] + # Subset and order metadata according to the order of samples in counts + metadata = metadata[ match(colnames(counts), rownames(metadata)) , ] - if( length(unique(as.vector(meta[,comparison]))) != 2) { - stop("Comparison column of meta does not have exactly two unique values.") + if( length(unique(as.vector(metadata[,dge_by]))) != 2) { + stop("Comparison column of metadata does not have exactly two unique values.") } - if( ! numerator %in% meta[,comparison] ) { - stop("Numerator value does not exist in comparison column of meta.") + if( ! case_group %in% metadata[,dge_by] ) { + stop("Case group does not exist in dge_by column of metadata.") } - if( ! denominator %in% meta[,comparison] ) { - stop("Denominator value does not exist in comparison column of meta.") + if( ! reference_group %in% metadata[,dge_by] ) { + stop("Reference group does not exist in dge_by column of metadata.") } - if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(meta)) ) { - stop("One or more fixed_effects do not exist in meta.") + if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(metadata)) ) { + stop("One or more fixed_effects do not exist in metadata.") } # Keep features that are expressed with CPM > 1 in more than 50% of samples in at least one comparison group. - # feats = names(which(apply(cpm(mat) > 1, 1, function(x) { - # any(tapply(x, meta[,comparison], mean) > 0.5) + # feats = names(which(apply(cpm(counts) > 1, 1, function(x) { + # any(tapply(x, metadata[,dge_by], mean) > 0.5) # } ) # ) # ) # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least 60% of min_gsize samples, regardless of the sample group identity. - min_gsize = min(table(meta[,comparison])) - feats = names(which(rowSums(cpm(mat) > 1) > (min_gsize * 0.6))) + min_gsize = min(table(metadata[,dge_by])) + feats = names(which(rowSums(cpm(counts) > 1) > (min_gsize * 0.6))) print(paste0("Number of final features: ", length(feats))) - mat = mat[feats,] + counts = counts[feats,] # Make formula string - formula_str = paste("~", comparison) + formula_str = paste("~", dge_by) if(! is.null(fixed_effects)) { for(fe in fixed_effects) { formula_str = paste(formula_str, "+", fe) @@ -50,11 +50,11 @@ run_deseq2 <- function(mat, meta, comparison, numerator, denominator, fixed_effe } # Perform DGE analysis - dds = DESeqDataSetFromMatrix(mat, colData = meta, design = as.formula( formula_str )) + dds = DESeqDataSetFromMatrix(counts, colData = metadata, design = as.formula( formula_str )) dds = DESeq(dds) # Extract and format DGE results - res = results(dds, contrast = c(comparison, numerator, denominator)) %>% as.data.frame() + res = results(dds, contrast = c(dge_by, case_group, reference_group)) %>% as.data.frame() res = arrange(res, padj) res = dplyr::rename(res, "logFC"="log2FoldChange", "AveExpr"="baseMean", "P.Value"="pvalue", "adj.P.Val"="padj") From 99ab50edd0686a2acf29e9467d62ec1e93a58fbe Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 27 Feb 2026 11:45:07 -0800 Subject: [PATCH 03/14] Modified deseq2 function to remove input checks and added an input_checks_function.R --- single_cell/dge/deseq2_function.R | 34 +++++++------------------ single_cell/dge/input_checks_function.R | 28 ++++++++++++++++++++ 2 files changed, 37 insertions(+), 25 deletions(-) create mode 100644 single_cell/dge/input_checks_function.R diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index 72778c1..c21e62a 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -1,7 +1,9 @@ library(edgeR) library(DESeq2) library(tidyverse) -##### DGE analysis using DESeq2 +### DGE analysis using DESeq2 +# +#' Run DGE between 2 groups using DESeq2 #' @param counts: feature (row) x sample (column) raw count matrix #' @param metadata: sample metadata, with rownames containing column names of \code{counts} #' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels @@ -9,31 +11,13 @@ library(tidyverse) #' @param reference_group: level to be used as denominator in log2FC calculation #' @param fixed_effects: a vector of \code{metadata} column names to use as fixed effects run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL) { - if(! all( colnames(counts) %in% rownames(metadata)) ) { - stop("One or more samples in counts are not present in metadata.") - } - # Subset and order metadata according to the order of samples in counts - metadata = metadata[ match(colnames(counts), rownames(metadata)) , ] - - if( length(unique(as.vector(metadata[,dge_by]))) != 2) { - stop("Comparison column of metadata does not have exactly two unique values.") - } - if( ! case_group %in% metadata[,dge_by] ) { - stop("Case group does not exist in dge_by column of metadata.") - } - if( ! reference_group %in% metadata[,dge_by] ) { - stop("Reference group does not exist in dge_by column of metadata.") - } - if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(metadata)) ) { - stop("One or more fixed_effects do not exist in metadata.") - } - # Keep features that are expressed with CPM > 1 in more than 50% of samples in at least one comparison group. - # feats = names(which(apply(cpm(counts) > 1, 1, function(x) { - # any(tapply(x, metadata[,dge_by], mean) > 0.5) - # } ) - # ) - # ) + + # Commong input checks + metadata = input_checks(counts, metadata, dge_by, case_group, reference_group, fixed_effects) + + # Method specific input checks + # None for DESeq2 # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least 60% of min_gsize samples, regardless of the sample group identity. min_gsize = min(table(metadata[,dge_by])) diff --git a/single_cell/dge/input_checks_function.R b/single_cell/dge/input_checks_function.R new file mode 100644 index 0000000..bc3394e --- /dev/null +++ b/single_cell/dge/input_checks_function.R @@ -0,0 +1,28 @@ +### Check inputs to DGE functions +# +#' Check inputs to DGE functions +#' @param counts: feature (row) x sample (column) raw count matrix. Users must include only those samples in \code{counts} that should be used in DGE analysis +#' @param metadata: sample metadata, with rownames containing column names of \code{counts} +#' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels +#' @param case_group: level to be used as numerator in log2FC calculation +#' @param reference_group: level to be used as denominator in log2FC calculation +#' @param fixed_effects: a vector of \code{metadata} column names to use as fixed effects +#' @param random_effects: a vector of \code{metadata} column names to use as random effects +#' @return \code{metadata} after matching the rows to the columns (samples) in \code{counts} +input_checks <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL, random_effects=NULL) { + if (!is.matrix(counts) & !is.data.frame(counts) ) stop("Error. Input data is not a matrix or a data.frame.") + if (ncol(counts)==0 | nrow(counts)==0) stop("Error. Count matrix does not have at least one column and row.") + if(! dge_by %in% colnames(metadata) ) stop("Error. dge_by colum does not exist in metadata.") + if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(metadata)) ) stop("Error. One or more fixed_effects do not exist in metadata.") + if(! is.null(random_effects) & ! all(random_effects %in% colnames(metadata)) ) stop("Error. One or more random_effects do not exist in metadata.") + if(! all( colnames(counts) %in% rownames(metadata)) ) stop("Error. One or more samples in counts are not present in metadata rownames.") + + # Subset and order metadata according to the order of samples in counts + metadata = metadata[ match(colnames(counts), rownames(metadata)) , ] + + if( length(unique(as.vector(metadata[,dge_by]))) != 2) stop("Error. dge_by column of metadata does not have exactly two levels (unique values).") + if( ! case_group %in% metadata[,dge_by] ) stop("Case group does not exist in dge_by column of metadata.") + if( ! reference_group %in% metadata[,dge_by] ) stop("Reference group does not exist in dge_by column of metadata.") + return(metadata) +} + From 2a4ddfd2b02911f4d8c07fd45aae1aaa1620bd7e Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 27 Feb 2026 11:55:50 -0800 Subject: [PATCH 04/14] Misc changes in the text --- single_cell/dge/deseq2_function.R | 2 +- single_cell/dge/input_checks_function.R | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index c21e62a..b8bbe1e 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -4,7 +4,7 @@ library(tidyverse) ### DGE analysis using DESeq2 # #' Run DGE between 2 groups using DESeq2 -#' @param counts: feature (row) x sample (column) raw count matrix +#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} #' @param metadata: sample metadata, with rownames containing column names of \code{counts} #' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels #' @param case_group: level to be used as numerator in log2FC calculation diff --git a/single_cell/dge/input_checks_function.R b/single_cell/dge/input_checks_function.R index bc3394e..69b5823 100644 --- a/single_cell/dge/input_checks_function.R +++ b/single_cell/dge/input_checks_function.R @@ -1,7 +1,7 @@ ### Check inputs to DGE functions # -#' Check inputs to DGE functions -#' @param counts: feature (row) x sample (column) raw count matrix. Users must include only those samples in \code{counts} that should be used in DGE analysis +#' Check inputs to DGE functions, match \code{metadata} rows to \code{counts} columns, and return \code{metadata}. +#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} #' @param metadata: sample metadata, with rownames containing column names of \code{counts} #' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels #' @param case_group: level to be used as numerator in log2FC calculation From 6aa13c735da0d4c5bbaec7b24d9d29c67de671c0 Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 27 Feb 2026 12:07:34 -0800 Subject: [PATCH 05/14] Add function to remove lowly expressed genes --- single_cell/dge/deseq2_function.R | 6 +----- .../dge/remove_low_expression_function.R | 17 +++++++++++++++++ 2 files changed, 18 insertions(+), 5 deletions(-) create mode 100644 single_cell/dge/remove_low_expression_function.R diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index b8bbe1e..05f3730 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -19,11 +19,7 @@ run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fi # Method specific input checks # None for DESeq2 - # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least 60% of min_gsize samples, regardless of the sample group identity. - min_gsize = min(table(metadata[,dge_by])) - feats = names(which(rowSums(cpm(counts) > 1) > (min_gsize * 0.6))) - print(paste0("Number of final features: ", length(feats))) - counts = counts[feats,] + counts = remove_low_expression(counts, metadata, dge_by, min_pct=60) # Make formula string formula_str = paste("~", dge_by) diff --git a/single_cell/dge/remove_low_expression_function.R b/single_cell/dge/remove_low_expression_function.R new file mode 100644 index 0000000..abd6616 --- /dev/null +++ b/single_cell/dge/remove_low_expression_function.R @@ -0,0 +1,17 @@ +### Remove lowly expressed genes from the count matrix. +# +#' Remove lowly expressed genes from the count matrix. +#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} +#' @param metadata: sample metadata, with rownames containing column names of \code{counts} +#' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels +#' @param min_pct: Genes with counts per million (CPM) more than 1 in at least \code{min_pct} % of samples of the smaller group are retained. Default 60% +#' @return \code{counts} after removing lowly expressed genes +remove_low_expression <- function(counts, metadata, dge_by, min_pct=60) { + # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least 60% of min_gsize samples, regardless of the sample group identity. + min_gsize = min(table(metadata[,dge_by])) + feats = names(which(rowSums(cpm(counts) > 1) > (min_gsize * min_pct/100))) + cat(paste0("Number of final features: ", length(feats)), "\n") + counts = counts[feats,] + return(counts) +} + From e86c01dd306888fb68070a0d7e2c07128ea5d9f6 Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 27 Feb 2026 17:00:55 -0800 Subject: [PATCH 06/14] replace min_pct by min_frac and move helper functions --- .../helper_functions/input_checks_function.R | 28 +++++++++++++++++++ .../remove_low_expression_function.R | 17 +++++++++++ 2 files changed, 45 insertions(+) create mode 100644 single_cell/dge/helper_functions/input_checks_function.R create mode 100644 single_cell/dge/helper_functions/remove_low_expression_function.R diff --git a/single_cell/dge/helper_functions/input_checks_function.R b/single_cell/dge/helper_functions/input_checks_function.R new file mode 100644 index 0000000..69b5823 --- /dev/null +++ b/single_cell/dge/helper_functions/input_checks_function.R @@ -0,0 +1,28 @@ +### Check inputs to DGE functions +# +#' Check inputs to DGE functions, match \code{metadata} rows to \code{counts} columns, and return \code{metadata}. +#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} +#' @param metadata: sample metadata, with rownames containing column names of \code{counts} +#' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels +#' @param case_group: level to be used as numerator in log2FC calculation +#' @param reference_group: level to be used as denominator in log2FC calculation +#' @param fixed_effects: a vector of \code{metadata} column names to use as fixed effects +#' @param random_effects: a vector of \code{metadata} column names to use as random effects +#' @return \code{metadata} after matching the rows to the columns (samples) in \code{counts} +input_checks <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL, random_effects=NULL) { + if (!is.matrix(counts) & !is.data.frame(counts) ) stop("Error. Input data is not a matrix or a data.frame.") + if (ncol(counts)==0 | nrow(counts)==0) stop("Error. Count matrix does not have at least one column and row.") + if(! dge_by %in% colnames(metadata) ) stop("Error. dge_by colum does not exist in metadata.") + if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(metadata)) ) stop("Error. One or more fixed_effects do not exist in metadata.") + if(! is.null(random_effects) & ! all(random_effects %in% colnames(metadata)) ) stop("Error. One or more random_effects do not exist in metadata.") + if(! all( colnames(counts) %in% rownames(metadata)) ) stop("Error. One or more samples in counts are not present in metadata rownames.") + + # Subset and order metadata according to the order of samples in counts + metadata = metadata[ match(colnames(counts), rownames(metadata)) , ] + + if( length(unique(as.vector(metadata[,dge_by]))) != 2) stop("Error. dge_by column of metadata does not have exactly two levels (unique values).") + if( ! case_group %in% metadata[,dge_by] ) stop("Case group does not exist in dge_by column of metadata.") + if( ! reference_group %in% metadata[,dge_by] ) stop("Reference group does not exist in dge_by column of metadata.") + return(metadata) +} + diff --git a/single_cell/dge/helper_functions/remove_low_expression_function.R b/single_cell/dge/helper_functions/remove_low_expression_function.R new file mode 100644 index 0000000..93718ec --- /dev/null +++ b/single_cell/dge/helper_functions/remove_low_expression_function.R @@ -0,0 +1,17 @@ +### Remove lowly expressed genes from the count matrix. +# +#' Remove lowly expressed genes from the count matrix. +#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} +#' @param metadata sample metadata, with rownames containing column names of \code{counts} +#' @param dge_by name of \code{metadata} column to use DGE comparion. Must have exactly two levels +#' @param min_frac Genes with counts per million (CPM) more than 1 in at least \code{min_frac} fraction of samples of the smaller group are retained. Default 0.6 (i.e., 60%) +#' @return \code{counts} after removing lowly expressed genes +remove_low_expression <- function(counts, metadata, dge_by, min_frac=0.6) { + # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least min_frac of min_gsize samples, regardless of the sample group identity. + min_gsize = min(table(metadata[ , dge_by])) + feats = names(which(rowSums(cpm(counts) > 1) > (min_gsize * min_frac))) + cat(paste0("Number of final features: ", length(feats)), "\n") + counts = counts[feats,] + return(counts) +} + From 21392bcff10ce4900d27cb280202ac3cf8ad2830 Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 27 Feb 2026 17:03:15 -0800 Subject: [PATCH 07/14] deleted helper function files from dge directory --- single_cell/dge/input_checks_function.R | 28 ------------------- .../dge/remove_low_expression_function.R | 17 ----------- 2 files changed, 45 deletions(-) delete mode 100644 single_cell/dge/input_checks_function.R delete mode 100644 single_cell/dge/remove_low_expression_function.R diff --git a/single_cell/dge/input_checks_function.R b/single_cell/dge/input_checks_function.R deleted file mode 100644 index 69b5823..0000000 --- a/single_cell/dge/input_checks_function.R +++ /dev/null @@ -1,28 +0,0 @@ -### Check inputs to DGE functions -# -#' Check inputs to DGE functions, match \code{metadata} rows to \code{counts} columns, and return \code{metadata}. -#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} -#' @param metadata: sample metadata, with rownames containing column names of \code{counts} -#' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels -#' @param case_group: level to be used as numerator in log2FC calculation -#' @param reference_group: level to be used as denominator in log2FC calculation -#' @param fixed_effects: a vector of \code{metadata} column names to use as fixed effects -#' @param random_effects: a vector of \code{metadata} column names to use as random effects -#' @return \code{metadata} after matching the rows to the columns (samples) in \code{counts} -input_checks <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL, random_effects=NULL) { - if (!is.matrix(counts) & !is.data.frame(counts) ) stop("Error. Input data is not a matrix or a data.frame.") - if (ncol(counts)==0 | nrow(counts)==0) stop("Error. Count matrix does not have at least one column and row.") - if(! dge_by %in% colnames(metadata) ) stop("Error. dge_by colum does not exist in metadata.") - if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(metadata)) ) stop("Error. One or more fixed_effects do not exist in metadata.") - if(! is.null(random_effects) & ! all(random_effects %in% colnames(metadata)) ) stop("Error. One or more random_effects do not exist in metadata.") - if(! all( colnames(counts) %in% rownames(metadata)) ) stop("Error. One or more samples in counts are not present in metadata rownames.") - - # Subset and order metadata according to the order of samples in counts - metadata = metadata[ match(colnames(counts), rownames(metadata)) , ] - - if( length(unique(as.vector(metadata[,dge_by]))) != 2) stop("Error. dge_by column of metadata does not have exactly two levels (unique values).") - if( ! case_group %in% metadata[,dge_by] ) stop("Case group does not exist in dge_by column of metadata.") - if( ! reference_group %in% metadata[,dge_by] ) stop("Reference group does not exist in dge_by column of metadata.") - return(metadata) -} - diff --git a/single_cell/dge/remove_low_expression_function.R b/single_cell/dge/remove_low_expression_function.R deleted file mode 100644 index abd6616..0000000 --- a/single_cell/dge/remove_low_expression_function.R +++ /dev/null @@ -1,17 +0,0 @@ -### Remove lowly expressed genes from the count matrix. -# -#' Remove lowly expressed genes from the count matrix. -#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} -#' @param metadata: sample metadata, with rownames containing column names of \code{counts} -#' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels -#' @param min_pct: Genes with counts per million (CPM) more than 1 in at least \code{min_pct} % of samples of the smaller group are retained. Default 60% -#' @return \code{counts} after removing lowly expressed genes -remove_low_expression <- function(counts, metadata, dge_by, min_pct=60) { - # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least 60% of min_gsize samples, regardless of the sample group identity. - min_gsize = min(table(metadata[,dge_by])) - feats = names(which(rowSums(cpm(counts) > 1) > (min_gsize * min_pct/100))) - cat(paste0("Number of final features: ", length(feats)), "\n") - counts = counts[feats,] - return(counts) -} - From 6da9ede4d4ab8bb0b965faffb3b2b254f7503833 Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 27 Feb 2026 17:37:14 -0800 Subject: [PATCH 08/14] using dge_groups in removing lowly expressed genes --- .../remove_low_expression_function.R | 27 ++++++++++++------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/single_cell/dge/helper_functions/remove_low_expression_function.R b/single_cell/dge/helper_functions/remove_low_expression_function.R index 93718ec..44b4777 100644 --- a/single_cell/dge/helper_functions/remove_low_expression_function.R +++ b/single_cell/dge/helper_functions/remove_low_expression_function.R @@ -1,17 +1,26 @@ ### Remove lowly expressed genes from the count matrix. # -#' Remove lowly expressed genes from the count matrix. -#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} -#' @param metadata sample metadata, with rownames containing column names of \code{counts} -#' @param dge_by name of \code{metadata} column to use DGE comparion. Must have exactly two levels -#' @param min_frac Genes with counts per million (CPM) more than 1 in at least \code{min_frac} fraction of samples of the smaller group are retained. Default 0.6 (i.e., 60%) -#' @return \code{counts} after removing lowly expressed genes -remove_low_expression <- function(counts, metadata, dge_by, min_frac=0.6) { +#' Remove lowly expressed genes from the count matrix. The function first determines the size of the smallest group (denoted min_gsize), where groups are defined by \code{dge_by}. A gene is retaind only if it is expressed (CPM > 1) in at least a \code{min_frac} proportion of min_gsize samples overall (i.e., not within each group separately). +#' @param counts Feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts}. +#' @param metadata Sample metadata, with rownames containing column names of \code{counts}. +#' @param dge_by Name of \code{metadata} column to use DGE comparion. Must have exactly two levels. +#' @param dge_groups A vector of group labels (from \code{dge_by}) to consider. Only samples belonging to these groups are used to determine the minimum number of samples required to have CPM > 1. If NULL, all groups in \code{dge_by} are considered. +#' @param min_frac Genes with counts per million (CPM) more than 1 in at least \code{min_frac} proportion of samples of the smaller group are retained, regardless of group membership. Default 0.6 (i.e., 60%). +#' @return \code{counts} after removing lowly expressed genes. +remove_low_expression <- function(counts, metadata, dge_by, dge_groups=NULL, min_frac=0.6) { # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least min_frac of min_gsize samples, regardless of the sample group identity. - min_gsize = min(table(metadata[ , dge_by])) - feats = names(which(rowSums(cpm(counts) > 1) > (min_gsize * min_frac))) + warning("Assuming that the order of samples (columns) in 'counts' matches the order of samples (rows) in 'metadata'.") + if(is.null(dge_groups)) { + warning("dge_groups is not provided. Using all groups in dge_by to determine the minimum number of samples required to have CPM > 1.") + sample_idx = rep(TRUE, nrow(metadata)) + } else { + sample_idx = metadata[, dge_by] %in% dge_groups + } + min_gsize = min(table(metadata[ sample_idx, dge_by])) + feats = names(which(rowSums(cpm(counts[, sample_idx]) > 1) > (min_gsize * min_frac))) cat(paste0("Number of final features: ", length(feats)), "\n") counts = counts[feats,] return(counts) } + From 263777718afd295f96a3524f3dd4418713061a10 Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 27 Feb 2026 17:38:20 -0800 Subject: [PATCH 09/14] replacing min_pct by min_frac in deseq2_function --- single_cell/dge/deseq2_function.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index 05f3730..8102f6b 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -12,14 +12,13 @@ library(tidyverse) #' @param fixed_effects: a vector of \code{metadata} column names to use as fixed effects run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL) { - # Commong input checks metadata = input_checks(counts, metadata, dge_by, case_group, reference_group, fixed_effects) # Method specific input checks # None for DESeq2 - counts = remove_low_expression(counts, metadata, dge_by, min_pct=60) + counts = remove_low_expression(counts, metadata, dge_by, min_frac=0.6) # Make formula string formula_str = paste("~", dge_by) From 59c56f11adcbafd0fbcff1153c9ecc279a8ddec4 Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Mon, 2 Mar 2026 14:13:28 -0800 Subject: [PATCH 10/14] Updated internal function names to start with a dot --- .../helper_functions/input_checks_function.R | 18 +++++++++--------- .../remove_low_expression_function.R | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/single_cell/dge/helper_functions/input_checks_function.R b/single_cell/dge/helper_functions/input_checks_function.R index 69b5823..1fffc57 100644 --- a/single_cell/dge/helper_functions/input_checks_function.R +++ b/single_cell/dge/helper_functions/input_checks_function.R @@ -1,15 +1,15 @@ ### Check inputs to DGE functions # #' Check inputs to DGE functions, match \code{metadata} rows to \code{counts} columns, and return \code{metadata}. -#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} -#' @param metadata: sample metadata, with rownames containing column names of \code{counts} -#' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels -#' @param case_group: level to be used as numerator in log2FC calculation -#' @param reference_group: level to be used as denominator in log2FC calculation -#' @param fixed_effects: a vector of \code{metadata} column names to use as fixed effects -#' @param random_effects: a vector of \code{metadata} column names to use as random effects -#' @return \code{metadata} after matching the rows to the columns (samples) in \code{counts} -input_checks <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL, random_effects=NULL) { +#' @param counts A feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts}. +#' @param metadata Sample metadata, with rownames containing column names of \code{counts}. +#' @param dge_by Name of \code{metadata} column to use DGE comparion. Must have exactly two levels. +#' @param case_group Group to be used as numerator in log2FC calculation. +#' @param reference_group Group to be used as denominator in log2FC calculation. +#' @param fixed_effects A vector of \code{metadata} column names to use as fixed effects. +#' @param random_effects A vector of \code{metadata} column names to use as random effects. +#' @return \code{metadata} after matching the rows to the columns (samples) in \code{counts}. +.input_checks <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL, random_effects=NULL) { if (!is.matrix(counts) & !is.data.frame(counts) ) stop("Error. Input data is not a matrix or a data.frame.") if (ncol(counts)==0 | nrow(counts)==0) stop("Error. Count matrix does not have at least one column and row.") if(! dge_by %in% colnames(metadata) ) stop("Error. dge_by colum does not exist in metadata.") diff --git a/single_cell/dge/helper_functions/remove_low_expression_function.R b/single_cell/dge/helper_functions/remove_low_expression_function.R index 44b4777..6d4566f 100644 --- a/single_cell/dge/helper_functions/remove_low_expression_function.R +++ b/single_cell/dge/helper_functions/remove_low_expression_function.R @@ -7,7 +7,7 @@ #' @param dge_groups A vector of group labels (from \code{dge_by}) to consider. Only samples belonging to these groups are used to determine the minimum number of samples required to have CPM > 1. If NULL, all groups in \code{dge_by} are considered. #' @param min_frac Genes with counts per million (CPM) more than 1 in at least \code{min_frac} proportion of samples of the smaller group are retained, regardless of group membership. Default 0.6 (i.e., 60%). #' @return \code{counts} after removing lowly expressed genes. -remove_low_expression <- function(counts, metadata, dge_by, dge_groups=NULL, min_frac=0.6) { +.remove_low_expression <- function(counts, metadata, dge_by, dge_groups=NULL, min_frac=0.6) { # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least min_frac of min_gsize samples, regardless of the sample group identity. warning("Assuming that the order of samples (columns) in 'counts' matches the order of samples (rows) in 'metadata'.") if(is.null(dge_groups)) { From 524c0e812acd0a017e4cd7a461dea94a8ac10bab Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 13 Mar 2026 11:54:14 -0700 Subject: [PATCH 11/14] allow contrasts as input: input_checks_function.R --- single_cell/dge/deseq2_function.R | 18 ++++++------- .../helper_functions/input_checks_function.R | 25 +++++++++++-------- 2 files changed, 24 insertions(+), 19 deletions(-) diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index 8102f6b..c68b6dd 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -3,22 +3,22 @@ library(DESeq2) library(tidyverse) ### DGE analysis using DESeq2 # -#' Run DGE between 2 groups using DESeq2 -#' @param counts: feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts} -#' @param metadata: sample metadata, with rownames containing column names of \code{counts} -#' @param dge_by: name of \code{metadata} column to use DGE comparion. Must have exactly two levels -#' @param case_group: level to be used as numerator in log2FC calculation -#' @param reference_group: level to be used as denominator in log2FC calculation -#' @param fixed_effects: a vector of \code{metadata} column names to use as fixed effects +#' Run DGE between 2 groups using DESeq2. +#' @param counts A feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts}. +#' @param metadata Sample metadata, with rownames containing column names of \code{counts}. +#' @param dge_by Name of \code{metadata} column to use for DGE comparion. Must have exactly two levels. +#' @param case_group Group to be used as numerator in log2FC calculation. +#' @param reference_group Group to be used as denominator in log2FC calculation. +#' @param fixed_effects A vector of \code{metadata} column names to use as fixed effects. run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL) { # Commong input checks - metadata = input_checks(counts, metadata, dge_by, case_group, reference_group, fixed_effects) + .input_checks(counts, metadata, dge_by, case_group, reference_group, fixed_effects) # Method specific input checks # None for DESeq2 - counts = remove_low_expression(counts, metadata, dge_by, min_frac=0.6) + counts = .remove_low_expression(counts=counts, metadata=metadata, dge_by=dge_by, dge_groups=c(case_group, reference_group), min_frac=0.6) # Make formula string formula_str = paste("~", dge_by) diff --git a/single_cell/dge/helper_functions/input_checks_function.R b/single_cell/dge/helper_functions/input_checks_function.R index 1fffc57..ee9b432 100644 --- a/single_cell/dge/helper_functions/input_checks_function.R +++ b/single_cell/dge/helper_functions/input_checks_function.R @@ -1,28 +1,33 @@ ### Check inputs to DGE functions # #' Check inputs to DGE functions, match \code{metadata} rows to \code{counts} columns, and return \code{metadata}. -#' @param counts A feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts}. +#' @param counts A feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts}. #' @param metadata Sample metadata, with rownames containing column names of \code{counts}. #' @param dge_by Name of \code{metadata} column to use DGE comparion. Must have exactly two levels. #' @param case_group Group to be used as numerator in log2FC calculation. #' @param reference_group Group to be used as denominator in log2FC calculation. #' @param fixed_effects A vector of \code{metadata} column names to use as fixed effects. #' @param random_effects A vector of \code{metadata} column names to use as random effects. -#' @return \code{metadata} after matching the rows to the columns (samples) in \code{counts}. -.input_checks <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL, random_effects=NULL) { +.input_checks <- function(counts, metadata, dge_by, case_group=NULL, reference_group=NULL, contrasts=NULL, fixed_effects=NULL, random_effects=NULL) { if (!is.matrix(counts) & !is.data.frame(counts) ) stop("Error. Input data is not a matrix or a data.frame.") if (ncol(counts)==0 | nrow(counts)==0) stop("Error. Count matrix does not have at least one column and row.") - if(! dge_by %in% colnames(metadata) ) stop("Error. dge_by colum does not exist in metadata.") + if(! dge_by %in% colnames(metadata) ) stop("Error. dge_by column does not exist in metadata.") if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(metadata)) ) stop("Error. One or more fixed_effects do not exist in metadata.") if(! is.null(random_effects) & ! all(random_effects %in% colnames(metadata)) ) stop("Error. One or more random_effects do not exist in metadata.") - if(! all( colnames(counts) %in% rownames(metadata)) ) stop("Error. One or more samples in counts are not present in metadata rownames.") + if(! all( colnames(counts) == rownames(metadata)) ) stop("Error. Column names of counts and rownames of metadata do not match.") # Subset and order metadata according to the order of samples in counts - metadata = metadata[ match(colnames(counts), rownames(metadata)) , ] + #metadata = metadata[ match(colnames(counts), rownames(metadata)), ] - if( length(unique(as.vector(metadata[,dge_by]))) != 2) stop("Error. dge_by column of metadata does not have exactly two levels (unique values).") - if( ! case_group %in% metadata[,dge_by] ) stop("Case group does not exist in dge_by column of metadata.") - if( ! reference_group %in% metadata[,dge_by] ) stop("Reference group does not exist in dge_by column of metadata.") - return(metadata) + case_ref = FALSE + cntrst = FALSE + if( !is.null(case_group) & !is.null(reference_group) ) { + case_ref = TRUE + if( ! case_group %in% metadata[,dge_by] ) stop("Error. Case group does not exist in dge_by column of metadata.") + if( ! reference_group %in% metadata[,dge_by] ) stop("Error. Reference group does not exist in dge_by column of metadata.") + } + if(! is.null(contrasts) ) + cntrst = TRUE + if(case_ref & cntrst) warning("Both contrasts and case-reference groups are given. The contrasts will be ignored for the analysis.") } From 432cc327fba29d6b3bfc278ba3713ac3eda7bb6e Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 8 May 2026 09:27:03 -0700 Subject: [PATCH 12/14] delete helper functions --- single_cell/dge/deseq2_function.R | 2 +- .../helper_functions/input_checks_function.R | 33 ------------------- .../remove_low_expression_function.R | 26 --------------- 3 files changed, 1 insertion(+), 60 deletions(-) delete mode 100644 single_cell/dge/helper_functions/input_checks_function.R delete mode 100644 single_cell/dge/helper_functions/remove_low_expression_function.R diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index c68b6dd..009c6fb 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -13,7 +13,7 @@ library(tidyverse) run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL) { # Commong input checks - .input_checks(counts, metadata, dge_by, case_group, reference_group, fixed_effects) + contrasts = .input_checks(counts, metadata, dge_by, case_group, reference_group, fixed_effects) # Method specific input checks # None for DESeq2 diff --git a/single_cell/dge/helper_functions/input_checks_function.R b/single_cell/dge/helper_functions/input_checks_function.R deleted file mode 100644 index ee9b432..0000000 --- a/single_cell/dge/helper_functions/input_checks_function.R +++ /dev/null @@ -1,33 +0,0 @@ -### Check inputs to DGE functions -# -#' Check inputs to DGE functions, match \code{metadata} rows to \code{counts} columns, and return \code{metadata}. -#' @param counts A feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts}. -#' @param metadata Sample metadata, with rownames containing column names of \code{counts}. -#' @param dge_by Name of \code{metadata} column to use DGE comparion. Must have exactly two levels. -#' @param case_group Group to be used as numerator in log2FC calculation. -#' @param reference_group Group to be used as denominator in log2FC calculation. -#' @param fixed_effects A vector of \code{metadata} column names to use as fixed effects. -#' @param random_effects A vector of \code{metadata} column names to use as random effects. -.input_checks <- function(counts, metadata, dge_by, case_group=NULL, reference_group=NULL, contrasts=NULL, fixed_effects=NULL, random_effects=NULL) { - if (!is.matrix(counts) & !is.data.frame(counts) ) stop("Error. Input data is not a matrix or a data.frame.") - if (ncol(counts)==0 | nrow(counts)==0) stop("Error. Count matrix does not have at least one column and row.") - if(! dge_by %in% colnames(metadata) ) stop("Error. dge_by column does not exist in metadata.") - if(! is.null(fixed_effects) & ! all(fixed_effects %in% colnames(metadata)) ) stop("Error. One or more fixed_effects do not exist in metadata.") - if(! is.null(random_effects) & ! all(random_effects %in% colnames(metadata)) ) stop("Error. One or more random_effects do not exist in metadata.") - if(! all( colnames(counts) == rownames(metadata)) ) stop("Error. Column names of counts and rownames of metadata do not match.") - - # Subset and order metadata according to the order of samples in counts - #metadata = metadata[ match(colnames(counts), rownames(metadata)), ] - - case_ref = FALSE - cntrst = FALSE - if( !is.null(case_group) & !is.null(reference_group) ) { - case_ref = TRUE - if( ! case_group %in% metadata[,dge_by] ) stop("Error. Case group does not exist in dge_by column of metadata.") - if( ! reference_group %in% metadata[,dge_by] ) stop("Error. Reference group does not exist in dge_by column of metadata.") - } - if(! is.null(contrasts) ) - cntrst = TRUE - if(case_ref & cntrst) warning("Both contrasts and case-reference groups are given. The contrasts will be ignored for the analysis.") -} - diff --git a/single_cell/dge/helper_functions/remove_low_expression_function.R b/single_cell/dge/helper_functions/remove_low_expression_function.R deleted file mode 100644 index 6d4566f..0000000 --- a/single_cell/dge/helper_functions/remove_low_expression_function.R +++ /dev/null @@ -1,26 +0,0 @@ -### Remove lowly expressed genes from the count matrix. -# -#' Remove lowly expressed genes from the count matrix. The function first determines the size of the smallest group (denoted min_gsize), where groups are defined by \code{dge_by}. A gene is retaind only if it is expressed (CPM > 1) in at least a \code{min_frac} proportion of min_gsize samples overall (i.e., not within each group separately). -#' @param counts Feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts}. -#' @param metadata Sample metadata, with rownames containing column names of \code{counts}. -#' @param dge_by Name of \code{metadata} column to use DGE comparion. Must have exactly two levels. -#' @param dge_groups A vector of group labels (from \code{dge_by}) to consider. Only samples belonging to these groups are used to determine the minimum number of samples required to have CPM > 1. If NULL, all groups in \code{dge_by} are considered. -#' @param min_frac Genes with counts per million (CPM) more than 1 in at least \code{min_frac} proportion of samples of the smaller group are retained, regardless of group membership. Default 0.6 (i.e., 60%). -#' @return \code{counts} after removing lowly expressed genes. -.remove_low_expression <- function(counts, metadata, dge_by, dge_groups=NULL, min_frac=0.6) { - # Determine the number of sample in the smaller group (min_gsize) and keep features that are expressed with CPM > 1 in at least min_frac of min_gsize samples, regardless of the sample group identity. - warning("Assuming that the order of samples (columns) in 'counts' matches the order of samples (rows) in 'metadata'.") - if(is.null(dge_groups)) { - warning("dge_groups is not provided. Using all groups in dge_by to determine the minimum number of samples required to have CPM > 1.") - sample_idx = rep(TRUE, nrow(metadata)) - } else { - sample_idx = metadata[, dge_by] %in% dge_groups - } - min_gsize = min(table(metadata[ sample_idx, dge_by])) - feats = names(which(rowSums(cpm(counts[, sample_idx]) > 1) > (min_gsize * min_frac))) - cat(paste0("Number of final features: ", length(feats)), "\n") - counts = counts[feats,] - return(counts) -} - - From 0633ba78130edc4d9952170c0889b94947f53f4b Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 8 May 2026 09:49:42 -0700 Subject: [PATCH 13/14] fixed output df column names and added return_model functionality --- single_cell/dge/deseq2_function.R | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index 009c6fb..d274642 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -10,7 +10,8 @@ library(tidyverse) #' @param case_group Group to be used as numerator in log2FC calculation. #' @param reference_group Group to be used as denominator in log2FC calculation. #' @param fixed_effects A vector of \code{metadata} column names to use as fixed effects. -run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL) { +#' @param return_model A boolean indicating whether the function should return the model (TRUE) or the data.frame of differential expression analysis results. +run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fixed_effects=NULL, return_model = FALSE) { # Commong input checks contrasts = .input_checks(counts, metadata, dge_by, case_group, reference_group, fixed_effects) @@ -36,8 +37,13 @@ run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fi res = results(dds, contrast = c(dge_by, case_group, reference_group)) %>% as.data.frame() res = arrange(res, padj) - res = dplyr::rename(res, "logFC"="log2FoldChange", "AveExpr"="baseMean", "P.Value"="pvalue", "adj.P.Val"="padj") - return(res) + res = dplyr::rename(res, "log2FC"="log2FoldChange", "aveExpr"="baseMean", "pval"="pvalue", "padj"="padj") + + if(return_model) { + return(dds) + } else { + return(res) + } } From d7de13f98d9b598137952577781a32ec7081be66 Mon Sep 17 00:00:00 2001 From: Ravi Patel Date: Fri, 8 May 2026 12:04:41 -0700 Subject: [PATCH 14/14] run_deseq2_within_cells function --- single_cell/dge/deseq2_function.R | 42 +++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/single_cell/dge/deseq2_function.R b/single_cell/dge/deseq2_function.R index d274642..1957148 100644 --- a/single_cell/dge/deseq2_function.R +++ b/single_cell/dge/deseq2_function.R @@ -47,3 +47,45 @@ run_deseq2 <- function(counts, metadata, dge_by, case_group, reference_group, fi } + + + + + + +### DGE analysis using DESeq2 across cell-types +# +#' Run DGE between 2 groups using DESeq2 across cell-types. +#' @param counts A feature (row) x sample (column) raw count matrix. Only samples intended for DGE analysis should be included in \code{counts}. +#' @param metadata Sample metadata, with rownames containing column names of \code{counts}. +#' @param dge_by Name of \code{metadata} column to use for DGE comparion. Must have exactly two levels. +#' @param case_group Group to be used as numerator in log2FC calculation. +#' @param reference_group Group to be used as denominator in log2FC calculation. +#' @param fixed_effects A vector of \code{metadata} column names to use as fixed effects. +#' @param return_model A boolean indicating whether the function should return the model (TRUE) or the data.frame of differential expression analysis results. +run_deseq2_within_cells <- function(counts, metadata, cell_by, cell_targets, dge_by, case_group, reference_group, fixed_effects=NULL, return_model = FALSE) { + + # Check inputs + cell_targets = .input_checks_within_cells(metadata, cell_by, cell_target, dge_by, case_group, reference_group) + + dge_res = list() + for(cell in cell_targets) { + # Assuming the rownames of metadata are in counts columns + cell_counts = cell_counts[ , rownames(metadata)[ metadata[,cell_by] == cell ] ] + cell_metadata = metadata[ metadata[,cell_by] == cell, ] + dge_res[[cell]] = run_deseq2(counts=cell_counts, + metadata=cell_metadata, + dge_by=dge_by, + case_group=case_group, + reference_group=reference_group, + fixed_effects=fixed_effects, + return_model = return_model) + if(! return_model) + dge_res[[cell]][, cell_by] = cell + } + + if(! return_model) + dge_res <- bind_rows(dge_res) + return(dge_res) +} +