Notebook

Differential gene expression with DESeq2¶

In this notebook, we'll perform differential gene expression tests using the DESeq2 package on our aggregated pseudobulk counts for each subject and cell type.

Load packages¶

In [1]:

quiet_library <- function(...) {suppressPackageStartupMessages(library(...))}
quiet_library(DESeq2)
quiet_library(dplyr)
quiet_library(hise)
quiet_library(purrr)
quiet_library(furrr)

In [2]:

plan(multicore, workers = 12)

In [3]:

if(!dir.exists("output")) {
    dir.create("output")
}
in_files <- c()
out_files <- c()

Helper functions¶

This function prepares data for analysis by setting factors, converting the matrix of data to a non-sparse, integer matrix, and assembling a DESeqDataSet with a given factor design.

In [4]:

prepare_deseq_data <- function(mat, meta, design) {
    meta$cohort.cohortGuid <- factor(meta$cohort.cohortGuid, levels = c("BR1", "BR2"))
    meta$subject.biologicalSex <- factor(meta$subject.biologicalSex, levels = c("Female", "Male"))
    meta$subject.cmv <- factor(meta$subject.cmv, levels = c("Negative", "Positive"))
    
    mat <- mat[,meta$barcodes]
    
    mat <- as.matrix(mat)
    mode(mat) <- "integer"
    
    DESeqDataSetFromMatrix(
        mat, colData = meta,
        design = design
    )
}

This function retrieves results from a DESeqDataSet after modeling and builds a data.frame with these results, cell types, and data grouping metadata.

In [5]:

format_deseq_results <- function(dds, cell_type, group_by, fg, bg) {
        meta <- as.data.frame(colData(dds))
        
        res <- results(dds, contrast = c(group_by, fg, bg))
        res <- data.frame(res)

        
        res <- res %>%
          arrange(padj) %>%
          mutate(direction = ifelse(log2FoldChange > 0,
                                    paste0("HigherIn", fg), 
                                    paste0("HigherIn", bg))) %>%
          mutate(gene = rownames(.),
                 cell_type = cell_type,
                 fg = fg, n_fg = sum(res[[group_by]] == fg),
                 bg = bg, n_bg = sum(res[[group_by]] == bg)) %>%
          select(cell_type, fg, n_fg, bg, n_bg, gene, log2FoldChange, padj, direction,
                 baseMean, lfcSE, stat, pvalue)

        res
        
    }

Single-factor comparisons¶

Here, we'll identify differences in comparisons between cohorts, sexes, and CMV status across all subjects while controlling for the other two factors.

Retrieve pseudobulk data from HISE¶

In [6]:

all_uuid <- "edf436fb-7063-4655-9055-6ad431fe7159"

In [7]:

res <- cacheFiles(list(all_uuid))
pb_file <- list.files(paste0("cache/", all_uuid), pattern = ".rds", full.names = TRUE)

In [8]:

pb_data <- readRDS(pb_file)

In [9]:

in_files = c(in_files, all_uuid)

YA vs OA; control for sex and CMV status¶

To test for cohort differences, we'll specify the design of the differential expression tests using a function:
design = ~ cohort.cohortGuid + subject.biologicalSex + subject.cmv

In [10]:

dds_data <- map2(
    pb_data$mats, pb_data$meta,
    prepare_deseq_data,
    design = ~ cohort.cohortGuid + subject.biologicalSex + subject.cmv
)

Run DESeq2¶

In [11]:

deseq_res <- future_map(
    dds_data,
    DESeq,
    parallel = FALSE, # parallelization handled by future_map
    quiet = TRUE,
    .options = furrr_options(seed = 3030L)
)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

Structure results per cell type¶

In [12]:

cohort_res <- map2_dfr(
    deseq_res, names(deseq_res),
    format_deseq_results,
    group_by = "cohort.cohortGuid",
    fg = "BR1", bg = "BR2"
)

Save results¶

In [13]:

out_group <- "cohort"
out_csv <- paste0(
    "output/diha_pbmc_", out_group, "_pseudobulk_deseq2_", Sys.Date(), ".csv"
)
write.csv(
    cohort_res, out_csv,
    quote = FALSE, row.names = FALSE
)
out_files <- c(out_files, out_csv)

CMV-Negative vs CMV-Positive; control for Cohort and sex¶

To test for CMV differences, we'll specify the design of the differential expression tests using a function:
design = ~ subject.cmv + cohort.cohortGuid + subject.biologicalSex

In [14]:

dds_data <- map2(
    pb_data$mats, pb_data$meta,
    prepare_deseq_data,
    design = ~ subject.cmv + cohort.cohortGuid + subject.biologicalSex
)

Run DESeq2¶

In [15]:

deseq_res <- future_map(
    dds_data,
    DESeq,
    parallel = FALSE, # parallelization handled by future_map
    quiet = TRUE,
    .options = furrr_options(seed = 3030L)
)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

Structure results per cell type¶

In [16]:

cmv_res <- map2_dfr(
    deseq_res, names(deseq_res),
    format_deseq_results,
    group_by = "subject.cmv",
    fg = "Positive", bg = "Negative"
)

Save results¶

In [17]:

out_group <- "cmv"
out_csv <- paste0(
    "output/diha_pbmc_", out_group, "_pseudobulk_deseq2_", Sys.Date(), ".csv"
)
write.csv(
    cmv_res, out_csv,
    quote = FALSE, row.names = FALSE
)
out_files <- c(out_files, out_csv)

Female vs Male; control for Cohort and CMV status¶

To test for CMV differences, we'll specify the design of the differential expression tests using a function:
design = ~ subject.biologicalSex + cohort.cohortGuid + subject.cmv

In [18]:

dds_data <- map2(
    pb_data$mats, pb_data$meta,
    prepare_deseq_data,
    design = ~ subject.biologicalSex + cohort.cohortGuid + subject.cmv
)

Run DESeq2¶

In [19]:

deseq_res <- future_map(
    dds_data,
    DESeq,
    parallel = FALSE, # parallelization handled by future_map
    quiet = TRUE,
    .options = furrr_options(seed = 3030L)
)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

Structure results per cell type¶

In [20]:

sex_res <- map2_dfr(
    deseq_res, names(deseq_res),
    format_deseq_results,
    group_by = "subject.biologicalSex",
    fg = "Female", bg = "Male"
)

Save results¶

In [21]:

out_group <- "sex"
out_csv <- paste0(
    "output/diha_pbmc_", out_group, "_pseudobulk_deseq2_", Sys.Date(), ".csv"
)
write.csv(
    sex_res, out_csv,
    quote = FALSE, row.names = FALSE
)
out_files <- c(out_files, out_csv)

Within-factor comparisons¶

Here, we'll identify differences between the cohorts within sex and CMV status.

To ensure gene selection was based on the data available within these groups, we previously generated pseudobulk datasets for the BR1 and BR2 cohorts, for CMV+ and CMV- subjects, and for Female and Male subjects separately.

Here, we're most interested in the effect of age - that is, comparisons between BR1 and BR2 - within the other two factors, CMV+ and CMV- or Male and Female.

We could perform every possible subsetting, but for now we'll stick with these 4 additional comparisons - BR1 vs BR2 within CMV+, within CMV-, within Female, and within Male subjects.

Cohort within CMV positive¶

Retrieve pseudobulk data from HISE¶

In [22]:

cmv_pos_uuid <- "c8e65e37-d11d-4baf-b450-6eda124ef102"

In [23]:

res <- cacheFiles(list(cmv_pos_uuid))
pb_file <- list.files(paste0("cache/", cmv_pos_uuid), pattern = ".rds", full.names = TRUE)

In [24]:

pb_data <- readRDS(pb_file)

In [25]:

in_files = c(in_files, cmv_pos_uuid)

Prep data for analysis¶

Here, we're comparing cohort and controlling for sex within our CMV-positive subjects, so our design formula is:

~ cohort.cohortGuid + subject.biologicalSex

In [26]:

table(pb_data$meta[[2]]$subject.cmv)

Positive 
      37

In [27]:

ncol(pb_data$mats[[2]])

In [28]:

dds_data <- map2(
    pb_data$mats, pb_data$meta,
    prepare_deseq_data,
    design = ~ cohort.cohortGuid + subject.biologicalSex
)

Run DESeq2 and extract results from tests¶

In [29]:

deseq_res <- future_map(
    dds_data,
    DESeq,
    parallel = FALSE, # parallelization handled by future_map
    quiet = TRUE,
    .options = furrr_options(seed = 3030L)
)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

Structure results per cell type¶

In [30]:

cmvpos_cohort_res <- map2_dfr(
    deseq_res, names(deseq_res),
    format_deseq_results,
    group_by = "cohort.cohortGuid",
    fg = "BR1", bg = "BR2"
)

Save results¶

In [31]:

out_group <- "cmvpos_cohort"
out_csv <- paste0(
    "output/diha_pbmc_", out_group, "_pseudobulk_deseq2_", Sys.Date(), ".csv"
)
write.csv(
    cmvpos_cohort_res, out_csv,
    quote = FALSE, row.names = FALSE
)
out_files <- c(out_files, out_csv)

Cohort within CMV negative¶

Retrieve pseudobulk data from HISE¶

In [32]:

cmv_neg_uuid <- "42de2197-4ce6-400e-a0e7-212f2b25e896"

In [33]:

res <- cacheFiles(list(cmv_neg_uuid))
pb_file <- list.files(paste0("cache/", cmv_neg_uuid), pattern = ".rds", full.names = TRUE)

In [34]:

pb_data <- readRDS(pb_file)

In [35]:

in_files <- c(in_files, cmv_neg_uuid)

Prep data for analysis¶

Here, we're comparing cohort and controlling for sex within our CMV-negative subjects, so our design formula is:

~ cohort.cohortGuid + subject.biologicalSex

In [36]:

dds_data <- map2(
    pb_data$mats, pb_data$meta,
    prepare_deseq_data,
    design = ~ cohort.cohortGuid + subject.biologicalSex
)

Run DESeq2 and extract results from tests¶

In [37]:

deseq_res <- future_map(
    dds_data,
    DESeq,
    parallel = FALSE, # parallelization handled by future_map
    quiet = TRUE,
    .options = furrr_options(seed = 3030L)
)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

Structure results per cell type¶

In [38]:

cmvneg_cohort_res <- map2_dfr(
    deseq_res, names(deseq_res),
    format_deseq_results,
    group_by = "cohort.cohortGuid",
    fg = "BR1", bg = "BR2"
)

Save results¶

In [39]:

out_group <- "cmvneg_cohort"
out_csv <- paste0(
    "output/diha_pbmc_", out_group, "_pseudobulk_deseq2_", Sys.Date(), ".csv"
)
write.csv(
    cmvneg_cohort_res, out_csv,
    quote = FALSE, row.names = FALSE
)
out_files <- c(out_files, out_csv)

Cohort within Female subjects¶

Retrieve pseudobulk data from HISE¶

In [40]:

female_uuid <- "56f09312-eb3f-40b7-b840-368d1ddbc22b"

In [41]:

res <- cacheFiles(list(female_uuid))
pb_file <- list.files(paste0("cache/", female_uuid), pattern = ".rds", full.names = TRUE)

In [42]:

pb_data <- readRDS(pb_file)

In [43]:

in_files = c(in_files, female_uuid)

Prep data for analysis¶

Here, we're comparing cohort and controlling for CMV status within our Female subjects, so our design formula is:

~ cohort.cohortGuid + subject.cmv

In [44]:

dds_data <- map2(
    pb_data$mats, pb_data$meta,
    prepare_deseq_data,
    design = ~ cohort.cohortGuid + subject.cmv
)

Run DESeq2 and extract results from tests¶

In [45]:

deseq_res <- future_map(
    dds_data,
    DESeq,
    parallel = FALSE, # parallelization handled by future_map
    quiet = TRUE,
    .options = furrr_options(seed = 3030L)
)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

Structure results per cell type¶

In [46]:

female_cohort_res <- map2_dfr(
    deseq_res, names(deseq_res),
    format_deseq_results,
    group_by = "cohort.cohortGuid",
    fg = "BR1", bg = "BR2"
)

Save results¶

In [47]:

out_group <- "female_cohort"
out_csv <- paste0(
    "output/diha_pbmc_", out_group, "_pseudobulk_deseq2_", Sys.Date(), ".csv"
)
write.csv(
    female_cohort_res, out_csv,
    quote = FALSE, row.names = FALSE
)
out_files <- c(out_files, out_csv)

Cohort within Male subjects¶

Retrieve pseudobulk data from HISE¶

In [48]:

male_uuid <- "4bf2bab1-b3c8-456b-aa90-566cd7b03577"

In [49]:

res <- cacheFiles(list(male_uuid))
pb_file <- list.files(paste0("cache/", male_uuid), pattern = ".rds", full.names = TRUE)

In [50]:

pb_data <- readRDS(pb_file)

In [51]:

in_files <- c(in_files, male_uuid)

Prep data for analysis¶

Here, we're comparing cohort and controlling for CMV status within our Female subjects, so our design formula is:

~ cohort.cohortGuid + subject.cmv

In [52]:

dds_data <- map2(
    pb_data$mats, pb_data$meta,
    prepare_deseq_data,
    design = ~ cohort.cohortGuid + subject.cmv
)

Run DESeq2 and extract results from tests¶

In [53]:

deseq_res <- future_map(
    dds_data,
    DESeq,
    parallel = FALSE, # parallelization handled by future_map
    quiet = TRUE,
    .options = furrr_options(seed = 3030L)
)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

Structure results per cell type¶

In [54]:

male_cohort_res <- map2_dfr(
    deseq_res, names(deseq_res),
    format_deseq_results,
    group_by = "cohort.cohortGuid",
    fg = "BR1", bg = "BR2"
)

Save results¶

In [55]:

out_group <- "male_cohort"
out_csv <- paste0(
    "output/diha_pbmc_", out_group, "_pseudobulk_deseq2_", Sys.Date(), ".csv"
)
write.csv(
    male_cohort_res, out_csv,
    quote = FALSE, row.names = FALSE
)
out_files <- c(out_files, out_csv)

Upload results to HISE¶

In [56]:

study_space_uuid <- "64097865-486d-43b3-8f94-74994e0a72e0"
title <- paste("PBMC Ref. Seurat Pseudobulk DESeq2 Res.", Sys.Date())

In [57]:

in_list <- as.list(in_files)
in_list

'edf436fb-7063-4655-9055-6ad431fe7159'
'c8e65e37-d11d-4baf-b450-6eda124ef102'
'42de2197-4ce6-400e-a0e7-212f2b25e896'
'56f09312-eb3f-40b7-b840-368d1ddbc22b'
'4bf2bab1-b3c8-456b-aa90-566cd7b03577'

In [58]:

out_list <- as.list(out_files)
out_list

'output/diha_pbmc_cohort_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_cmv_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_sex_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_cmvpos_cohort_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_cmvneg_cohort_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_female_cohort_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_male_cohort_pseudobulk_deseq2_2024-04-05.csv'

In [59]:

uploadFiles(
    files = out_list,
    studySpaceId = study_space_uuid,
    title = title,
    inputFileIds = in_list,
    store = "project",
    destination = "pseudobulk_deseq2",
    doPrompt = FALSE
)

[1] "Cannot determine the current notebook."
[1] "1) /home/jupyter/scRNA-Reference-IH-A/08-Differential_expression/36-R_pseudobulk_deseq2.ipynb"
[1] "2) /home/jupyter/scRNA-Reference-IH-A/08-Differential_expression/36a-Python_pseudobulk_deseq2_vis_format.ipynb"
[1] "3) /home/jupyter/scRNA-Reference-IH-A/07-Pseudobulk/35-R_pseudobulk_seurat.ipynb"

$files

'output/diha_pbmc_cohort_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_cmv_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_sex_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_cmvpos_cohort_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_cmvneg_cohort_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_female_cohort_pseudobulk_deseq2_2024-04-05.csv'
'output/diha_pbmc_male_cohort_pseudobulk_deseq2_2024-04-05.csv'

$traceId

'4f461fd6-3cff-468a-88e5-c1621734eaba'

In [60]:

sessionInfo()

R version 4.3.2 (2023-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS/LAPACK: /opt/conda/lib/libopenblasp-r0.3.25.so;  LAPACK version 3.11.0

Random number generation:
 RNG:     L'Ecuyer-CMRG 
 Normal:  Inversion 
 Sample:  Rejection 
 
locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] furrr_0.3.1                 future_1.33.1              
 [3] purrr_1.0.2                 hise_2.16.0                
 [5] dplyr_1.1.4                 DESeq2_1.42.0              
 [7] SummarizedExperiment_1.32.0 Biobase_2.62.0             
 [9] MatrixGenerics_1.14.0       matrixStats_1.2.0          
[11] GenomicRanges_1.54.1        GenomeInfoDb_1.38.5        
[13] IRanges_2.36.0              S4Vectors_0.40.2           
[15] BiocGenerics_0.48.1        

loaded via a namespace (and not attached):
 [1] gtable_0.3.4            ggplot2_3.4.4           lattice_0.22-5         
 [4] vctrs_0.6.5             tools_4.3.2             bitops_1.0-7           
 [7] generics_0.1.3          curl_5.1.0              parallel_4.3.2         
[10] tibble_3.2.1            fansi_1.0.6             pkgconfig_2.0.3        
[13] Matrix_1.6-4            assertthat_0.2.1        uuid_1.2-0             
[16] lifecycle_1.0.4         GenomeInfoDbData_1.2.11 stringr_1.5.1          
[19] compiler_4.3.2          munsell_0.5.0           repr_1.1.6.9000        
[22] codetools_0.2-19        htmltools_0.5.7         RCurl_1.98-1.14        
[25] pillar_1.9.0            crayon_1.5.2            BiocParallel_1.36.0    
[28] DelayedArray_0.28.0     abind_1.4-5             parallelly_1.36.0      
[31] tidyselect_1.2.0        locfit_1.5-9.8          digest_0.6.34          
[34] stringi_1.8.3           listenv_0.9.0           fastmap_1.1.1          
[37] grid_4.3.2              colorspace_2.1-0        cli_3.6.2              
[40] SparseArray_1.2.3       magrittr_2.0.3          S4Arrays_1.2.0         
[43] base64enc_0.1-3         utf8_1.2.4              IRdisplay_1.1          
[46] withr_3.0.0             scales_1.3.0            IRkernel_1.3.2         
[49] httr_1.4.7              XVector_0.42.0          globals_0.16.2         
[52] pbdZMQ_0.3-10           evaluate_0.23           rlang_1.1.3            
[55] Rcpp_1.0.12             glue_1.7.0              jsonlite_1.8.8         
[58] R6_2.5.1                zlibbioc_1.48.0

In [ ]: