TMRC3 Comprehensive Data Analysis: 202106

samplesheet <- "sample_sheets/tmrc3_samples_20210620.xlsx"

1 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")

2 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.

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.

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.
sanitize_columns <- c("visitnumber", "clinicaloutcome", "donor",
                      "typeofcells", "clinicalpresentation",
                      "condition", "batch")
hs_expt <- 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", meta_column = "ncrna_lost") %>%
  sanitize_expt_metadata(columns = sanitize_columns) %>%
  set_expt_factors(columns = sanitize_columns, class = "factor")

## Reading the sample metadata.

## Dropped 71 rows from the sample metadata because they were blank.

## The sample definitions comprises: 173 rows(samples) and 76 columns(metadata fields).

## Warning in create_expt(samplesheet, file_column = "hg38100hisatfile", savefile =
## glue::glue("rda/hs_expt_all-v{ver}.rda"), : Some samples were removed when cross
## referencing the samples against the count data.

## Matched 21452 annotations and counts.

## Bringing together the count matrix and gene information.

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

## The final expressionset has 21481 rows and 171 columns.

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

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

## TMRC30015 TMRC30017 TMRC30019 TMRC30044 TMRC30045 TMRC30154 TMRC30097 TMRC30075 
##     79.24     85.72     89.75     80.34     73.33     83.20     89.90     86.97 
## TMRC30087 TMRC30101 TMRC30104 TMRC30114 TMRC30127 TMRC30120 TMRC30128 TMRC30131 
##     83.63     88.41     80.29     87.62     89.49     79.16     82.53     86.82 
## TMRC30073 
##     89.26

levels(pData(hs_expt[["expressionset"]])[["visitnumber"]]) <- list(
    '0' = "notapplicable", '1' = 1, '2' = 2, '3' = 3)

Split this data into CDS and lncRNA. Oh crap in order to do that I need to recount the data. Running now (20210518)

## lnc_expt <- create_expt(samplesheet,
##                         file_column = "hg38100lncfile",
##                         gene_info = hs_annot)

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: 143 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ncrna_lost_df <- as.data.frame(pData(hs_expt)[["ncrna_lost"]])
rownames(ncrna_lost_df) <- rownames(pData(hs_expt))
colnames(ncrna_lost_df) <- "ncrna_lost"

tmpdf <- merge(nonzero$table, ncrna_lost_df, by = "row.names")
rownames(tmpdf) <- tmpdf[["Row.names"]]
tmpdf[["Row.names"]] <- NULL

ggplot(tmpdf, aes(x=ncrna_lost, y=nonzero_genes)) +
  ggplot2::geom_point() +
  ggplot2::ggtitle("Nonzero genes with respect to percent counts 
lost when ncRNA was removed.")

Najib doesn’t want this plot, but I am using it to check new samples, so will hide it from general use.

libsize <- plot_libsize(hs_expt)
libsize$plot

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

## subset_expt(): There were 171, now there are 168 samples.

valid_write <- sm(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.

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

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", size_column = "visitnumber")
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’.

hs_clinical <- hs_valid %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells") %>%
  subset_expt(subset = "typeofcells!='pbmcs'&typeofcells!='macrophages'")

## subset_expt(): There were 168, now there are 148 samples.

chosen_colors <- c("#D95F02", "#7570B3", "#1B9E77", "#FF0000", "#FF0000")
names(chosen_colors) <- c("cure", "failure", "lost", "null", "notapplicable")
hs_clinical <- set_expt_colors(hs_clinical, colors = chosen_colors)

## Warning in set_expt_colors(hs_clinical, colors = chosen_colors): Colors for the
## following categories are not being used: nullnotapplicable.

newnames <- make.names(pData(hs_clinical)[["samplename"]], unique = TRUE)
hs_clinical <- set_expt_samplenames(hs_clinical, newnames = newnames)

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,
                         size_column = "visitnumber", 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'")

## subset_expt(): There were 148, now there are 98 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 5255 low-count genes (14686 remaining).

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

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

## Setting 22846 low elements to zero.

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

clinical_batch_pca <- plot_pca(hs_clinical_nb, plot_labels = FALSE, cis = NULL,
                               size_column = "visitnumber", 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)

test <- plot_pca(hs_clinical_nobiop_nb, size_column = "visitnumber",
                 plot_title = "PCA - clinical samples without biopsies",
                 plot_labels = FALSE)
test$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

4.3.2.1 Look at remaining variance with variancePartition

test <- simple_varpart(hs_clinical_nobiop)

## 
## Total:104 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'")

## subset_expt(): There were 98, now there are 83 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_V1 limma_V2 edger_V1 edger_V2 deseq_V1 deseq_V2 ebseq_V1 ebseq_V2
## 1      142      110      207      174      201      189      100      166
##   basic_V1 basic_V2
## 1       57       19

hs_clinic_de[["comparison"]][["heat"]]

## NULL

4.4.1 Perform LRT with the clinical samples

I am not sure if we have enough samples across the three visit to completely work as well as we would like, but there is only 1 way to find out! Now that I think about it, one thing which might be awesome is to use cell type as an interacting factor…

4.4.1.1 With biopsy samples

I figure this might be a place where the biopsy samples might prove useful.

clinical_nolost <- subset_expt(hs_clinical, subset="condition!='lost'")

## subset_expt(): There were 148, now there are 131 samples.

lrt_visit_clinical_test <- sm(deseq_lrt(clinical_nolost, transform = "vst",
                                        interactor_column = "visitnumber",
                                        interest_column = "clinicaloutcome"))

## Error in nrow(lrg_significant): object 'lrg_significant' not found

lrt_visit_clinical_test[["favorite_genes"]]

## Error in eval(expr, envir, enclos): object 'lrt_visit_clinical_test' not found

lrt_celltype_clinical_test <- sm(deseq_lrt(clinical_nolost, transform = "vst",
                                           interactor_column = "typeofcells",
                                           interest_column = "clinicaloutcome"))

## Error in nrow(lrg_significant): object 'lrg_significant' not found

summary(lrt_celltype_clinical_test[["favorite_genes"]])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'lrt_celltype_clinical_test' not found

4.4.2 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 98 samples which kept less than 90 percent counts.

##  X1017n1  X1017m1  X1034n1  X1034n2  X1034m2 X1034m2.  X2052e1  X2052m2 
##   0.7276   0.4162   2.7013   3.0884   1.2799   1.2513   0.5507   0.5749 
##  X2052n2  X2052m3  X2052n3  X2065m1  X2065n1  X2066m1  X2066n1  X2065m2 
##   1.8012   0.6464   1.5493   1.0218   2.9441   0.3853   1.0215   0.5777 
##  X2065n2  X2065e2  X2066m2  X2066n2  X2066e2  X2068m1  X2068n1  X2068e1 
##   0.7291   1.2382   0.5166   0.8269   0.9474   0.5516   0.7353   1.1808 
##  X2072m1  X2072n1  X2072e1  X2071m1  X2071n1  X2071e1  X2073m1  X2073n1 
##   0.5416   0.7864   2.1099   0.7446   1.8134   0.9283   0.6518   1.9317 
##  X2073e1  X2068m2  X2068n2  X2068e2  X2072m2  X2072n2  X2072e2  X2073m2 
##   0.8120   0.3497   0.6980   0.7661   0.4201   1.0475   0.6019   0.9693 
##  X2073n2  X2073e2  X2066m3  X2066n3  X2066e3  X2065m3  X2065n3  X2065e3 
##   2.8379   0.7924   0.4418   0.9547   0.6221   0.3589   0.8420   0.5483 
##  X2068m3  X2068n3  X2068e3  X2072m3  X2072n3  X2072e3  X2073m3  X2073n3 
##   0.4487   0.7081   0.6811   0.7139   2.2444   0.7535   0.4998   0.8912 
##  X2073e3  X2162m1  X2162n1  X2162e1  X2162m2  X2162n2  X2162e2  X2162n3 
##   0.4817   0.4687   0.6579   0.8271   0.5095   0.7904   0.8386   0.8464 
##  X2162e3  X2167m1  X2167n1  X2167e1  X2168m1  X2168n1  X2168e1  X2168m2 
##   0.6290   0.6206   0.9449   1.0535   0.9953   2.3022   1.0149   0.9023 
##  X2168n2  X2168e2  X2167m2  X2167n2  X2167e2  X2167m3  X2167n3  X2167e3 
##   2.3564   0.8392   0.3654   0.9109   0.7656   0.7490   1.3369   1.2115 
##  X2168m3  X2168n3  X2168e3  X2172m1  X2172n1  X2172e1  X2173m1  X2173n1 
##   1.2749   3.6074   0.9863   0.3770   0.6949   0.4562   0.3371   0.4632 
##  X2172m3  X2172e3  X1168m1  X1168n1  X1168m2  X1168e2  X1168m3  X1168n3 
##   0.4955   0.7815   0.7637   1.6990   0.7457   0.5287   0.7347   1.2161 
##  X1167m3  X1167n3 
##   0.5319   1.1020

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

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

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, 12 entries are x==0: 0%.

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

## Setting 184 low elements to zero.

## transform_counts: Found 184 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.3 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 201 genes against hsapiens.

## GO search found 77 hits.

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

## KEGG search found 6 hits.

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

## REAC search found 6 hits.

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

## MI search found 0 hits.

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

## TF search found 49 hits.

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

## CORUM search found 0 hits.

## Performing gProfiler HP search of 201 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 189 genes against hsapiens.

## GO search found 71 hits.

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

## KEGG search found 3 hits.

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

## REAC search found 7 hits.

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

## MI search found 0 hits.

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

## TF search found 9 hits.

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

## CORUM search found 4 hits.

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

## HP search found 0 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)

## subset_expt(): There were 168, now there are 36 samples.

## Warning in set_expt_colors(., colors = chosen_colors): Colors for the following
## categories are not being used: nullnotapplicable.

## 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 8861 low-count genes (11080 remaining).

## transform_counts: Found 7 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 8861 low-count genes (11080 remaining).

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

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

## Setting 701 low elements to zero.

## transform_counts: Found 701 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 116 genes against hsapiens.

## GO search found 75 hits.

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

## KEGG search found 9 hits.

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

## REAC search found 9 hits.

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

## MI search found 0 hits.

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

## TF search found 29 hits.

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

## CORUM search found 1 hits.

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

## HP search found 38 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 207 genes against hsapiens.

## GO search found 78 hits.

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

## KEGG search found 2 hits.

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

## REAC search found 1 hits.

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

## MI search found 0 hits.

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

## TF search found 46 hits.

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

## CORUM search found 2 hits.

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

## HP search found 1 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: 116, after conversion: 116.

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

## Found 113 go_db genes and 116 length_db genes out of 116.

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

## Before conversion: 207, after conversion: 206.

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

## Found 197 go_db genes and 206 length_db genes out of 206.

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)

## subset_expt(): There were 168, now there are 35 samples.

## Warning in set_expt_colors(., colors = chosen_colors): Colors for the following
## categories are not being used: nullnotapplicable.

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 149 genes against hsapiens.

## GO search found 57 hits.

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

## KEGG search found 6 hits.

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

## REAC search found 5 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 46 hits.

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

## CORUM search found 0 hits.

## Performing gProfiler HP search of 149 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 103 genes against hsapiens.

## GO search found 0 hits.

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

## KEGG search found 0 hits.

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

## REAC search found 0 hits.

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

## MI search found 0 hits.

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

## TF search found 4 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

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

## NULL

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

## NULL

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

## NULL

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: 149, after conversion: 149.

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

## Found 144 go_db genes and 149 length_db genes out of 149.

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

## Before conversion: 103, after conversion: 100.

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

## Found 95 go_db genes and 100 length_db genes out of 100.

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)

## subset_expt(): There were 168, now there are 27 samples.

## Warning in set_expt_colors(., colors = chosen_colors): Colors for the following
## categories are not being used: nullnotapplicable.

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 114 genes against hsapiens.

## GO search found 101 hits.

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

## KEGG search found 6 hits.

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

## REAC search found 8 hits.

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

## MI search found 0 hits.

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

## TF search found 52 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 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 40 genes against hsapiens.

## GO search found 1 hits.

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

## KEGG search found 1 hits.

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

## REAC search found 0 hits.

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

## MI search found 0 hits.

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

## TF search found 1 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

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

## NULL

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: 114, after conversion: 114.

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

## Found 111 go_db genes and 114 length_db genes out of 114.

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

## Before conversion: 40, after conversion: 39.

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

## Found 35 go_db genes and 39 length_db genes out of 39.

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)

## subset_expt(): There were 168, now there are 50 samples.

## Warning in set_expt_colors(., colors = chosen_colors): Colors for the following
## categories are not being used: nullnotapplicable.

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

## Removing 5752 low-count genes (14189 remaining).

## transform_counts: Found 16 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-v202106.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

6 Look for shared genes among Monocytes/Neutrophils/Eosinophils

We have three variables containing the ‘significant’ DE genes for the three cell types. For this I am choosing (for the moment) to use the sva data.

## mono_sig_sva, neut_sig_sva, eo_sig_sva
sig_vectors <- list(
    "monocytes" = c(rownames(mono_sig_sva[["deseq"]][["ups"]][["fail_vs_cure"]]),
                    rownames(mono_sig_sva[["deseq"]][["downs"]][["fail_vs_cure"]])),
    "neutrophils" = c(rownames(neut_sig_sva[["deseq"]][["ups"]][["fail_vs_cure"]]),
                      rownames(neut_sig_sva[["deseq"]][["downs"]][["fail_vs_cure"]])),
    "eosinophils" =  c(rownames(eo_sig_sva[["deseq"]][["ups"]][["fail_vs_cure"]]),
                       rownames(eo_sig_sva[["deseq"]][["downs"]][["fail_vs_cure"]])))

shared_vector <- Vennerable::Venn(Sets = sig_vectors)
Vennerable::plot(shared_vector, doWeights = FALSE)

shared_ids <- shared_vector@IntersectionSets[["111"]]
shared_expt <- exclude_genes_expt(hs_clinical, ids = shared_ids, method = "keep")

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

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

##  X1017n1  X1017m1  X1034n1 X1034bp1  X1034n2  X1034m2 X1034m2. X2050bp1 
##  0.09618  0.09574  1.36052  0.45728  1.66473  0.57028  0.54732  0.19514 
## X2052bp1  X2052e1  X2052m2  X2052n2  X2052m3  X2052n3  X2065m1  X2065n1 
##  0.46300  0.04695  0.15612  0.32061  0.20004  0.28746  0.24780  0.99756 
## X2065bp1  X2066m1  X2066n1 X2066bp1  X2065m2  X2065n2  X2065e2  X2066m2 
##  0.50189  0.09884  0.18213  0.37438  0.09035  0.06769  0.04269  0.10083 
##  X2066n2  X2066e2  X2068m1  X2068n1  X2068e1 X2068bp1  X2072m1  X2072n1 
##  0.12348  0.04596  0.06585  0.08253  0.04199  0.22971  0.06952  0.08924 
##  X2072e1 X2072bp1  X2071m1  X2071n1  X2071e1 X2071bp1  X2073m1  X2073n1 
##  0.03975  0.22404  0.16805  0.64400  0.07191  0.29254  0.17698  0.73538 
##  X2073e1 X2073bp1  X2068m2  X2068n2  X2068e2  X2072m2  X2072n2  X2072e2 
##  0.10784  0.34703  0.08602  0.11732  0.05140  0.10111  0.13478  0.05578 
##  X2073m2  X2073n2  X2073e2  X2066m3  X2066n3  X2066e3  X2065m3  X2065n3 
##  0.36766  1.52085  0.24028  0.10646  0.12149  0.04346  0.09801  0.10366 
##  X2065e3  X2068m3  X2068n3  X2068e3  X2072m3  X2072n3  X2072e3 X2159bp1 
##  0.03510  0.07402  0.08672  0.04494  0.22587  0.85090  0.10602  0.28181 
##  X2073m3  X2073n3  X2073e3  X2162m1  X2162n1  X2162e1 X2162bp1  X2162m2 
##  0.16457  0.31538  0.08631  0.08197  0.09055  0.05417  0.26041  0.11265 
##  X2162n2  X2162e2  X2162n3  X2162e3  X2167m1  X2167n1  X2167e1 X2167bp1 
##  0.12397  0.07943  0.18645  0.05671  0.14654  0.07029  0.12166  0.25303 
##  X2168m1  X2168n1  X2168e1 X2168bp1  X2168m2  X2168n2  X2168e2  X2167m2 
##  0.38955  1.16962  0.24932  0.42659  0.34443  1.21392  0.21625  0.08277 
##  X2167n2  X2167e2  X2167m3  X2167n3  X2167e3  X2168m3  X2168n3  X2168e3 
##  0.08672  0.02541  0.14818  0.08081  0.08518  0.47671  1.74663  0.35040 
##  X2172m1  X2172n1  X2172e1 X2172bp1  X2173m1  X2173n1 X2173bp1  X2172m3 
##  0.08611  0.07582  0.04167  0.50228  0.08284  0.07457  0.27043  0.14335 
##  X2172e3 X1154bp1  X2008.1  X2008.2  X2008.3  X1029.1  X1029.2  X1029.3 
##  0.09091  0.28539  0.23131  0.15783  0.22867  0.22776  0.21855  0.18686 
##  X1036.1  X1036.2  X1036.3  X1037.1  X1037.2  X1037.3  X1031.1  X1031.2 
##  0.45938  0.30700  0.22476  0.36593  0.16479  0.15420  0.38822  0.21053 
##  X1031.3  X2002.1  X2002.3  X2009.1  X2009.2  X2010.1  X2010.2  X2010.3 
##  0.24431  0.33992  0.13768  0.23662  0.15738  0.29477  0.28700  0.17601 
##  X1019.1  X1019.2  X1019.3  X2004.1  X2004.2  X2004.3  X2001.1  X2001.2 
##  0.52085  0.36943  0.28509  0.25365  0.21867  0.18717  0.27900  0.18098 
##  X2001.3  X2003.1  X2003.2  X2003.3  X1168m1  X1168n1  X1168m2  X1168e2 
##  0.13415  0.53556  0.18122  0.18324  0.26484  0.52693  0.25146  0.08427 
##  X1168m3  X1168n3  X1167m3  X1167n3 
##  0.28475  0.39542  0.13990  0.13891

shared_written <- write_expt(shared_expt,
                             excel = glue::glue("excel/shared_across_celltypes-v{ver}.xlsx"))

## Deleting the file excel/shared_across_celltypes-v202106.xlsx before writing the tables.

## Writing the first sheet, containing a legend and some summary data.

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

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

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## plot labels was not set and there are more than 100 samples, disabling it.
## plot labels was not set and there are more than 100 samples, disabling it.

## 
## Total:0 s

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

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

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## 
## Total:0 s

7 Monocytes by visit

1. Can you please share with us a PCA (SVA and non-SVA) of the monocytes of the TMRC3 project, but labeling them based on the visit (V1, V2, V3)?
2. Can you please share DE lists of V1 vs V2, V1 vs V3, V1 vs. V2+V3 and V2 vs V3?

visit_colors <- chosen_colors <- c("#D95F02", "#7570B3", "#1B9E77")
names(visit_colors) <- c(1, 2, 3)
mono_visit <- subset_expt(hs_valid, subset = "typeofcells=='monocytes'") %>%
  set_expt_conditions(fact = "visitnumber") %>%
  set_expt_batches(fact = "clinicaloutcome") %>%
  set_expt_colors(colors = chosen_colors)

## subset_expt(): There were 168, now there are 36 samples.

mono_visit_norm <- normalize_expt(mono_visit, filter = TRUE, norm = "quant", convert = "cpm",
                                  transform = "log2")

## Removing 8861 low-count genes (11080 remaining).

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

mono_visit_pca <- plot_pca(mono_visit_norm)
pp(file = "images/monocyte_by_visit.png", image = mono_visit_pca$plot)

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

## Removing 8861 low-count genes (11080 remaining).

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

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

## Setting 563 low elements to zero.

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

mono_visit_nb_pca <- plot_pca(mono_visit_nb)
pp(file = "images/monocyte_by_visit_nb.png", image = mono_visit_nb_pca$plot)

table(pData(mono_visit_norm)$batch)

## 
##    cure failure    lost 
##      15      15       6

keepers <- list(
    "second_vs_first" = c("c2", "c1"),
    "third_vs_second" = c("c3", "c2"),
    "third_vs_first" = c("c3", "c1"))
mono_visit_de <- all_pairwise(mono_visit, model_batch = "svaseq", filter = TRUE)

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

## Plotting a PCA before surrogate/batch inclusion.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

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

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

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

## Setting 563 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

mono_visit_tables <- combine_de_tables(
    mono_visit_de,
    keepers = keepers,
    excel = glue::glue("excel/mono_visit_tables-v{ver}.xlsx"))

## Deleting the file excel/mono_visit_tables-v202106.xlsx before writing the tables.

new_factor <- as.character(pData(mono_visit)[["visitnumber"]])
not_one_idx <- new_factor != 1
new_factor[not_one_idx] <- "not_1"
mono_one_vs <- set_expt_conditions(mono_visit, new_factor)

mono_one_vs_de <- all_pairwise(mono_one_vs, model_batch = "svaseq", filter = TRUE)

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

## Plotting a PCA before surrogate/batch inclusion.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

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

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

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

## Setting 569 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

mono_one_vs_tables <- combine_de_tables(
    mono_one_vs_de,
    excel = glue::glue("excel/mono_one_vs_tables-v{ver}.xlsx"))

## Deleting the file excel/mono_one_vs_tables-v202106.xlsx before writing the tables.

8 Test TSP

In writing the following, I quickly realized that tspair was not joking when it said it is intended for small numbers of genes. For a full expressionset of human data it is struggling. I like the idea, it may prove worth while to spend some time optimizing the package so that it is more usable.

expt <- hs_clinical_nobiop

simple_tsp <- function(expt, column = "condition") {
  facts <- levels(as.factor(pData(expt)[[column]]))
  retlist <- list()
  if (length(facts) < 2) {
    stop("This requires factors with at least 2 levels.")
  } else if (length(facts) == 2) {
    retlist <- simple_tsp_pair(expt, column = column)
  } else {
    for (first in 1:(length(facts) - 1)) {
      for (second in 2:(length(facts))) {
        if (first < second) {
          name <- glue::glue("{facts[first]}_vs_{facts[second]}")
          message("Starting ", name, ".")
          substring <- glue::glue("{column}=='{facts[first]}'|{column}=='{facts[second]}'")
          subby <- subset_expt(expt, subset=as.character(substring))
          retlist[[name]] <- simple_tsp_pair(subby, column = column)
        }
      }
    }
  }
}

simple_tsp_pair <- function(subby, column = "condition", repetitions = 50) {
  tsp_input <- subby[["expressionset"]]
  tsp_output <- tspcalc(tsp_input, column)
  tsp_scores <- tspsig(tsp_input, column, B = repetitions)
}

tsp1 <- tspcalc(tsp_input, "condition")

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 f80b9b77ce0f5b0c13e9172e5cf51a94eaadaa9e

## This is hpgltools commit: Mon Jun 28 14:14:28 2021 -0400: f80b9b77ce0f5b0c13e9172e5cf51a94eaadaa9e

## Saving to tmrc3_02sample_estimation_v202106.rda.xz

tmp <- loadme(filename = savefile)