Run DESeq2 differential expression analysis on RNA-seq count data

Run the DESeq2 differential gene expression analysis workflow from input file counts matrix and conditional labels.

The workflow will read in standard output count matrix from featureCounts, run DGE analysis between the two conditions provided in the formula string and conditions label file.

The counts matrix should contain a column with either the gene symbols (i.e. HUGO) or Ensembl IDs (i.e. ENSG###). The final output gene set will be mapped to Ensembl v108. See gene_body_hg38 for more info.

The formula_string should be a character string of the desired comparison in standard y ~ x format. Both x and y should be present in the conditions label file. This formula also handle covariates which should be included before x. For example, y ~ age - BMI + x. The covariates should also be present in the conditions label file as separate columns.

The counts will be filtered to remove genes with little to no expression across all samples using a minimum threshold for mean expression of a gene. This filter can be deactivated by setting min_transcripts = 0. The results of this applied filter will be captured in QC plots.

The final output will be a list that includes: the full DESeq2 object, a quilt plot of QC diagnostics after filtering, the DGE results as a tidy data.table, and a transcripts per million (TPM) counts data.table.

Usage

get_deseq2_diff_expr(
  counts_file_path,
  condition_file_path,
  gene_universe = gene_body_hg38,
  gene_symbol_column = NULL,
  ensembl_id_column = NULL,
  formula_string,
  min_transcripts = 10,
  verbose = T
)

Arguments

counts_file_path

Path to the counts matrix file

condition_file_path

Path to the conditions file

gene_universe

GenomicRanges object with genes to include for DGE analysis, default: gene_body_hg38

gene_symbol_column

The column name that contains HUGO gene symbols

ensembl_id_column

The column name that contains Ensembl gene IDs

formula_string

A character string of the formula of desired comparison

min_transcripts

The minimum threshold for mean transcript counts per gene, default: 10

verbose

Output CLI during workflow, default: TRUE

Value

List object with Full DESeq2 object, diagnostic plots, DGE log fold change results data.table, per gene TPM results data.table

Examples

# dge_demo <- get_deseq2_diff_expr(
# counts_file_path = fs::path_package("extdata", "hmcl_counts.hg38.txt.gz", package = "devgru"),
# condition_file_path = fs::path_package("extdata", "hmcl_conditions.hg38.txt", package = "devgru"),
# gene_universe = gene_body_hg38,
# gene_symbol_column = "#GENE",
# ensembl_id_column = NULL,
# formula_string = "U266 ~ RPMI8226",
# min_transcripts = 10)

# Look at the output
# Full DESeq2 object
# dge_demo$deseq2_obj

# Diagnostic plots after simple filtering
# dge_demo$diagnostic_plots

# DGE results as DT
# dge_demo$deseq2_results_dt

# Per gene TPM as DT
# dge_demo$deseq2_tpm_dt