Revisiting a macrogen experiment.

• TODO: Box, create ‘plasma’, print out plots, name xlsx appropriately. Check out olps further.

1 Introduction

2 Gene annotations

hg_annot <- load_biomart_annotations()

## The biomart annotations file already exists, loading from it.

hg_df <- hg_annot[["gene_annotations"]]

3 Sample sheet annotation

We have an older revision of the sample sheet for this dataset. I added some samples from Dr. Mosser in order to compare against M1/M2 activation states. These extra samples are not likely to be the most appropriate because they are not U937 samples.

hg38_se <- create_se("sample_sheets/macrogen_samples.xlsx",
                     file_column = "hisat_hg38", gene_info = hg_df)

## Reading the sample metadata.

## Checking the state of the condition column.

## Warning in extract_metadata(metadata, id_column = id_column, condition_column =
## condition_column, : There were NA values in the condition column, setting them
## to 'undefined'.

## Checking the state of the batch column.

## Warning in extract_metadata(metadata, id_column = id_column, condition_column =
## condition_column, : There were NA values in the condition column, setting them
## to 'undefined'.

## Checking the condition factor.

## The sample definitions comprises: 54 rows(samples) and 21 columns(metadata fields).

## Warning in create_se("sample_sheets/macrogen_samples.xlsx", file_column =
## "hisat_hg38", : Some samples were removed when cross referencing the samples
## against the count data.

## Matched 21405 annotations and counts.

## Some annotations were lost in merging, setting them to 'undefined'.

## The final summarized experiment has 21481 rows and 21 columns.

hg38_sampletype <- set_conditions(hg38_se, fact = "sampletype")

## The numbers of samples by condition are:

## 
## Asymptomatic      Chronic      control      Healthy 
##            5            5            3            2

4 Plots

hg38_norm <- normalize(hg38_sampletype, convert = "cpm", norm = "quant",
                       transform = "log2", filter = TRUE)

## Removing 9903 low-count genes (11578 remaining).

## transform_counts: Found 26 values equal to 0, adding 1 to the matrix.

plot_quantreads(hg38_sampletype)

## Library sizes of 15 samples, 
## ranging from 1,872,219 to 7,805,747.

plot_nonzero(hg38_sampletype, y_intercept = 0.65)

## The following samples have less than 13962.65 genes.

##  [1] "m1"      "m2"      "a_20179" "c_10036" "a_20187" "c_10046" "a_20132"
##  [8] "c_10063" "a_20133" "c_10093" "su1160"  "a_20134" "c_10073"

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the hpgltools package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## A non-zero genes plot of 15 samples.
## These samples have an average 3.382 CPM coverage and 13550 genes observed, ranging from 12975 to
## 14728.

plot_corheat(hg38_norm)

## A heatmap of pairwise sample correlations ranging from: 
## 0.932974671519701 to 0.990409066934894.

plot_pca(hg38_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Asymptomatic, Chronic, control, Healthy
## Shapes are defined by undefined.

hg38_nb <- normalize(hg38_sampletype, convert = "cpm", batch = "svaseq",
                     filter = TRUE, transform = "log2")

## Removing 9903 low-count genes (11578 remaining).
## transform_counts: Found 169 values less than 0.
## transform_counts: Found 169 values equal to 0, adding 1 to the matrix.

plot_pca(hg38_nb)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Asymptomatic, Chronic, control, Healthy
## Shapes are defined by undefined.

5 DE

keepers <- list(
  "chr_asy" = c("Chronic", "Asymptomatic"),
  "chr_hea" = c("Chronic", "Healthy"),
  "chr_con" = c("Chronic", "control"),
  "asy_hea" = c("Asymptomatic", "Healthy"),
  "asy_con" = c("Asymptomatic", "control"),
  "hea_con" = c("Healthy", "control"))

hg38_de <- all_pairwise(hg38_sampletype, keepers = keepers, model_fstring = "~ 0 + condition",
                        model_svs = "svaseq", filter = TRUE)

## Asymptomatic      Chronic      control      Healthy 
##            5            5            3            2

## Removing 9903 low-count genes (11578 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 4687 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## conditions
## Asymptomatic      Chronic      control      Healthy 
##            5            5            3            2

## conditions
## Asymptomatic      Chronic      control      Healthy 
##            5            5            3            2

## conditions
## Asymptomatic      Chronic      control      Healthy 
##            5            5            3            2

hg38_de

## A pairwise differential expression with results from: basic, deseq, ebseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 6 comparisons.

hg38_tables <- combine_de_tables(hg38_de, excel = "excel/hg38_tables.xlsx")

## Deleting the file excel/hg38_tables.xlsx before writing the tables.

## Looking for subscript invalid names, start of extract_keepers.

## Looking for subscript invalid names, end of extract_keepers.

hg38_tables

## A set of combined differential expression results.

##                     table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 Chronic_vs_Asymptomatic          31            41          43            73
## 2      Chronic_vs_Healthy          59            43          98            77
## 3      Chronic_vs_control         680           458         748           551
## 4 Asymptomatic_vs_Healthy         112            83         183           127
## 5 Asymptomatic_vs_control         759           492         837           575
## 6      Healthy_vs_control         427           311         488           421
##   limma_sigup limma_sigdown
## 1          46            28
## 2          34            46
## 3         611           554
## 4          92           113
## 5         646           626
## 6         414           392

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the UpSetR package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the UpSetR package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Plot describing unique/shared genes in a differential expression table.

hg38_sig <- extract_significant_genes(hg38_tables, excel = "excel/hg38_sig.xlsx",
                                      according_to = "deseq")

## Deleting the file excel/hg38_sig.xlsx before writing the tables.

hg38_sig

## A set of genes deemed significant according to deseq.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                         deseq_up deseq_down
## Chronic_vs_Asymptomatic       31         41
## Chronic_vs_Healthy            59         43
## Chronic_vs_control           680        458
## Asymptomatic_vs_Healthy      112         83
## Asymptomatic_vs_control      759        492
## Healthy_vs_control           427        311

6 Repeat without the controls

no_controls <- subset_se(hg38_sampletype, subset = "condition!='control'")

no_controls_norm <- normalize(no_controls, transform = "log2", convert = "cpm",
                              norm = "quant", filter = TRUE)

## Removing 10032 low-count genes (11449 remaining).

## transform_counts: Found 6 values equal to 0, adding 1 to the matrix.

plot_pca(no_controls_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Asymptomatic, Chronic, Healthy
## Shapes are defined by undefined.

no_control_keepers <- list(
  "chr_asy" = c("Chronic", "Asymptomatic"),
  "chr_hea" = c("Chronic", "Healthy"),
  "asy_hea" = c("Asymptomatic", "Healthy"))
no_control_de <- all_pairwise(no_controls, keepers = no_control_keepers,
                              model_fstring = "~ 0 + condition", model_svs = "svaseq",
                              filter = TRUE)

## Asymptomatic      Chronic      Healthy 
##            5            5            2

## Removing 10032 low-count genes (11449 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 2623 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## conditions
## Asymptomatic      Chronic      Healthy 
##            5            5            2

## conditions
## Asymptomatic      Chronic      Healthy 
##            5            5            2

## conditions
## Asymptomatic      Chronic      Healthy 
##            5            5            2

no_control_table <- combine_de_tables(no_control_de, excel = "excel/hg38_nocontrols_tables.xlsx")

## Deleting the file excel/hg38_nocontrols_tables.xlsx before writing the tables.

## Looking for subscript invalid names, start of extract_keepers.

## Looking for subscript invalid names, end of extract_keepers.

no_control_sig <- extract_significant_genes(
  no_control_table, according_to = "deseq", excel = "excel/hg38_nocontrol_sig.xlsx")

## Deleting the file excel/hg38_nocontrol_sig.xlsx before writing the tables.

pp(file = "images/no_control_sig_barplot.png", image = no_control_sig[["sig_bar_plots"]][["deseq"]])

7 Repeat GSE(A) analyses

all_gp <- all_gprofiler(hg38_sig)
all_cp <- all_cprofiler(hg38_sig, hg38_tables)

## Error in simple_clusterprofiler(sig_genes = structure(list(ensembl_transcript_id = c("ENST00000492807",  : 
##   unused argument (internal = FALSE)
## Error in simple_clusterprofiler(sig_genes = structure(list(ensembl_transcript_id = c("ENST00000256104",  : 
##   unused argument (internal = FALSE)

## Error in `simple_cl[["kegg_universe"]]`:
## ! subscript out of bounds

## cp_test <- simple_clusterprofiler(ups, de_table = table, orgdb = pombe_orgdb,
##                                  orgdb_to = "GID", orgdb_from = "GID")

pander::pander(sessionInfo())
message(paste0("This is hpgltools commit: ", get_git_commit()))
message(paste0("Saving to ", savefile))
tmp <- sm(saveme(filename = savefile))

tmp <- loadme(filename = savefile)