L. panamensis 202011: Reorganizing TMRC3 Analyses

1 Introduction

The set of analyses performed in tmrc3 sample estimation has become unorganized and difficult to follow. This is intended to start at the level of broad organization before stepping into the analyses.

2 Goals

These samples are from patients who either successfully cleared a Leishmania panamensis infection following treatment, or did not. They include biopsies from each patient along with purifications for Monocytes, Neutrophils, and Eosinophils. When possible, this process was repeated over three visits; but some patients did not return for the second or third visit.

The over-arching goal is to look for attributes(most likely genes) which distinguish patients who do and do not cure the infection after treatment. If possible, these will be apparent on the first visit.

2.1 Relevant Metadata

The metadata factors (in no particular order) of the experiment are:

1. Visit. 1, 2, or 3.
2. Patient.
3. Clinical outcome: cure or fail. This includes the population of patients who did not return and are labeled ‘lost’. This is, by definition, confounded with patient.
4. Drug treatment. Most of the patients were treated with an antimonial, but a few were treated with miltefosine.
5. Cell type or biopsy.

Metadata was also collected for the patients and there are many factors which have potential to affect the outcome. These analyses will generally remain agnostic for these factors, which include (among other things):

1. sex
2. ethnicity
3. age
4. ulcer/lesion attributes (by visit)
5. amount of time spent with infection

2.2 Libraries, sequencing, mapping, quantification

The samples used in these analyses were all collected, purified, and libraries generated by the scientists/doctors at CIDEIM. The sequencing libraries were generated via the TruSeq non-stranded library kit and sequenced either at JHU or UMD; earlier samples were single-ended, but most were paired.

All samples were trimmed with trimomatic using the same set of parameters. All mapping was performed with hisat2 version 2.2.1. Quantifications were also performed with salmon version 1.2.0. The reference genome used was hg38 100 (released 202004). When samples were mapped against L.panamensis, the TriTrypDB version 36 reference was used. The set of annotations were therefore limited to ensembl’s 2020 release. The early samples were actually first mapped with hg38 91 and later redone.

2.3 Types of analyses included

This document will limit itself to a set of canonical RNAseq analyses. It must therefore create files containing the raw data to facilitate sharing the data. It will plot metrics of the data to demonstrate the sequencing quality and clustering of the samples under the various conditions examined and normalizations employed. It will perform differential expression analyses for the metadata factors of interest alongside likelihood ratio tests for factors like celltype and time. Given the sets of over/under expressed genes observed in the various DE methods, it will perform the likely gene set tests for over represented gene ontology groups, reactome, etc. Simultaneously, the raw data will be passed to gene set variance analyses to see if there are groups of genes overrepresented in other experiments.

3 Metadata collection

There are two metadata sources:

1. The online sample sheet, which I periodically update and download into the ‘sample_sheets/’ directory.
2. The crf metadata, describing the individual patients.

samplesheet <- "sample_sheets/tmrc3_samples_20211223.xlsx"
crf_metadata <- "sample_sheets/20210825_EXP_ESPECIAL_TMRC3_VERSION_2.xlsx"

4 Annotation Collection

The primary annotation sources are:

1. Ensembl’s biomart archive from 2020 for human annotations.
2. The TriTrypDB release 36 for parasite annotations.

Both provide GO data. They also provide helpful links to other data sources. For the moment, we are focusing on the human annotations.

4.1 Gene annotations

These analyses have focused on gene-level abundances/differences. Thus, when htseq-count was invoked against the hisat2-based mappings, parameters were chosen to count genes rather than transcripts. Similarly, when salmon counts were used via tximport, a mapping of genes to transcripts was used to collapse the matrix to gene-level abundances. This decision may be revisited.

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  
##                    
##                    
##                    
## 

4.2 Gene ontology data

The set of GO categories has not been limited to the 2020 data at the time of this writing.

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

5 The complete set of data

Before we do any of the following subsets/analyses of the data, we need to collect it all in one place. Let’s do that here and name it ‘hs_valid’ meaning it is the set of valid data for Homo sapiens.

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") %>%
  set_expt_conditions(fact="clinicaloutcome") %>%
  set_expt_batches(fact="visitnumber")

## Reading the sample metadata.

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

## The sample definitions comprises: 238 rows(samples) and 82 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 224 columns.

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

## There are 19 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 TMRC30126 TMRC30127 TMRC30120 TMRC30128 
##     83.63     88.41     80.29     87.62     84.52     89.49     79.16     82.53 
## TMRC30141 TMRC30131 TMRC30073 
##     89.40     86.82     89.26

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

5.1 Add the CRF patient metadata

Let us also merge in the clinician’s metadata. I worry a little that this might not be allowed for dbGap data, but if it is a problem I suspect we can just remove the bad columns from it. Also note that I rarely use the join function, but it is somewhat required here because I do not want to risk shuffling the metadata when I add the new metadata, which comes from a spreadsheet sorted by patient, not sample. In doing this I therefore just created a new column ‘join’ which contains the shared information, e.g. the patient ID from the existing metadata and the same ID from the CRF file which has been coerced into lowercase.

hs_pd <- pData(hs_expt)
start <- rownames(hs_pd)
hs_crf <- openxlsx::read.xlsx(crf_metadata)
hs_crf[["join"]] <- tolower(hs_crf[["codigo_paciente"]])
hs_pd[["join"]] <- hs_pd[["tubelabelorigin"]]
test <- plyr::join(hs_pd, hs_crf, by="join")
test[["join"]] <- NULL
rownames(test) <- rownames(hs_pd)
na_idx <- is.na(test)
test[na_idx] <- "undefined"
pData(hs_expt) <- test

5.2 Set our initial coverage goal

There exists a baseline coverage below which we do not wish to fall. One likely way to approach it heuristically is to assume we should observe some number of genes. With that in mind, I arbitrarily chose 11,000 non-zero genes as the minimum.

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 224, now there are 221 samples.

5.3 Set up initial data subsets of interest

One of the first global metrics I would like to provide is the set of library sizes. Unfortunately, we have too many samples to fit on a screen now. Therefore, I am going to do an early split of the data by cell type in the hopes that doing so will make it possible to visualize the library sizes/nonzero genes.

The initial factor for this is ‘typeofcells’.

table(pData(hs_valid)[["typeofcells"]])

## 
##      biopsy eosinophils macrophages   monocytes neutrophils       pbmcs 
##          52          35          14          58          56           6

biopsy_samples <- subset_expt(hs_valid, subset="typeofcells=='biopsy'")

## subset_expt(): There were 221, now there are 52 samples.

eosinophil_samples <- subset_expt(hs_valid, subset="typeofcells=='eosinophils'")

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

macrophage_samples <- subset_expt(hs_valid, subset="typeofcells=='macrophages'")

## subset_expt(): There were 221, now there are 14 samples.

monocyte_samples <- subset_expt(hs_valid, subset="typeofcells=='monocytes'")

## subset_expt(): There were 221, now there are 58 samples.

neutrophil_samples <- subset_expt(hs_valid, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 221, now there are 56 samples.

pbmc_samples <- subset_expt(hs_valid, subset="typeofcells=='pbmcs'")

## subset_expt(): There were 221, now there are 6 samples.

## Currently, these are not used, but instead I pulled the samples from the hs_clinical
## which means the biopsies are not included.
v1_samples <- subset_expt(hs_valid, subset="visitnumber=='1'")

## subset_expt(): There were 221, now there are 88 samples.

v1_biopsies <- subset_expt(v1_samples, subset="typeofcells=='biopsy'")

## subset_expt(): There were 88, now there are 28 samples.

v1_monocytes <- subset_expt(v1_samples, subset="typeofcells=='monocytes'")

## subset_expt(): There were 88, now there are 21 samples.

v1_neutrophils <- subset_expt(v1_samples, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 88, now there are 21 samples.

v1_eosinophils <- subset_expt(v1_samples, subset="typeofcells=='eosinophils'")

## subset_expt(): There were 88, now there are 12 samples.

v2_samples <- subset_expt(hs_valid, subset="visitnumber=='2'")

## subset_expt(): There were 221, now there are 60 samples.

v2_biopsies <- subset_expt(v2_samples, subset="typeofcells=='biopsy'")

## subset_expt(): There were 60, now there are 12 samples.

v2_monocytes <- subset_expt(v2_samples, subset="typeofcells=='monocytes'")

## subset_expt(): There were 60, now there are 19 samples.

v2_neutrophils <- subset_expt(v2_samples, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 60, now there are 18 samples.

v2_eosinophils <- subset_expt(v2_samples, subset="typeofcells=='eosinophils'")

## subset_expt(): There were 60, now there are 11 samples.

v3_samples <- subset_expt(hs_valid, subset="visitnumber=='3'")

## subset_expt(): There were 221, now there are 59 samples.

v3_biopsies <- subset_expt(v3_samples, subset="typeofcells=='biopsy'")

## subset_expt(): There were 59, now there are 12 samples.

v3_monocytes <- subset_expt(v3_samples, subset="typeofcells=='monocytes'")

## subset_expt(): There were 59, now there are 18 samples.

v3_neutrophils <- subset_expt(v3_samples, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 59, now there are 17 samples.

v3_eosinophils <- subset_expt(v3_samples, subset="typeofcells=='eosinophils'")

## subset_expt(): There were 59, now there are 12 samples.

6 Contrasts of interest

This might be a bit early to consider the contrasts, but I think we should consider this question immediately. The main thing we will be comparing is of course cure vs. fail; but we may also look at the visits and compare cell types.

cf_contrasts <- list(
    "fail_vs_cure" = c("failure", "cure"))
visit_contrasts <- list(
    "v2v1" = c("c2", "c1"),
    "v3v1" = c("c3", "c1"),
    "v3v2" = c("c3", "c2"))
type_contrasts <- list(
    "mono_biopsy" = c("monocytes", "biopsy"),
    "eosinophil_biopsy" = c("eosinophils", "biopsy"),
    "neutrophil_biopsy" = c("neutrophils", "biopsy"))

7 Parasite reads

Let us see if we can make an expressionset of the parasite reads in the TMRC3 samples and distinguish between the faux and real reads. E.g: Are there samples which definitively contain parasites?

lp_expt <- create_expt("sample_sheets/tmrc3_samples_20211207.xlsx",
                       file_column="lpanamensisv36hisatfile", gene_info = NULL) %>%
  subset_expt(coverage=20) %>%
  set_expt_conditions(fact="typeofcells")

## Reading the sample metadata.

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

## The sample definitions comprises: 238 rows(samples) and 75 columns(metadata fields).

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

## Matched 8778 annotations and counts.

## Bringing together the count matrix and gene information.

## Warning in create_expt("sample_sheets/tmrc3_samples_20211207.xlsx", file_column
## = "lpanamensisv36hisatfile", : The following samples have no counts!
## TMRC30010TMRC30144TMRC30145TMRC30146TMRC30185TMRC30207TMRC30217TMRC30219TMRC30220TMRC30209TMRC30210TMRC30211TMRC30212TMRC30213TMRC30214TMRC30216TMRC30221TMRC30222TMRC30223TMRC30225TMRC30226TMRC30227TMRC30229TMRC30230TMRC30231TMRC30233TMRC30234TMRC30235

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 8778 rows and 223 columns.

## The samples removed (and read coverage) when filtering samples with less than 20 reads are:

## TMRC30001 TMRC30002 TMRC30003 TMRC30004 TMRC30007 TMRC30008 TMRC30009 TMRC30010 
##        12         8         9        16         3        16        16         0 
## TMRC30011 TMRC30012 TMRC30013 TMRC30014 TMRC30030 TMRC30031 TMRC30037 TMRC30036 
##         5         9        13         4         3         8         9        11 
## TMRC30142 TMRC30143 TMRC30144 TMRC30145 TMRC30146 TMRC30147 TMRC30185 TMRC30124 
##         4         1         0         0         0         1         0         3 
## TMRC30207 TMRC30208 TMRC30217 TMRC30218 TMRC30219 TMRC30220 TMRC30209 TMRC30210 
##         0         1         0         2         0         0         0         0 
## TMRC30211 TMRC30212 TMRC30213 TMRC30214 TMRC30215 TMRC30216 TMRC30221 TMRC30222 
##         0         0         0         0         1         0         0         0 
## TMRC30223 TMRC30224 TMRC30225 TMRC30226 TMRC30227 TMRC30229 TMRC30230 TMRC30231 
##         0         2         0         0         0         0         0         0 
## TMRC30232 TMRC30233 TMRC30234 TMRC30235 
##         1         0         0         0

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

visit_fact <- pData(lp_expt)[["visitnumber"]]
batch_na <- is.na(visit_fact)
visit_fact[batch_na] <- "undefined"
lp_expt <- set_expt_batches(lp_expt, fact = visit_fact)

lp_norm <- normalize_expt(lp_expt, filter="simple", norm="quant",
                          convert="cpm", transform="log2")

## Removing 75 low-count genes (8703 remaining).

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

plotted <- plot_pca(lp_norm, plot_labels=FALSE)
plotted$plot

plotted_3d <- plot_3d_pca(plotted)

8 Distributions/Visualizations of interest

The sets of samples used to visualize the data will also comprise the sets used when later performing the various differential expression analyses.

8.1 Global metrics

Start out with some initial metrics of all samples. The most obvious are plots of the numbers of non-zero genes observed, heatmaps showing the relative relationships among the samples, the relative library sizes, and some PCA. It might be smart to split the library sizes up across subsets of the data, because they have expanded too far to see well on a computer screen.

The most likely factors to query when considering the entire dataset are cure/fail, visit, and cell type. This is the level at which we will choose samples to exclude from future analyses.

plot_legend(biopsy_samples)$plot

plot_libsize(biopsy_samples)$plot

plot_nonzero(biopsy_samples)$plot

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

## Minimum number of biopsy genes: ~ 14,000

plot_libsize(eosinophil_samples)$plot

plot_nonzero(eosinophil_samples)$plot

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

## Minimum number of eosinophil genes: ~ 13,500

plot_libsize(macrophage_samples)$plot

plot_nonzero(macrophage_samples)$plot

## Minimum number of macrophage genes: ~ 14,000

plot_libsize(monocyte_samples)$plot

plot_nonzero(monocyte_samples)$plot

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

## Minimum number of monocyte genes: ~ 7,500 before setting the minimum.

plot_libsize(neutrophil_samples)$plot

plot_nonzero(neutrophil_samples)$plot

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

## Minimum number of neutrophil genes: ~ 10,000 before setting minimum coverage.

plot_libsize(pbmc_samples)$plot

plot_nonzero(pbmc_samples)$plot

## Minimum number of pbmc genes: ~ 14,800

Now that those ‘global’ metrics are out of the way, lets look at some global metrics of the data following normalization; the most likely plots are of course PCA but also a couple of heatmaps.

type_valid <- set_expt_conditions(hs_valid, fact="typeofcells")
all_norm <- sm(normalize_expt(type_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

Continue looking, but switch the conditions/colors so that the clinical outcome becomes the focus and get rid of a few samples which are not actually a part of the TMRC3 focus (e.g. the PBMC and macrophage samples, which are all from the Wellcome Trust).

8.1.1 Clinically relevant samples

Included in this group will be the samples from patients who were 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 221, now there are 201 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: null.

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)

table(pData(hs_clinical)[["condition"]])

## 
##          cure       failure          lost notapplicable 
##           106            78            15             2

8.1.1.1 Remove the lost samples

Now let us focus purely on the cures and fails.

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

## subset_expt(): There were 201, now there are 184 samples.

The biopsy samples are a bit peculiar, therefore let us repeat this without them.

8.1.2 Repeat without biopsy samples

hs_clinical_nobiop <- hs_clinical %>%
  subset_expt(subset="typeofcells!='biopsy'") %>%
  subset_expt(subset="condition=='lost'|condition=='cure'|condition=='failure'")

## subset_expt(): There were 201, now there are 149 samples.

## subset_expt(): There were 149, now there are 147 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)

8.1.3 Visualize again with sva

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

## Removing 5207 low-count genes (14734 remaining).

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

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

## Setting 30774 low elements to zero.

## transform_counts: Found 30774 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

8.1.4 Samples separated by visit

Separate the samples by visit in order to more easily see what patterns emerge across cell type and clinical outcome.

I have a couple of likely starting points for this. The data sets: hs_clinical and clinical_nolost are the most likely. Given that this, at least in theory, the lost samples are not relevant.

8.1.4.1 All visits together

visit_expt <- set_expt_conditions(hs_clinical, fact = "visitnumber") %>%
  set_expt_batches(fact = "clinicaloutcome")

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

## Removing 5207 low-count genes (14734 remaining).

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

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

## Setting 35107 low elements to zero.

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

plot_pca(visit_nb)$plot

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

8.1.4.2 All visits, separated by cell type

visit_biopsy <- subset_expt(visit_expt, subset = "typeofcells=='biopsy'")

## subset_expt(): There were 201, now there are 52 samples.

visit_monocyte <- subset_expt(visit_expt, subset = "typeofcells=='monocytes'")

## subset_expt(): There were 201, now there are 58 samples.

visit_neutrophil <- subset_expt(visit_expt, subset = "typeofcells=='neutrophils'")

## subset_expt(): There were 201, now there are 56 samples.

visit_eosinophil <- subset_expt(visit_expt, subset = "typeofcells=='eosinophils'")

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

8.1.4.3 Visit 1 alone

v1_clinical <- subset_expt(visit_expt, subset = "visitnumber=='1'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells")

## subset_expt(): There were 201, now there are 82 samples.

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

## Removing 5627 low-count genes (14314 remaining).

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

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

## Setting 10928 low elements to zero.

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

plot_pca(v1_nb)$plot

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

8.1.4.4 Visit 2 alone

v2_clinical <- subset_expt(visit_expt, subset="visitnumber=='2'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells")

## subset_expt(): There were 201, now there are 60 samples.

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

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

## Removing 5839 low-count genes (14102 remaining).

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

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

## Setting 6254 low elements to zero.

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

plot_pca(v2_nb)$plot

8.1.4.5 Visit 3 alone

v3_clinical <- subset_expt(visit_expt, subset="visitnumber=='3'") %>%
  set_expt_conditions(fact = "clinicaloutcome") %>%
  set_expt_batches(fact = "typeofcells")

## subset_expt(): There were 201, now there are 59 samples.

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

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

## Removing 5774 low-count genes (14167 remaining).

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

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

## Setting 2511 low elements to zero.

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

plot_pca(v3_nb)$plot

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

8.1.5 Samples separated by cell type

Separate the samples by cell type in order to more easily observe patterns with respect to visit and clinical outcome.

8.1.5.1 Biopsies

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

## Removing 5738 low-count genes (14203 remaining).

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

plot_pca(biopsy_norm)$plot

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

## Removing 5738 low-count genes (14203 remaining).

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

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

## Setting 1876 low elements to zero.

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

plot_pca(biopsy_nb)$plot

8.1.5.2 Monocytes

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

## Removing 8744 low-count genes (11197 remaining).

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

plot_pca(monocyte_norm, plot_labels = FALSE)$plot

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

## Removing 8744 low-count genes (11197 remaining).

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

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

## Setting 1470 low elements to zero.

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

plot_pca(monocyte_nb, plot_labels = FALSE)$plot

8.1.5.3 Neutrophils

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

## Removing 10575 low-count genes (9366 remaining).

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

plot_pca(neutrophil_norm, plot_labels = FALSE)$plot

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

## Removing 10575 low-count genes (9366 remaining).

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

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

## Setting 1336 low elements to zero.

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

plot_pca(neutrophil_nb, plot_labels = FALSE)$plot

8.1.5.4 Eosinophils

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

## Removing 9113 low-count genes (10828 remaining).

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

plot_pca(eosinophil_norm, plot_labels = FALSE)$plot

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

## Removing 9113 low-count genes (10828 remaining).

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

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

## Setting 689 low elements to zero.

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

plot_pca(eosinophil_nb, plot_labels = FALSE)$plot

8.1.5.5 Monocytes, Neutrophils, and Eosinophils

9 Differential expression analyses

The primary goal is to learn about cure vs. fail.

9.1 Cure/Fail all samples

cf_clinical_de <- all_pairwise(hs_clinical, model_batch = "svaseq", filter = TRUE)

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

## 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 (14734 remaining).

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

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

## Setting 30774 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_clinical_tables <- combine_de_tables(
    cf_clinical_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_clinical_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/cf_clinical_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_clinical_tables-v202011.xlsx before writing the tables.

## Deleting the file excel/cf_clinical_sig-v202011.xlsx before writing the tables.

9.1.1 By cell type

9.1.1.1 Cure/Fail Biopsies

cf_biopsy_de <- all_pairwise(biopsy_samples, model_batch = "svaseq", filter = TRUE)

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

## 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 (14203 remaining).

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

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

## Setting 1876 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_biopsy_tables <- combine_de_tables(
    cf_biopsy_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_biopsy_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/cf_biopsy_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_biopsy_tables-v202011.xlsx before writing the tables.

9.1.1.2 Cure/Fail Monocytes

cf_monocyte_de <- all_pairwise(monocyte_samples, model_batch = "svaseq", filter = TRUE)

## batch_counts: Before batch/surrogate estimation, 2764 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 (11197 remaining).

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

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

## Setting 1470 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_monocyte_tables <- combine_de_tables(
    cf_monocyte_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_monocyte_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/cf_monocyte_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_monocyte_tables-v202011.xlsx before writing the tables.

## Deleting the file excel/cf_monocyte_sig-v202011.xlsx before writing the tables.

9.1.1.3 Cure/Fail Neutrophils

cf_neutrophil_de <- all_pairwise(neutrophil_samples, model_batch = "svaseq", filter = TRUE)

## batch_counts: Before batch/surrogate estimation, 1959 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 (9366 remaining).

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

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

## Setting 1336 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_neutrophil_tables <- combine_de_tables(
    cf_neutrophil_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_neutrophil_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/cf_neutrophil_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_neutrophil_tables-v202011.xlsx before writing the tables.

## Deleting the file excel/cf_neutrophil_sig-v202011.xlsx before writing the tables.

9.1.1.4 Cure/Fail Eosinophils

cf_eosinophil_de <- all_pairwise(eosinophil_samples, model_batch = "svaseq", filter = TRUE)

## batch_counts: Before batch/surrogate estimation, 973 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 (10828 remaining).

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

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

## Setting 689 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_eosinophil_tables <- combine_de_tables(
    cf_eosinophil_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_eosinophil_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/cf_eosinophil_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_eosinophil_tables-v202011.xlsx before writing the tables.

## Deleting the file excel/cf_eosinophil_sig-v202011.xlsx before writing the tables.

9.1.1.5 Cure/Fail clinical (Not biopsies)

cf_nobiopsy_de <- all_pairwise(hs_clinical_nobiop, model_batch = "svaseq", filter = TRUE)

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

## 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 (12168 remaining).

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

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

## Setting 14152 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_nobiopsy_tables <- combine_de_tables(
    cf_nobiopsy_de,
    keepers = cf_contrasts,
    excel = glue::glue("excel/cf_nobiopsy_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/cf_nobiopsy_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_nobiopsy_tables-v202011.xlsx before writing the tables.

## Deleting the file excel/cf_nobiopsy_sig-v202011.xlsx before writing the tables.

9.1.2 By visit

For these contrasts, we want to see fail_v1 vs. cure_v1, fail_v2 vs. cure_v2 etc. As a result, we will need to juggle the data slightly and add another set of contrasts.

9.1.3 Setting up

visit_cf_contrasts <- list(
    "v1fail_vs_cure" = c("v1failure", "v1cure"),
    "v2fail_vs_cure" = c("v2failure", "v2cure"),
    "v3fail_vs_cure" = c("v3failure", "v3cure"))
visit_cf_expt_factor <- paste0("v", pData(hs_clinical)[["visitnumber"]],
                               pData(hs_clinical)[["condition"]])
visit_cf_expt <- set_expt_conditions(hs_clinical, fact = visit_cf_expt_factor)

visit_cf_biopsy <- subset_expt(visit_cf_expt, subset="typeofcells=='biopsy'")

## subset_expt(): There were 201, now there are 52 samples.

visit_cf_eosinophil <- subset_expt(visit_cf_expt, subset="typeofcells=='eosinophils'")

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

visit_cf_monocyte <- subset_expt(visit_cf_expt, subset="typeofcells=='monocytes'")

## subset_expt(): There were 201, now there are 58 samples.

visit_cf_neutrophil <- subset_expt(visit_cf_expt, subset="typeofcells=='neutrophils'")

## subset_expt(): There were 201, now there are 56 samples.

9.1.3.1 Cure/Fail by visits, all cell types

visit_cf_all_de <- all_pairwise(visit_cf_expt, model_batch = "svaseq", filter = TRUE)

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

## 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 (14734 remaining).

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

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

## Setting 32744 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_all_tables <- combine_de_tables(
    visit_cf_all_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_all_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_cf_all_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_all_tables-v202011.xlsx before writing the tables.

9.1.3.2 Cure/Fail by visit, Biopsies

visit_cf_biopsy_de <- all_pairwise(visit_cf_expt, model_batch = "svaseq", filter = TRUE)

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

## 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 (14734 remaining).

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

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

## Setting 32744 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_biopsy_tables <- combine_de_tables(
    visit_cf_biopsy_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_biopsy_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_cf_biopsy_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_biopsy_tables-v202011.xlsx before writing the tables.

9.1.3.3 Cure/Fail by visit, Monocytes

visit_cf_monocyte_de <- all_pairwise(visit_cf_expt, model_batch = "svaseq", filter = TRUE)

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

## 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 (14734 remaining).

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

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

## Setting 32744 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_monocyte_tables <- combine_de_tables(
    visit_cf_monocyte_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_monocyte_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_cf_monocyte_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_monocyte_tables-v202011.xlsx before writing the tables.

9.1.3.4 Cure/Fail by visit, Neutrophils

visit_cf_neutrophil_de <- all_pairwise(visit_cf_expt, model_batch = "svaseq", filter = TRUE)

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

## 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 (14734 remaining).

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

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

## Setting 32744 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_neutrophil_tables <- combine_de_tables(
    visit_cf_neutrophil_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_neutrophil_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_cf_neutrophil_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_neutrophil_tables-v202011.xlsx before writing the tables.

9.1.3.5 Cure/Fail by visit, Eosinophils

visit_cf_eosinophil_de <- all_pairwise(visit_cf_expt, model_batch = "svaseq", filter = TRUE)

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

## 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 (14734 remaining).

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

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

## Setting 32744 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_cf_eosinophil_tables <- combine_de_tables(
    visit_cf_eosinophil_de,
    keepers = visit_cf_contrasts,
    excel = glue::glue("excel/visit_cf_eosinophil_tables-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_cf_eosinophil_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_cf_eosinophil_tables-v202011.xlsx before writing the tables.

9.2 Persistence in visit 3

Having put some SL read mapping information in the sample sheet, Maria Adelaida added a new column using it with the putative persistence state on a per-sample basis. One question which arised from that: what differences are observable between the persistent yes vs. no samples on a per-cell-type basis among the visit 3 samples.

9.2.1 Setting up

First things first, create the datasets.

persistence_expt <- subset_expt(hs_clinical, subset = "persistence=='Y'|persistence=='N'") %>%
  subset_expt(subset = 'visitnumber==3') %>%
  set_expt_conditions(fact = 'persistence')

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

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

persistence_biopsy <- subset_expt(persistence_expt, subset = "typeofcells=='biopsy'")

## subset_expt(): There were 49, now there are 9 samples.

persistence_monocyte <- subset_expt(persistence_expt, subset = "typeofcells=='monocytes'")

## subset_expt(): There were 49, now there are 16 samples.

persistence_neutrophil <- subset_expt(persistence_expt, subset = "typeofcells=='neutrophils'")

## subset_expt(): There were 49, now there are 14 samples.

persistence_eosinophil <- subset_expt(persistence_expt, subset = "typeofcells=='eosinophils'")

## subset_expt(): There were 49, now there are 10 samples.

9.2.2 Take a look

See if there are any patterns which look usable.

## All
persistence_norm <- normalize_expt(persistence_expt, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 5891 low-count genes (14050 remaining).

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

plot_pca(persistence_norm)$plot

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

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

## Removing 5891 low-count genes (14050 remaining).

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

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

## Setting 4937 low elements to zero.

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

plot_pca(persistence_nb)$plot

## Biopsies
persistence_biopsy_norm <- normalize_expt(persistence_biopsy, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 6483 low-count genes (13458 remaining).

plot_pca(persistence_biopsy_norm)$plot

## Insufficient data

## Monocytes
persistence_monocyte_norm <- normalize_expt(persistence_monocyte, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 9403 low-count genes (10538 remaining).

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

plot_pca(persistence_monocyte_norm)$plot

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

## Removing 9403 low-count genes (10538 remaining).

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

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

## Setting 106 low elements to zero.

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

plot_pca(persistence_monocyte_nb)$plot

## Neutrophils
persistence_neutrophil_norm <- normalize_expt(persistence_neutrophil, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 11381 low-count genes (8560 remaining).

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

plot_pca(persistence_neutrophil_norm)$plot

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

## Removing 11381 low-count genes (8560 remaining).

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

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

## Setting 136 low elements to zero.

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

plot_pca(persistence_neutrophil_nb)$plot

## Eosinophils
persistence_eosinophil_norm <- normalize_expt(persistence_eosinophil, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 9766 low-count genes (10175 remaining).

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

plot_pca(persistence_eosinophil_norm)$plot

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

## Removing 9766 low-count genes (10175 remaining).

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

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

## Setting 60 low elements to zero.

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

plot_pca(persistence_eosinophil_nb)$plot

I think this is a bust. But, I can do the de anyhow and see what happens…

9.2.3 persistence DE

persistence_de <- all_pairwise(persistence_expt, filter = TRUE, model_batch = "svaseq")

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

## 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 (14050 remaining).

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

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

## Setting 4937 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

persistence_table <- combine_de_tables(
    persistence_de,
    excel = glue::glue("excel/persistence_all_de-v{ver}.xlsx"))

## Deleting the file excel/persistence_all_de-v202011.xlsx before writing the tables.

persistence_monocyte_de <- all_pairwise(persistence_monocyte, filter = TRUE, model_batch = "svaseq")

## batch_counts: Before batch/surrogate estimation, 92 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 (10538 remaining).

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

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

## Setting 106 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/persistence_monocyte_de-v202011.xlsx before writing the tables.

persistence_neutrophil_de <- all_pairwise(persistence_neutrophil, filter = TRUE, model_batch = "svaseq")

## batch_counts: Before batch/surrogate estimation, 92 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 (8560 remaining).

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

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

## Setting 136 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/persistence_neutrophil_de-v202011.xlsx before writing the tables.

persistence_eosinophil_de <- all_pairwise(persistence_eosinophil, filter = TRUE, model_batch = "svaseq")

## batch_counts: Before batch/surrogate estimation, 88 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 (10175 remaining).

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

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

## Setting 60 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

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

## Deleting the file excel/persistence_eosinophil_de-v202011.xlsx before writing the tables.

9.3 Comparing visits without regard to cure/fail

9.3.1 All cell types

visit_all_de <- all_pairwise(visit_expt, filter = TRUE, model_batch = "svaseq")

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

## 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 (14734 remaining).

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

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

## Setting 35107 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_all_table <- combine_de_tables(
    visit_all_de,
    excel = glue::glue("excel/visit_all_de-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_all_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_all_de-v202011.xlsx before writing the tables.

## Deleting the file excel/visit_all_sig-v202011.xlsx before writing the tables.

9.3.2 Biopsy samples

visit_biopsy_de <- all_pairwise(visit_biopsy, filter = TRUE, model_batch = "svaseq")

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

## 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 (14203 remaining).

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

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

## Setting 1749 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_biopsy_table <- combine_de_tables(
    visit_biopsy_de,
    keepers = visit_contrasts,
    excel = glue::glue("excel/visit_biopsy_de-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_biopsy_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_biopsy_de-v202011.xlsx before writing the tables.

## Deleting the file excel/visit_biopsy_sig-v202011.xlsx before writing the tables.

9.3.3 Monocyte samples

visit_monocyte_de <- all_pairwise(visit_monocyte, filter = TRUE, model_batch = "svaseq")

## batch_counts: Before batch/surrogate estimation, 2764 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 (11197 remaining).

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

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

## Setting 1162 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_monocyte_table <- combine_de_tables(
    visit_monocyte_de,
    keepers = visit_contrasts,
    excel = glue::glue("excel/visit_monocyte_de-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_monocyte_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_monocyte_de-v202011.xlsx before writing the tables.

## Deleting the file excel/visit_monocyte_sig-v202011.xlsx before writing the tables.

9.3.4 Neutrophil samples

visit_neutrophil_de <- all_pairwise(visit_neutrophil, filter = TRUE, model_batch = "svaseq")

## batch_counts: Before batch/surrogate estimation, 1959 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 (9366 remaining).

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

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

## Setting 992 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_neutrophil_table <- combine_de_tables(
    visit_neutrophil_de,
    keepers = visit_contrasts,
    excel = glue::glue("excel/visit_neutrophil_de-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_neutrophil_sig-v{ver}.xlsx"))

## Deleting the file excel/visit_neutrophil_de-v202011.xlsx before writing the tables.

9.3.5 Eosinophil samples

visit_eosinophil_de <- all_pairwise(visit_eosinophil, filter = TRUE, model_batch = "svaseq")

## batch_counts: Before batch/surrogate estimation, 973 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 (10828 remaining).

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

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

## Setting 490 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

visit_eosinophil_table <- combine_de_tables(
    visit_eosinophil_de,
    keepers = visit_contrasts,
    excel = glue::glue("excel/visit_eosinophil_de-v{ver}.xlsx"),
    sig_excel = glue::glue("excel/visit_eosinophil_sig-v{ver}.xlsx"))

10 Likelihood Ratio Testing

10.1 Shared patterns across visits

lrt_visit <- deseq_lrt(hs_clinical, transform = "vst",
                       interactor_column = "visitnumber",
                       interest_column = "clinicaloutcome")

## converting counts to integer mode

## Error in checkFullRank(modelMatrix): the model matrix is not full rank, so the model cannot be fit as specified.
##   Levels or combinations of levels without any samples have resulted in
##   column(s) of zeros in the model matrix.
## 
##   Please read the vignette section 'Model matrix not full rank':
## 
##   vignette('DESeq2')

10.2 Shared patterns across cell types

11 Gene Set Analyses

11.1 GSEA

11.1.1 Cure/Fail groups

11.1.2 Cell type groups

11.1.3 Visit groups

11.2 GSVA