TMRC3 Comprehensive Data Analysis: 202104

1 Introduction

This document is intended to provide an overview of TMRC3 samples which have been sequenced. It includes some plots and analyses showing the relationships among the samples as well as some differential analyses when possible.

2 Annotation

We take the annotation data from ensembl’s biomart instance. The genome which was used to map the data was hg38 revision 100. My default when using biomart is to load the data from 1 year before the current date.

hs_annot <- sm(load_biomart_annotations(year = "2020"))
hs_annot <- hs_annot[["annotation"]]
hs_annot[["transcript"]] <- paste0(rownames(hs_annot), ".", hs_annot[["version"]])
rownames(hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique = TRUE)
tx_gene_map <- hs_annot[, c("transcript", "ensembl_gene_id")]

summary(hs_annot)

##  ensembl_transcript_id ensembl_gene_id       version     transcript_version
##  Length:227921         Length:227921      Min.   : 1.0   Min.   : 1.00     
##  Class :character      Class :character   1st Qu.: 6.0   1st Qu.: 1.00     
##  Mode  :character      Mode  :character   Median :12.0   Median : 1.00     
##                                           Mean   :10.7   Mean   : 3.08     
##                                           3rd Qu.:16.0   3rd Qu.: 5.00     
##                                           Max.   :29.0   Max.   :17.00     
##                                                                            
##  hgnc_symbol        description        gene_biotype         cds_length    
##  Length:227921      Length:227921      Length:227921      Min.   :     3  
##  Class :character   Class :character   Class :character   1st Qu.:   357  
##  Mode  :character   Mode  :character   Mode  :character   Median :   694  
##                                                           Mean   :  1139  
##                                                           3rd Qu.:  1446  
##                                                           Max.   :107976  
##                                                           NA's   :127343  
##  chromosome_name       strand          start_position      end_position     
##  Length:227921      Length:227921      Min.   :5.77e+02   Min.   :6.47e+02  
##  Class :character   Class :character   1st Qu.:3.11e+07   1st Qu.:3.12e+07  
##  Mode  :character   Mode  :character   Median :6.04e+07   Median :6.06e+07  
##                                        Mean   :7.41e+07   Mean   :7.42e+07  
##                                        3rd Qu.:1.09e+08   3rd Qu.:1.09e+08  
##                                        Max.   :2.49e+08   Max.   :2.49e+08  
##                                                                             
##   transcript       
##  Length:227921     
##  Class :character  
##  Mode  :character  
##                    
##                    
##                    
## 

hs_go <- sm(load_biomart_go()[["go"]])
hs_length <- hs_annot[, c("ensembl_gene_id", "cds_length")]
colnames(hs_length) <- c("ID", "length")

3 Sample Estimation

3.1 Generate expressionsets

The sample sheet is copied from our shared online sheet and updated with each release of sequencing data.

samplesheet <- "sample_sheets/tmrc3_samples_202103.xlsx"

3.1.1 Hisat2 expressionsets

The first thing to note is the large range in coverage. There are multiple samples with coverage which is too low to use. These will be removed shortly.

In the following block I immediately exclude any non-coding reads as well.

## Create the expressionset and immediately pass it to a filter
## removing the non protein coding genes.
hs_expt <- sm(create_expt(samplesheet,
                          file_column = "hg38100hisatfile",
                          savefile = glue::glue("rda/hs_expt_all-v{ver}.rda"),
                          gene_info = hs_annot)) %>%
  exclude_genes_expt(column = "gene_biotype", method = "keep", pattern = "protein_coding")

## Before removal, there were 21481 genes, now there are 19941.

## There are 11 samples which kept less than 90 percent counts.

## TMRC30015 TMRC30017 TMRC30019 TMRC30044 TMRC30045 TMRC30097 TMRC30075 TMRC30087 
##     79.24     85.72     89.75     80.34     73.33     89.90     86.97     83.63 
## TMRC30101 TMRC30104 TMRC30073 
##     88.41     80.29     89.26

3.1.1.1 Initial metrics

Once the data was loaded, there are a couple of metrics which may be plotted immediately.

nonzero <- plot_nonzero(hs_expt)
nonzero$plot

## Warning: ggrepel: 80 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

3.2 Minimum coverage sample filtering

I arbitrarily chose 11,000 non-zero genes as a minimum. We may want this to be higher.

hs_valid <- subset_expt(hs_expt, nonzero = 11000)

## The samples (and read coverage) removed when filtering 11000 non-zero genes are:

## TMRC30010 TMRC30050 TMRC30052 
##     52471    808149   3087347

## There were 108, now there are 105 samples.

## valid_write <- write_expt(hs_valid, excel = glue("excel/hs_valid-v{ver}.xlsx"))

4 Project Aims

The project seeks to determine the relationship of the innate immune response and inflammatory signaling to the clinical outcome of antileishmanial drug treatment. We will test the hypothesis that the profile of innate immune cell activation and their dynamics through the course of treatment differ between CL patients with prospectively determined therapeutic cure or failure.

This will be achieved through the characterization of the in vivo dynamics of blood-derived monocyte, neutrophil and eosinophil transcriptome before, during and at the end of treatment in CL patients. Cell-type specific transcriptomes, composite signatures and time-response expression profiles will be contrasted among patients with therapeutic cure or failure.

4.1 Preparation

To address these, I added to the end of the sample sheet columns named ‘condition’, ‘batch’, ‘donor’, and ‘time’. These are filled in with shorthand values according to the above.

4.2 Global view

Before addressing the questions explicitly by subsetting the data, I want to get a look at the samples as they are.

hs_valid <- hs_valid %>%
  set_expt_batches(fact = "cellssource") %>%
  set_expt_conditions(fact = "typeofcells") %>%
  set_expt_samplenames(newnames = pData(hs_valid)[["samplename"]])

all_norm <- sm(normalize_expt(hs_valid, transform = "log2", norm = "quant",
                              convert = "cpm", filter = TRUE))

all_pca <- plot_pca(all_norm, plot_labels = FALSE, plot_title = "PCA - Cell type")
pp(file = glue("images/tmrc3_pca_nolabels-v{ver}.png"), image = all_pca$plot)

write.csv(all_pca$table, file = "coords/hs_donor_pca_coords.csv")
plot_corheat(all_norm, plot_title = "Heirarchical clustering - cell types")$plot

4.3 Examine samples relevant to clinical outcome

Now let us consider only the samples for which we have a clinical outcome. These fall primarily into either ‘cured’ or ‘failed’, but some people have not yet returned to the clinic after the first or second visit. These are deemed ‘lost’.

chosen_colors <- c("#D95F02", "#7570B3", "#1B9E77", "#FF0000")
names(chosen_colors) <- c("cure", "failure", "lost", "null")
newnames <- pData(hs_valid)[["samplename"]]

hs_clinical <- hs_valid %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells") %>%
  set_expt_colors(colors = chosen_colors) %>%
  set_expt_samplenames(newnames = newnames) %>%
  subset_expt(subset = "typeofcells!='PBMCs'&typeofcells!='Macrophages'")

## There were 105, now there are 87 samples.

hs_clinical_norm <- sm(normalize_expt(hs_clinical, filter = TRUE, transform = "log2",
                                      convert = "cpm", norm = "quant"))
clinical_pca <- plot_pca(hs_clinical_norm, plot_labels = FALSE, cis = NULL,
                         plot_title = "PCA - clinical samples")
pp(file = glue("images/all_clinical_nobatch_pca-v{ver}.png"), image = clinical_pca$plot,
   height = 8, width = 20)

4.3.1 Repeat without the biopsy samples

hs_clinical_nobiop <- hs_clinical %>%
  subset_expt(subset = "typeofcells!='Biopsy'")

## There were 87, now there are 55 samples.

hs_clinical_nobiop_norm <- sm(normalize_expt(hs_clinical_nobiop, filter = TRUE, transform = "log2",
                                             convert = "cpm", norm = "quant"))
clinical_nobiop_pca <- plot_pca(hs_clinical_nobiop_norm, plot_labels = FALSE, cis = NULL,
                                plot_title = "PCA - clinical samples without biopsies")
pp(file = glue("images/all_clinical_nobiop_nobatch_pca-v{ver}.png"),
   image = clinical_nobiop_pca$plot)

4.3.2 Attempt to correct for the surrogate variables

At this time we have two primary data structures of interest: hs_clinical and hs_clinical_nobiop

hs_clinical_nb <- normalize_expt(hs_clinical, filter = TRUE, batch = "svaseq",
                                 transform = "log2", convert = "cpm")

## Removing 5457 low-count genes (14484 remaining).

## batch_counts: Before batch/surrogate estimation, 72813 entries are x==0: 6%.

## batch_counts: Before batch/surrogate estimation, 220887 entries are 0<x<1: 18%.

## Setting 14890 low elements to zero.

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

clinical_batch_pca <- plot_pca(hs_clinical_nb, plot_labels = FALSE, cis = NULL,
                               plot_title = "PCA - clinical samples")
clinical_batch_pca$plot

hs_clinical_nobiop_nb <- sm(normalize_expt(hs_clinical_nobiop, filter = TRUE, batch = "svaseq",
                                           transform = "log2", convert = "cpm"))
clinical_nobiop_batch_pca <- plot_pca(hs_clinical_nobiop_nb,
                                      plot_title = "PCA - clinical samples without biopsies",
                                      plot_labels = FALSE)
pp(file = "images/clinical_batch.png", image = clinical_nobiop_batch_pca$plot)

clinical_nobiop_batch_tsne <- plot_tsne(hs_clinical_nobiop_nb,
                                        plot_title = "tSNE - clinical samples without biopsies",
                                        plot_labels = FALSE)
clinical_nobiop_batch_tsne$plot

test <- simple_varpart(hs_clinical_nobiop)

## Loading required package: Matrix

## 
## Total:109 s

test$partition_plot

4.4 Perform DE of the clinical samples cure vs. fail

individual_celltypes <- subset_expt(hs_clinical_nobiop, subset="condition!='lost'")

## There were 55, now there are 40 samples.

hs_clinic_de <- sm(all_pairwise(individual_celltypes, model_batch = "svaseq", filter = TRUE))
hs_clinic_table <- sm(combine_de_tables(
    hs_clinic_de,
    excel = glue::glue("excel/individual_celltypes_table-v{ver}.xlsx")))

hs_clinic_sig <- sm(extract_significant_genes(
    hs_clinic_table,
    excel = glue::glue("excel/individual_celltypes_sig-v{ver}.xlsx")))

hs_clinic_sig[["summary_df"]]

##                 limma_change_counts_up limma_change_counts_down
## failure_vs_cure                    230                      369
##                 edger_change_counts_up edger_change_counts_down
## failure_vs_cure                    347                      322
##                 deseq_change_counts_up deseq_change_counts_down
## failure_vs_cure                    290                      371
##                 ebseq_change_counts_up ebseq_change_counts_down
## failure_vs_cure                     84                      268
##                 basic_change_counts_up basic_change_counts_down
## failure_vs_cure                     55                       69

4.4.1 Look at only the differential genes

A good suggestion from Theresa was to examine only the most variant genes from failure vs. cure and see how they change the clustering/etc results. This is my attempt to address this query.

hs_clinic_topn <- sm(extract_significant_genes(hs_clinic_table, n = 100))
table <- "failure_vs_cure"
wanted <- rbind(hs_clinic_topn[["deseq"]][["ups"]][[table]],
                hs_clinic_topn[["deseq"]][["downs"]][[table]])

small_expt <- exclude_genes_expt(hs_clinical_nobiop, ids = rownames(wanted), method = "keep")

## Before removal, there were 19941 genes, now there are 200.

## There are 55 samples which kept less than 90 percent counts.

##  1017n1  1017m1  1034n1  1034n2  1034m2 1034m2-  2052e1  2052m2  2052n2  2052m3 
##  0.1412  0.2526  1.1251  1.3798  0.8174  0.8127  0.4081  0.4147  0.3047  0.5205 
##  2052n3  2065m1  2065n1  2066m1  2066n1  2065m2  2065n2  2066m2  2068m1  2068n1 
##  0.3258  0.8297  1.1489  0.3545  0.2052  0.6134  0.5544  0.5485  0.5715  0.5106 
##  2068e1  2072m1  2072n1  2072e1  2073m1  2073n1  2073e1  2068m2  2068n2  2068e2 
##  1.0353  0.5206  0.3950  0.6709  0.5844  0.7572  0.5837  0.2638  0.2737  0.3495 
##  2072m2  2072n2  2072e2  2073m2  2073n2  2073e2  2068m3  2068n3  2068e3  2072m3 
##  0.2881  0.2474  0.3760  0.7659  1.3507  0.4852  0.3715  0.3559  0.4352  0.5425 
##  2072n3  2072e3  2073m3  2073n3  2073e3  2162m1  2162n1  2162e1  2168e1  2168m2 
##  0.8864  0.5136  0.4433  0.3778  0.3168  0.4533  0.3589  0.6662  0.7565  0.7862 
##  2168n2  2168e2  2168m3  2168n3  2168e3 
##  1.2432  0.6534  0.8916  1.4272  0.7247

small_norm <- sm(normalize_expt(small_expt, transform = "log2", convert = "cpm",
                                norm = "quant", filter = TRUE))
plot_pca(small_norm)$plot

small_nb <- normalize_expt(small_expt, transform = "log2", convert = "cpm",
                           batch = "svaseq", norm = "quant", filter = TRUE)

## Warning in normalize_expt(small_expt, transform = "log2", convert = "cpm", :
## Quantile normalization and sva do not always play well together.

## Removing 0 low-count genes (200 remaining).

## batch_counts: Before batch/surrogate estimation, 0 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 131 entries are 0<x<1: 1%.

## Setting 8 low elements to zero.

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

plot_pca(small_nb)$plot

## DESeq2 MA plot of failure / cure
hs_clinic_table[["plots"]][["failure_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
hs_clinic_table[["plots"]][["failure_vs_cure"]][["deseq_vol_plots"]]$plot

4.4.2 g:Profiler results using the significant up and down genes

ups <- hs_clinic_sig[["deseq"]][["ups"]][[1]]
downs <- hs_clinic_sig[["deseq"]][["downs"]][[1]]

hs_clinic_gprofiler_ups <- simple_gprofiler(ups)

## Performing gProfiler GO search of 290 genes against hsapiens.

## GO search found 22 hits.

## Performing gProfiler KEGG search of 290 genes against hsapiens.

## KEGG search found 2 hits.

## Performing gProfiler REAC search of 290 genes against hsapiens.

## REAC search found 6 hits.

## Performing gProfiler MI search of 290 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 290 genes against hsapiens.

## TF search found 28 hits.

## Performing gProfiler CORUM search of 290 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 290 genes against hsapiens.

## HP search found 0 hits.

hs_clinic_gprofiler_ups[["pvalue_plots"]][["bpp_plot_over"]]

hs_clinic_gprofiler_ups[["pvalue_plots"]][["mfp_plot_over"]]

hs_clinic_gprofiler_ups[["pvalue_plots"]][["reactome_plot_over"]]

##hs_try2 <- simple_gprofiler2(ups)

hs_clinic_gprofiler_downs <- simple_gprofiler(downs)

## Performing gProfiler GO search of 371 genes against hsapiens.

## GO search found 62 hits.

## Performing gProfiler KEGG search of 371 genes against hsapiens.

## KEGG search found 4 hits.

## Performing gProfiler REAC search of 371 genes against hsapiens.

## REAC search found 4 hits.

## Performing gProfiler MI search of 371 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 371 genes against hsapiens.

## TF search found 21 hits.

## Performing gProfiler CORUM search of 371 genes against hsapiens.

## CORUM search found 3 hits.

## Performing gProfiler HP search of 371 genes against hsapiens.

## HP search found 6 hits.

hs_clinic_gprofiler_downs[["pvalue_plots"]][["bpp_plot_over"]]

hs_clinic_gprofiler_downs[["pvalue_plots"]][["mfp_plot_over"]]

hs_clinic_gprofiler_downs[["pvalue_plots"]][["reactome_plot_over"]]

4.5 Perform GSVA on the clinical samples

hs_celltype_gsva_c2 <- sm(simple_gsva(individual_celltypes))
hs_celltype_gsva_c2_sig <- sm(get_sig_gsva_categories(
    hs_celltype_gsva_c2,
    excel = "excel/individual_celltypes_gsva_c2.xlsx"))

broad_c7 <- GSEABase::getGmt("reference/msigdb/c7.all.v7.2.entrez.gmt",
                             collectionType = GSEABase::BroadCollection(category = "c7"),
                             geneIdType = GSEABase::EntrezIdentifier())
hs_celltype_gsva_c7 <- sm(simple_gsva(individual_celltypes, signatures = broad_c7,
                                      msig_xml = "reference/msigdb_v7.2.xml", cores = 10))
hs_celltype_gsva_c7_sig <- sm(get_sig_gsva_categories(
    hs_celltype_gsva_c7,
    excel = "excel/individual_celltypes_gsva_c7.xlsx"))

4.5.1 Print some plots of the GSVA outputs

## The raw heatmap of the C2 values
hs_celltype_gsva_c2_sig[["raw_plot"]]

## The 'significance' scores of the C2 values
hs_celltype_gsva_c2_sig[["score_plot"]]

## The subset of scores for categories deemed significantly different.
hs_celltype_gsva_c2_sig[["subset_plot"]]

## The raw heatmap of the C7 values
hs_celltype_gsva_c7_sig[["raw_plot"]]

## The 'significance' scores of the C7 values
hs_celltype_gsva_c7_sig[["score_plot"]]

## The subset of scores for categories deemed significantly different.
hs_celltype_gsva_c7_sig[["subset_plot"]]

5 Individual Cell types

The following blocks split the samples into a few groups by sample type and look at the distributions between them.

5.1 Implementation details

Get top/bottom n genes for each cell type, using clinical outcome as the factor of interest. For the moment, use sva for the DE analysis. Provide cpms for the top/bottom n genes.

Start with top/bottom 200. Perform default logFC and p-value as well.

5.1.1 Shared contrasts

Here is the contrast we will use throughput, I am leaving open the option to add more.

keepers <- list(
  "fail_vs_cure" = c("failure", "cure"))

5.2 Monocytes

5.2.1 Evaluate Monocyte samples

mono <- subset_expt(hs_valid, subset = "typeofcells=='Monocytes'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "donor") %>%
  set_expt_colors(colors = chosen_colors)

## There were 105, now there are 21 samples.

## FIXME set_expt_colors should speak up if there are mismatches here!!!

save_result <- save(mono, file = "rda/monocyte_expt.rda")
mono_norm <- normalize_expt(mono, convert = "cpm", filter = TRUE,
                            transform = "log2", norm = "quant")

## Removing 8968 low-count genes (10973 remaining).

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

plt <- plot_pca(mono_norm, plot_labels = FALSE)$plot
pp(file = glue("images/mono_pca_normalized-v{ver}.pdf"), image = plt)

mono_nb <- normalize_expt(mono, convert = "cpm", filter = TRUE,
                          transform = "log2", batch = "svaseq")

## Removing 8968 low-count genes (10973 remaining).

## batch_counts: Before batch/surrogate estimation, 1144 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 12729 entries are 0<x<1: 6%.

## Setting 368 low elements to zero.

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

plt <- plot_pca(mono_nb, plot_labels = FALSE)$plot
pp(file = glue("images/mono_pca_normalized_batch-v{ver}.pdf"), image = plt)

5.2.2 DE of Monocyte samples

5.2.2.1 Without sva

mono_de <- sm(all_pairwise(mono, model_batch = FALSE, filter = TRUE))

mono_tables <- sm(combine_de_tables(
    mono_de, keepers = keepers,
    excel = glue::glue("excel/monocyte_clinical_all_tables-v{ver}.xlsx")))

written <- write_xlsx(data = mono_tables[["data"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_table-v{ver}.xlsx"))
mono_sig <- sm(extract_significant_genes(mono_tables, according_to = "deseq"))
written <- write_xlsx(data = mono_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_sigup-v{ver}.xlsx"))
written <- write_xlsx(data = mono_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_sigdown-v{ver}.xlsx"))

mono_pct_sig <- sm(extract_significant_genes(mono_tables, n = 200,
                                             lfc = NULL, p = NULL, according_to = "deseq"))
written <- write_xlsx(data = mono_pct_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_sigup_pct-v{ver}.xlsx"))
written <- write_xlsx(data = mono_pct_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/monocyte_clinical_sigdown_pct-v{ver}.xlsx"))
mono_sig$summary_df

## data frame with 0 columns and 1 row

## Print out a table of the cpm values for other explorations.
mono_cpm <- sm(normalize_expt(mono, convert = "cpm"))
written <- write_xlsx(data = exprs(mono_cpm),
                      excel = glue::glue("excel/monocyte_cpm_before_batch-v{ver}.xlsx"))
mono_bcpm <- sm(normalize_expt(mono, filter = TRUE, convert = "cpm", batch = "svaseq"))
written <- write_xlsx(data = exprs(mono_bcpm),
                      excel = glue::glue("excel/monocyte_cpm_after_batch-v{ver}.xlsx"))

5.2.2.2 With sva

mono_de_sva <- sm(all_pairwise(mono, model_batch = "svaseq", filter = TRUE))

mono_tables_sva <- sm(combine_de_tables(
    mono_de_sva, keepers = keepers,
    excel = glue::glue("excel/monocyte_clinical_all_tables_sva-v{ver}.xlsx")))

mono_sig_sva <- sm(extract_significant_genes(
    mono_tables_sva,
    excel = glue::glue("excel/monocyte_clinical_sig_tables_sva-v{ver}.xlsx"),
    according_to = "deseq"))

5.2.2.3 Monocyte DE plots

First print out the DE plots without and then with sva estimates.

## DESeq2 MA plot of failure / cure
mono_tables[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
mono_tables[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## DESeq2 MA plot of failure / cure with svaseq
mono_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure with svaseq
mono_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

5.2.2.4 Monocyte ontology search

ups <- mono_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
downs <- mono_sig[["deseq"]][["downs"]][["fail_vs_cure"]]

mono_up_gp <- simple_gprofiler(sig_genes = ups)

## Performing gProfiler GO search of 149 genes against hsapiens.

## GO search found 61 hits.

## Performing gProfiler KEGG search of 149 genes against hsapiens.

## KEGG search found 3 hits.

## Performing gProfiler REAC search of 149 genes against hsapiens.

## REAC search found 1 hits.

## Performing gProfiler MI search of 149 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 149 genes against hsapiens.

## TF search found 5 hits.

## Performing gProfiler CORUM search of 149 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 149 genes against hsapiens.

## HP search found 0 hits.

mono_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

mono_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

mono_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

mono_down_gp <- simple_gprofiler(sig_genes = downs)

## Performing gProfiler GO search of 262 genes against hsapiens.

## GO search found 53 hits.

## Performing gProfiler KEGG search of 262 genes against hsapiens.

## KEGG search found 0 hits.

## Performing gProfiler REAC search of 262 genes against hsapiens.

## REAC search found 1 hits.

## Performing gProfiler MI search of 262 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 262 genes against hsapiens.

## TF search found 32 hits.

## Performing gProfiler CORUM search of 262 genes against hsapiens.

## CORUM search found 2 hits.

## Performing gProfiler HP search of 262 genes against hsapiens.

## HP search found 0 hits.

mono_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

mono_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

## NULL

mono_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

5.2.2.5 Monocyte MSigDB query

broad_c7 <- GSEABase::getGmt("reference/msigdb/c7.all.v7.2.entrez.gmt",
                             collectionType = GSEABase::BroadCollection(category = "c7"),
                             geneIdType = GSEABase::EntrezIdentifier())

mono_up_goseq_msig <- goseq_msigdb(sig_genes = ups, signatures = broad_c7,
                                   signature_category = "c7", length_db = hs_length)

## Before conversion: 149, after conversion: 147.

## Before conversion: 227921, after conversion: 35341.

## Found 132 go_db genes and 147 length_db genes out of 147.

mono_down_goseq_msig <- goseq_msigdb(sig_genes = downs, signatures = broad_c7,
                                     signature_category = "c7", length_db = hs_length)

## Before conversion: 262, after conversion: 258.

## Before conversion: 227921, after conversion: 35341.

## Found 249 go_db genes and 258 length_db genes out of 258.

5.2.2.6 Plot of similar experiments

## Monocyte genes with increased expression in the failed samples
## share genes with the following experiments
mono_up_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

## Monocyte genes with increased expression in the cured samples
## share genes with the following experiments
mono_down_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

5.2.3 Evaluate Neutrophil samples

neut <- subset_expt(hs_valid, subset = "typeofcells=='Neutrophils'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "donor") %>%
  set_expt_colors(colors = chosen_colors)

## There were 105, now there are 20 samples.

save_result <- save(neut, file = "rda/neutrophil_expt.rda")
neut_norm <- sm(normalize_expt(neut, convert = "cpm", filter = TRUE, transform = "log2"))
plt <- plot_pca(neut_norm, plot_labels = FALSE)$plot
pp(file = glue("images/neut_pca_normalized-v{ver}.pdf"), image = plt)

neut_nb <- sm(normalize_expt(neut, convert = "cpm", filter = TRUE,
                             transform = "log2", batch = "svaseq"))
plt <- plot_pca(neut_nb, plot_labels = FALSE)$plot
pp(file = glue("images/neut_pca_normalized_svaseq-v{ver}.pdf"), image = plt)

5.2.4 DE of Netrophil samples

5.2.4.1 Without sva

neut_de <- sm(all_pairwise(neut, model_batch = FALSE, filter = TRUE))

neut_tables <- sm(combine_de_tables(
    neut_de, keepers = keepers,
    excel = glue::glue("excel/neutrophil_clinical_all_tables-v{ver}.xlsx")))

written <- write_xlsx(data = neut_tables[["data"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_table-v{ver}.xlsx"))
neut_sig <- sm(extract_significant_genes(neut_tables, according_to = "deseq"))
written <- write_xlsx(data = neut_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_sigup-v{ver}.xlsx"))
written <- write_xlsx(data = neut_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_sigdown-v{ver}.xlsx"))

neut_pct_sig <- sm(extract_significant_genes(neut_tables, n = 200, lfc = NULL,
                                             p = NULL, according_to = "deseq"))
written <- write_xlsx(data = neut_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_sigup_pct-v{ver}.xlsx"))
written <- write_xlsx(data = neut_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/neutrophil_clinical_sigdown_pct-v{ver}.xlsx"))
neut_cpm <- sm(normalize_expt(neut, convert = "cpm"))
written <- write_xlsx(data = exprs(neut_cpm),
                      excel = glue::glue("excel/neutrophil_cpm_before_batch-v{ver}.xlsx"))
neut_bcpm <- sm(normalize_expt(neut, filter = TRUE, batch = "svaseq", convert = "cpm"))
written <- write_xlsx(data = exprs(neut_bcpm),
                      excel = glue::glue("excel/neutrophil_cpm_after_batch-v{ver}.xlsx"))

5.2.4.2 With sva

neut_de_sva <- sm(all_pairwise(neut, model_batch = "svaseq", filter = TRUE))

neut_tables_sva <- sm(combine_de_tables(
    neut_de_sva, keepers = keepers,
    excel = glue::glue("excel/neutrophil_clinical_all_tables_sva-v{ver}.xlsx")))

neut_sig_sva <- sm(extract_significant_genes(
    neut_tables_sva,
    excel = glue::glue("excel/neutrophil_clinical_sig_tables_sva-v{ver}.xlsx"),
    according_to = "deseq"))

5.2.4.3 Neutrophil DE plots

## DESeq2 MA plot of failure / cure
neut_tables[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
neut_tables[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## DESeq2 MA plot of failure / cure with sva
neut_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure with sva
neut_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

5.2.4.4 Neutrophil ontology search

ups <- neut_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
downs <- neut_sig[["deseq"]][["downs"]][["fail_vs_cure"]]

neut_up_gp <- simple_gprofiler(sig_genes = ups)

## Performing gProfiler GO search of 271 genes against hsapiens.

## GO search found 5 hits.

## Performing gProfiler KEGG search of 271 genes against hsapiens.

## KEGG search found 3 hits.

## Performing gProfiler REAC search of 271 genes against hsapiens.

## REAC search found 36 hits.

## Performing gProfiler MI search of 271 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 271 genes against hsapiens.

## TF search found 41 hits.

## Performing gProfiler CORUM search of 271 genes against hsapiens.

## CORUM search found 6 hits.

## Performing gProfiler HP search of 271 genes against hsapiens.

## HP search found 0 hits.

neut_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

neut_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

neut_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

neut_down_gp <- simple_gprofiler(downs)

## Performing gProfiler GO search of 260 genes against hsapiens.

## GO search found 3 hits.

## Performing gProfiler KEGG search of 260 genes against hsapiens.

## KEGG search found 0 hits.

## Performing gProfiler REAC search of 260 genes against hsapiens.

## REAC search found 1 hits.

## Performing gProfiler MI search of 260 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 260 genes against hsapiens.

## TF search found 40 hits.

## Performing gProfiler CORUM search of 260 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 260 genes against hsapiens.

## HP search found 2 hits.

neut_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

neut_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

## NULL

neut_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

5.2.4.5 Neutrophil GSVA query

neut_up_goseq_msig <- goseq_msigdb(sig_genes = ups, signatures = broad_c7,
                                   signature_category = "c7", length_db = hs_length)

## Before conversion: 271, after conversion: 267.

## Before conversion: 227921, after conversion: 35341.

## Found 251 go_db genes and 267 length_db genes out of 267.

neut_down_goseq_msig <- goseq_msigdb(sig_genes = downs, signatures = broad_c7,
                                     signature_category = "c7", length_db = hs_length)

## Before conversion: 260, after conversion: 255.

## Before conversion: 227921, after conversion: 35341.

## Found 247 go_db genes and 255 length_db genes out of 255.

5.2.4.6 Plot of similar experiments

## Neutrophil genes with increased expression in the failed samples
## share genes with the following experiments
neut_up_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

## Neutrophil genes with increased expression in the cured samples
## share genes with the following experiments
neut_down_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

5.3 Eosinophils

5.3.1 Evaluate Eosinophil samples

eo <- subset_expt(hs_valid, subset = "typeofcells=='Eosinophils'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "donor") %>%
  set_expt_colors(colors = chosen_colors)

## There were 105, now there are 14 samples.

save_result <- save(eo, file = "rda/eosinophil_expt.rda")
eo_norm <- sm(normalize_expt(eo, convert = "cpm", transform = "log2",
                             norm = "quant", filter = TRUE))
plt <- plot_pca(eo_norm, plot_labels = FALSE)$plot
pp(file = glue("images/eo_pca_normalized-v{ver}.pdf"), image = plt)

eo_nb <- sm(normalize_expt(eo, convert = "cpm", transform = "log2",
                           filter = TRUE, batch = "svaseq"))
plot_pca(eo_nb)$plot

5.3.2 DE of Eosinophil samples

5.3.2.1 Withouth sva

eo_de <- sm(all_pairwise(eo, model_batch = FALSE, filter = TRUE))

eo_tables <- sm(combine_de_tables(
    eo_de, keepers = keepers,
    excel = glue::glue("excel/eosinophil_clinical_all_tables-v{ver}.xlsx")))

written <- write_xlsx(data = eo_tables[["data"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_table-v{ver}.xlsx"))
eo_sig <- sm(extract_significant_genes(eo_tables, according_to = "deseq"))
written <- write_xlsx(data = eo_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_sigup-v{ver}.xlsx"))
written <- write_xlsx(data = eo_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_sigdown-v{ver}.xlsx"))

eo_pct_sig <- sm(extract_significant_genes(eo_tables, n = 200,
                                           lfc = NULL, p = NULL, according_to = "deseq"))
written <- write_xlsx(data = eo_pct_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_sigup_pct-v{ver}.xlsx"))
written <- write_xlsx(data = eo_pct_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/eosinophil_clinical_sigdown_pct-v{ver}.xlsx"))

eo_cpm <- sm(normalize_expt(eo, convert = "cpm"))
written <- write_xlsx(data = exprs(eo_cpm),
                      excel = glue::glue("excel/eosinophil_cpm_before_batch-v{ver}.xlsx"))
eo_bcpm <- sm(normalize_expt(eo, filter = TRUE, batch = "svaseq", convert = "cpm"))
written <- write_xlsx(data = exprs(eo_bcpm),
                      excel = glue::glue("excel/eosinophil_cpm_after_batch-v{ver}.xlsx"))

5.3.2.2 With sva

eo_de_sva <- sm(all_pairwise(eo, model_batch = "svaseq", filter = TRUE))

eo_tables_sva <- sm(combine_de_tables(
    eo_de_sva, keepers = keepers,
    excel = glue::glue("excel/eosinophil_clinical_all_tables_sva-v{ver}.xlsx")))

eo_sig_sva <- sm(extract_significant_genes(
    eo_tables_sva,
    excel = glue::glue("excel/eosinophil_clinical_sig_tables_sva-v{ver}.xlsx"),
    according_to = "deseq"))

5.3.2.3 Eosinophil DE plots

## DESeq2 MA plot of failure / cure
eo_tables[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
eo_tables[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## DESeq2 MA plot of failure / cure with sva
eo_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure with sva
eo_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

5.3.2.4 Eosinophil ontology search

ups <- eo_sig[["deseq"]][["ups"]][["fail_vs_cure"]]
downs <- eo_sig[["deseq"]][["downs"]][["fail_vs_cure"]]

eo_up_gp <- simple_gprofiler(sig_genes = ups)

## Performing gProfiler GO search of 212 genes against hsapiens.

## GO search found 53 hits.

## Performing gProfiler KEGG search of 212 genes against hsapiens.

## KEGG search found 8 hits.

## Performing gProfiler REAC search of 212 genes against hsapiens.

## REAC search found 5 hits.

## Performing gProfiler MI search of 212 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 212 genes against hsapiens.

## TF search found 44 hits.

## Performing gProfiler CORUM search of 212 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 212 genes against hsapiens.

## HP search found 1 hits.

eo_up_gp[["pvalue_plots"]][["bpp_plot_over"]]

eo_up_gp[["pvalue_plots"]][["mfp_plot_over"]]

eo_up_gp[["pvalue_plots"]][["reactome_plot_over"]]

eo_down_gp <- simple_gprofiler(downs)

## Performing gProfiler GO search of 154 genes against hsapiens.

## GO search found 1 hits.

## Performing gProfiler KEGG search of 154 genes against hsapiens.

## KEGG search found 1 hits.

## Performing gProfiler REAC search of 154 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 154 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 154 genes against hsapiens.

## TF search found 0 hits.

## Performing gProfiler CORUM search of 154 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 154 genes against hsapiens.

## HP search found 0 hits.

eo_down_gp[["pvalue_plots"]][["bpp_plot_over"]]

eo_down_gp[["pvalue_plots"]][["mfp_plot_over"]]

## NULL

eo_down_gp[["pvalue_plots"]][["reactome_plot_over"]]

## NULL

5.3.2.5 Eosinophil MSigDB query

eo_up_goseq_msig <- goseq_msigdb(sig_genes = ups, signatures = broad_c7,
                                 signature_category = "c7", length_db = hs_length)

## Before conversion: 212, after conversion: 211.

## Before conversion: 227921, after conversion: 35341.

## Found 183 go_db genes and 211 length_db genes out of 211.

eo_down_goseq_msig <- goseq_msigdb(sig_genes = downs, signatures = broad_c7,
                                   signature_category = "c7", length_db = hs_length)

## Before conversion: 154, after conversion: 152.

## Before conversion: 227921, after conversion: 35341.

## Found 139 go_db genes and 152 length_db genes out of 152.

5.3.2.6 Plot of similar experiments

## Eosinophil genes with increased expression in the failed samples
## share genes with the following experiments
eo_up_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

## Eosinophil genes with increased expression in the cured samples
## share genes with the following experiments
eo_down_goseq_msig[["pvalue_plots"]][["mfp_plot_over"]]

5.4 Biopsies

5.4.1 Evaluate Biopsy samples

biop <- subset_expt(hs_valid, subset = "typeofcells=='Biopsy'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "donor") %>%
  set_expt_colors(colors = chosen_colors)

## There were 105, now there are 32 samples.

save_result <- save(biop, file = "rda/biopsy_expt.rda")
biop_norm <- normalize_expt(biop, filter = TRUE, convert = "cpm",
                            transform = "log2", norm = "quant")

## Removing 5950 low-count genes (13991 remaining).

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

plt <- plot_pca(biop_norm, plot_labels = FALSE)$plot
pp(file = glue("images/biop_pca_normalized-v{ver}.pdf"), image = plt)

biop_nb <- sm(normalize_expt(biop, convert = "cpm", filter = TRUE,
                             transform = "log2", batch = "svaseq"))
plt <- plot_pca(biop_nb, plot_labels = FALSE)$plot
pp(file = glue("images/biop_pca_normalized_svaseq-v{ver}.pdf"), image = plt)

5.4.2 DE of Biopsy samples

5.4.2.1 Without sva

biop_de <- sm(all_pairwise(biop, model_batch = FALSE, filter = TRUE))

biop_tables <- combine_de_tables(biop_de, keepers = keepers,
                                 excel = glue::glue("excel/biopsy_clinical_all_tables-v{ver}.xlsx"))

## Deleting the file excel/biopsy_clinical_all_tables-v202104.xlsx before writing the tables.

written <- write_xlsx(data = biop_tables[["data"]][[1]],
                      excel = glue::glue("excel/biopsy_clinical_table-v{ver}.xlsx"))
biop_sig <- extract_significant_genes(biop_tables, according_to = "deseq")
##written <- write_xlsx(data = biop_sig[["deseq"]][["ups"]][[1]],
##                      excel = glue::glue("excel/biopsy_clinical_sigup-v{ver}.xlsx"))
written <- write_xlsx(data = biop_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/biopsy_clinical_sigdown-v{ver}.xlsx"))
biop_pct_sig <- extract_significant_genes(biop_tables, n = 200, lfc = NULL, p = NULL, according_to = "deseq")

## Getting the top and bottom 200 genes.

written <- write_xlsx(data = biop_pct_sig[["deseq"]][["ups"]][[1]],
                      excel = glue::glue("excel/biopsy_clinical_sigup_pct-v{ver}.xlsx"))
written <- write_xlsx(data = biop_pct_sig[["deseq"]][["downs"]][[1]],
                      excel = glue::glue("excel/biopsy_clinical_sigdown_pct-v{ver}.xlsx"))

biop_cpm <- sm(normalize_expt(biop, convert = "cpm"))
written <- write_xlsx(data = exprs(biop_cpm),
                      excel = glue::glue("excel/biopsy_cpm_before_batch-v{ver}.xlsx"))
biop_bcpm <- sm(normalize_expt(biop, filter = TRUE, batch = "svaseq", convert = "cpm"))
written <- write_xlsx(data = exprs(biop_bcpm),
                      excel = glue::glue("excel/biopsy_cpm_after_batch-v{ver}.xlsx"))

5.4.2.2 with sva

biop_de_sva <- sm(all_pairwise(biop, model_batch = "svaseq", filter = TRUE))

biop_tables_sva <- sm(combine_de_tables(
    biop_de_sva, keepers = keepers,
    excel = glue::glue("excel/biopsy_clinical_all_tables_sva-v{ver}.xlsx")))

biop_sig_sva <- sm(extract_significant_genes(
    biop_tables_sva,
    excel = glue::glue("excel/biopsy_clinical_sig_tables_sva-v{ver}.xlsx"),
    according_to = "deseq"))

5.4.2.3 Biopsy DE plots

## DESeq2 MA plot of failure / cure
biop_tables[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
biop_tables[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

## DESeq2 MA plot of failure / cure
biop_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]]$plot

## DESeq2 Volcano plot of failure / cure
biop_tables_sva[["plots"]][["fail_vs_cure"]][["deseq_vol_plots"]]$plot

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  message(paste0("Saving to ", savefile))
  tmp <- sm(saveme(filename = savefile))
}

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 75c439759cfc03b66dda979814171e43e69506ca

## This is hpgltools commit: Thu Apr 22 12:37:06 2021 -0400: 75c439759cfc03b66dda979814171e43e69506ca

## Saving to tmrc3_02sample_estimation_v202104.rda.xz

tmp <- loadme(filename = savefile)