TMRC3 202411: Likelihood Ratio Tests and GSVA

1 Introduction

Though there is no fundamental reason to separate the likelihood ratio testing and GSVA from the differential expression analyses, they are both tasks which reach back to the primary expression data, unlike GSEA. Therefore I chose to separate these two tasks into a their own document.

2 LRT

EdgeR (McCarthy, Chen, and Smyth (2012)) and DESeq2 (Love, Huber, and Anders (2014)) both provide easily accessible LRT analyses. If I remember correctly, the function I wrote to simplify these analyses is aware of both, but this document primarily uses the results from DESeq2. I am reasonably certain they end up basically identical…

The likelihood ratio testing performed on the TMRC3 data is intended to look for shared patterns of expression across some known, constant factor. This may be time (are there patterns across the visits of each person), or cell type (do some genes act consistently across type). We could also ask this question of our two clinics, or indeed across the cure/fail samples.

2.1 Patterns across clinic

tc_clinical_filt <- normalize_expt(tc_clinical, filter = TRUE)

## Removing 5654 low-count genes (14298 remaining).

tc_lrt_clinic <- deseq_lrt(tc_clinical_filt, transform = "vst", interaction = FALSE,
                           interactor_column = "visitnumber",
                           interest_column = "clinic")

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 169 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## A large number of genes was given-- please, make sure this is not an error. Normally, only DE genes will be useful for this function.

## Working with 6279 genes.

## Working with 6279 genes after filtering: minc > 3

## Joining with `by = join_by(merge)`
## Joining with `by = join_by(merge)`

## Warning in `labels<-.dendrogram`(dend, value = value, ...): The lengths of the
## new labels is shorter than the number of leaves in the dendrogram - labels are
## recycled.

tc_lrt_clinic[["cluster_data"]][["plot"]]

2.2 Patterns across visits, only Tumaco

2.2.1 All clinical cell types

t_clinical_filt <- normalize_expt(t_clinical, filter=TRUE)

## Removing 5796 low-count genes (14156 remaining).

lrt_visit <- deseq_lrt(t_clinical_filt, transform = "vst", interaction = FALSE,
                       interactor_column = "visitnumber",
                       interest_column = "finaloutcome")

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 120 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## A large number of genes was given-- please, make sure this is not an error. Normally, only DE genes will be useful for this function.

## Working with 5445 genes.

## Working with 5444 genes after filtering: minc > 3

## Joining with `by = join_by(merge)`
## Joining with `by = join_by(merge)`

## Warning in `labels<-.dendrogram`(dend, value = value, ...): The lengths of the
## new labels is shorter than the number of leaves in the dendrogram - labels are
## recycled.

lrt_visit$cluster_data$plot

summary(lrt_visit[["favorite_genes"]])

##     genes              cluster     
##  Length:5444        Min.   : 1.00  
##  Class :character   1st Qu.: 2.00  
##  Mode  :character   Median : 3.00  
##                     Mean   : 3.65  
##                     3rd Qu.: 4.00  
##                     Max.   :11.00

written <- write_xlsx(data = as.data.frame(lrt_visit[["deseq_table"]]),
                      excel = glue("excel/lrt_clinical_visit-v{ver}.xlsx"))

2.2.2 Monocytes, only Tumaco

lrt_monocyte_visit <- deseq_lrt(t_visitcf_monocyte, transform = "vst",
                                interaction = FALSE,
                                interactor_column = "visitnumber",
                                interest_column = "finaloutcome")

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 68 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## Working with 12 genes.

## Working with 12 genes after filtering: minc > 3

## Joining with `by = join_by(merge)`
## Joining with `by = join_by(merge)`

## Warning in `labels<-.dendrogram`(dend, value = value, ...): The lengths of the
## new labels is shorter than the number of leaves in the dendrogram - labels are
## recycled.

lrt_monocyte_visit[["cluster_data"]][["plot"]]

2.2.3 Neutrophils, only Tumaco

lrt_neutrophil_visit <- deseq_lrt(t_visitcf_neutrophil, transform = "vst",
                                  interaction = FALSE,
                                  interactor_column = "visitnumber",
                                  interest_column = "finaloutcome")

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 49 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## Working with 952 genes.

## Working with 952 genes after filtering: minc > 3

## Joining with `by = join_by(merge)`
## Joining with `by = join_by(merge)`

## Warning in `labels<-.dendrogram`(dend, value = value, ...): The lengths of the
## new labels is shorter than the number of leaves in the dendrogram - labels are
## recycled.

lrt_neutrophil_visit[["cluster_data"]][["plot"]]

2.2.4 Eosinophils, only Tumaco

lrt_eosinophil_visit <- deseq_lrt(t_visitcf_eosinophil, transform = "vst",
                                  interaction = FALSE,
                                  interactor_column = "visitnumber",
                                  interest_column = "finaloutcome")

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## Warning in deseq_lrt(t_visitcf_eosinophil, transform = "vst", interaction =
## FALSE, : There are no significant differences given the 0.05 adjusted p-value.

## Returning the full LRT table just so that you have something to look at.

2.3 Shared patterns across cell types

lrt_celltype_clinical_test <- deseq_lrt(tc_clinical, transform = "vst",
                                        interactor_column = "typeofcells",
                                        interest_column = "finaloutcome")

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 15 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## Working with 535 genes.

## Working with 534 genes after filtering: minc > 3

## Joining with `by = join_by(merge)`
## Joining with `by = join_by(merge)`

## Warning in `labels<-.dendrogram`(dend, value = value, ...): The lengths of the
## new labels is shorter than the number of leaves in the dendrogram - labels are
## recycled.

hs_annot <- fData(hs_expt)
deseq_lrt_df <- merge(hs_annot, as.data.frame(lrt_celltype_clinical_test[["deseq_table"]]), all.y=TRUE,
                      by="row.names")
rownames(deseq_lrt_df) <- deseq_lrt_df[["Row.names"]]
deseq_lrt_df[["Row.names"]] <- NULL
written <- write_xlsx(data=deseq_lrt_df,
                      excel=glue("excel/lrt_clinical_celltype-v{ver}.xlsx"))

2.4 Using LRT to examine donor effects

Suggestion from Neal: Get residual deviance of donor+visit, visit, donor.

Oh oh, we are actually comparing the results of these models: Define terms: exprs: expression matrix with rows=genes, column=samples This is synonymous with ‘outcome’. Conversely, finaloutcome is the term we wish to explain via the metadata and expression matrix.

• exprs ~ visit, data = meta+exprs
• exprs ~ donor, data = meta+exprs (this does not support the argument)
• exprs ~ donor + visit (caveat regarding donors with <=2, ergo will need to restrict sample set for other contrasts as well.)

1. Extract from deseq/limma/edger the residual deviance folowing the fit of each model. In the case of DESeq2, this may be acquired via a LRT test.
2. This is expected to return the deviance between the various models applied.
3. ok I think I get it, invoke deseq_lrt with donor and visit+donor and examine the results.

2.4.1 Setting up

I tested a few models and it seems that the LRT test accepts a model of ~visit+donor.

lrt_eo_vd <- deseq_lrt(t_eosinophils, transform = "vst",
                       interactor_column = "visitnumber",
                       interest_column = "donor", interaction = FALSE)

## Warning in deseq_lrt(t_eosinophils, transform = "vst", interactor_column =
## "visitnumber", : The donor should probably be a factor, set it with the
## 'factors' arg.

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## Working with 44 genes.

## Working with 39 genes after filtering: minc > 3

## Joining with `by = join_by(merge)`
## Joining with `by = join_by(merge)`

## Warning: Failed to fit group 2.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 3.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 5.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 10.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 11.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 2.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 3.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 5.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 10.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 11.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 2.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 3.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 5.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 10.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 11.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 2.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 3.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 5.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 10.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning: Failed to fit group 11.
## Caused by error in `poly()`:
## ! 'degree' must be less than number of unique points

## Warning in `labels<-.dendrogram`(dend, value = value, ...): The lengths of the
## new labels is shorter than the number of leaves in the dendrogram - labels are
## recycled.

## Oh wait, I think that is the opposite of our goal, lets do the other way around.

## This creates a dataset from ~ donor + visitnumber
## and performs the LRT with ~ visitnumber
lrt_eo_dv <- deseq_lrt(t_eosinophils, transform = "vst",
                       interactor_column = "donor",
                       interest_column = "visitnumber", interaction = FALSE)

## Warning in deseq_lrt(t_eosinophils, transform = "vst", interactor_column =
## "donor", : The donor should probably be a factor, set it with the 'factors'
## arg.

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## A large number of genes was given-- please, make sure this is not an error. Normally, only DE genes will be useful for this function.
## Working with 6086 genes.
## Working with 6028 genes after filtering: minc > 3
## Joining with `by = join_by(merge)`Joining with `by = join_by(merge)`

## Warning in `labels<-.dendrogram`(dend, value = value, ...): The lengths of the
## new labels is shorter than the number of leaves in the dendrogram - labels are
## recycled.

Let us consider the rank of the eosinophil data. We are hypothesizing that each donor must have >= 2 samples for ~ donor + finaloutcome to work. Thus I performed table() on the experimental design and see that donor 2052 is unsuitable.

3 GSVA

GSVA(Hänzelmann, Castelo, and Guinney (2013)) is a completely different family of analysis, and at least in my hands is used as a way to explore and look for publications which may be of interest. The general idea is that it performs a specific set of normalizations on the raw data, rank orders the results; then cross references them against extant analyses looking for over represented gene sets in specific papers or categories (e.g. reactome/GO/etc). This may either use a built-in set of mSigDB (Liberzon et al. (2011)) gene sets, or you may pull in manually downloaded (newer) data. On my workstation, I do the latter because I can leverage the downloaded gmt/xml/etc files in order to get the full annotations for the gene sets and accompanying papers. This container just uses the pre-downloaded mSigDB data (“GSVAdata” (n.d.)), which is a bit older.

Reminder: I disabled these for the moment in order to consider how to handle the (non)inclusion of the mSigDB gmt files.

3.1 Load some signatures

If one chooses, simple_gsva() can either load signatures from the GSVAdata package, which is a little old, or load an arbitrary set. load_gmt_signatures() provides a quick way to extract them from a gmt file.

broad_c7 <- load_gmt_signatures(signatures = "reference/msigdb/c7.all.v7.5.1.entrez.gmt",
                                signature_category = "c7")
broad_c2 <- load_gmt_signatures(signatures = "reference/msigdb/c2.all.v7.5.1.entrez.gmt",
                                signature_category = "c2")
broad_h <- load_gmt_signatures(signatures = "reference/msigdb/h.all.v7.5.1.entrez.gmt",
                               signature_category = "h")

Actually, I am going to try to semi-change my mind: GSVA comes with a dataset comprised of some older mSigDB data. Thus I will repeat each of the following blocks with an invocation which uses the older data and leave the previous invocations here so you can see what I actually ran to make pretty pictures. In the best case scenario, the results should be nearly identical with only a few categories missing between the new copy I downloaded and loaded above and the loaded-by-gsva version.

3.1.1 Clinical samples

3.1.1.1 Clinical C2

tc_celltype_gsva_c2 <- simple_gsva(
    tc_valid, signatures = broad_c2,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml")
tc_celltype_gsva_c2_sig <- get_sig_gsva_categories(
    tc_celltype_gsva_c2,
    excel = "analyses/3_cali_and_tumaco/GSVA/tc_valid_gsva_c2.xlsx")
tc_celltype_gsva_c2_sig$subset_plot
tc_celltype_gsva_c2_sig$score_plot

There are two other important differences between these results and the results obtained when using the above invocations and the following:

1. The above invocations were done before mSigDB changed the xml format so that it no longer successfully parses. I therefore changed my function to use the newer sqlite versions for versions of mSigDB > 7.2. Thus I should change the above to reflect this.
2. The R data package version of the mSigDB (at least the last time I checked) does not have all the fun/interesting experiment annotation information included in the xml/sqlite. Thus, the following results will be sparser (unless I am mistaken about the GSVAdata) vis a vis the annotations.

heh, only 1 way to find out! Let us try!

Reminder to self, these are the current mSigDB human categories, I am not certain which (I assume all) are included in GSVAdata.

• H: Hallmark genes, 50 gene sets for well defined states/processes.
• C1: Genes by chromosomal position. E.g. using the karyogram to make groups.
• C2: Curated gene sets: Various other datasets brought together: BioCarta pathways, KEGG_MEDICUS, PID pathways, reactome pathways, wikipathways, and legacy KEGG.
• C3: Regulatory targets sets (e.g. expected to be affected by stuff like miRNA): MIR targets, miRDB targets, MIR_LEGACY, TFT genes, and the GTRD subset of TFT.
• C4: Computational sets in cancer: 3CA from the cell atlas, GCN (cancer gene neighborhoods), CM: Cancer modules.
• C5: Gene ontology – this is split between the canonical MF/BP/CC and the human phenotype pathology dataset.
• C6: oncogenic signature gene sets: pulled from NCBI GEO and/or Broad experiments.
• C7: immunologic signature gene sets: Perturbations of the immune system: ImmuneSigDB, VAX (this one is our favorite for obvious (I think) reasons)
• C8: cell type signatures

tc_celltype_gsva_c2 <- simple_gsva(
  tc_valid, signature_category = "c2")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

tc_celltype_gsva_c2_sig <- get_sig_gsva_categories(
    tc_celltype_gsva_c2,
    excel = "analyses/3_cali_and_tumaco/GSVA/tc_valid_gsva_c2.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: failure_vs_cure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: cure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: failure.  Adjust = BH

## The factor cure has 122 rows.

## The factor failure has 62 rows.

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

tc_celltype_gsva_c2_sig$subset_plot

tc_celltype_gsva_c2_sig$score_plot

3.1.1.2 Valid C7

Ditto for C7

tc_celltype_gsva_c7 <- simple_gsva(
    tc_valid, signatures = broad_c7,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category="c7")
tc_celltype_gsva_c7_sig <- get_sig_gsva_categories(
    tc_celltype_gsva_c7,
    excel = "analyses/3_cali_and_tumaco/GSVA/tc_valid_gsva_c7.xlsx")

tc_celltype_gsva_c7 <- simple_gsva(
  tc_valid, signature_category = "c7")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

tc_celltype_gsva_c7_sig <- get_sig_gsva_categories(
    tc_celltype_gsva_c7,
    excel = "analyses/3_cali_and_tumaco/GSVA/tc_valid_gsva_c7.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: failure_vs_cure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: cure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: failure.  Adjust = BH

## The factor cure has 122 rows.

## The factor failure has 62 rows.

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

tc_celltype_gsva_c7_sig$subset_plot

tc_celltype_gsva_c7_sig$score_plot

I am not too fussed about the hallmark genes.

3.1.1.3 Valid H

tc_celltype_gsva_h <- simple_gsva(
    tc_valid,
    signatures = broad_h,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "h")
tc_celltype_gsva_h_sig <- get_sig_gsva_categories(
    tc_celltype_gsva_h,
    excel = "analyses/3_cali_and_tumaco/GSVA/tc_valid_gsva_h.xlsx")

3.1.2 Tumaco samples

3.1.2.1 Clinical C2

The C2 set includes a broad array of databases and studies. Thus when we look at the various gene sets which are deemed to be poor scoring, we see things like ‘the retinoid cycle in cones’ or ‘Nicotine metabolism’.

t_clinical_gsva_c2 <- simple_gsva(
    t_clinical,
    signatures = broad_c2,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c2")
t_clinical_gsva_c2_sig <- get_sig_gsva_categories(
    t_clinical_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_clinical_gsva_c2.xlsx")
t_clinical_gsva_c2_sig$subset_plot
t_clinical_gsva_c2_sig$score_plot

t_clinical_gsva_c2 <- simple_gsva(
    t_clinical, signature_category = "c2")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_clinical_gsva_c2_sig <- get_sig_gsva_categories(
    t_clinical_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_clinical_gsva_c2.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: failure_vs_cure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: cure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: failure.  Adjust = BH

## The factor cure has 67 rows.

## The factor failure has 56 rows.

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

t_clinical_gsva_c2_sig$subset_plot

t_clinical_gsva_c2_sig$score_plot

The get_sig_gsva_categories() function uses a few metrics to attempt to find msigDB categories of interest from the GSVA results. One of the methods employed which I hope will prove useful is to test the scores observed for one condition against all values in order to look for categories which are the most likely to be interesting.

With that in mind, here are the brief descriptions for a few of the most highly scored C2 groups and their associated PMIDs. Note that this are taken from the set of failed samples vs. all; there are some differences among the cure samples, but they are generally very similar.

• IL22 Soluble Receptor Signaling Pathway
• LRR FLII-interacting protein 1 (LRRFIP1) activates type I IFN production
• Stat3 Signaling Pathway
• DEx/H-box helicases activate type I IFN and inflammatory cytokines production
• Maturation of SARS-CoV-1 spike protein
• 18285459: Proteins significantly induced by oxidative stress (hydrogen peroxide [PubChem=784] in 786-O cells (renal clear cell carcinoma, RCC) expressing VHL [GeneID=7428].
• Negative feedback regulation of MAPK pathway
• 16424014: Genes from the 12p region that were up-regulated in choriocarcinoma cells compared to normal testis.
• Dicer Pathway
• 17440165: Proteins with reduced expression in mulignant glioma cell line (A172) which bears loss of heterozygosity (LOH) in the 1p region.
• Interleukin-6 signaling
• Human Cytomegalovirus and Map Kinase
• Signaling by cytosolic FGFR1 fusion mutants
• 15710396: Genes down-regulated during transition from G2 (moderately differentiated tumor, infected with HCV) to G3 (poorly differentiated tumor, infected with HCV) in the development of hepatocellular carcinoma.
• 15288478: Genes down-regulated in hepatocellular carcinoma (HCC) with early recurrence.
• Endosomal/Vacuolar pathway
• IkBA variant leads to EDA-ID
• Canonical NF-kB pathway
• 17234770: Genes involved in cell cycle regulation which were up-regulated in MCF-7 cells (breast cancer) by tretinoin (retinoic acid) [PubChem=444795].
• 15710396: Genes down-regulated during transition from L0 (non-tumor, not infected with HCV) to L1 (non-tumor, infected with HCV) in the development of hepatocellular carcinoma.
• Maturation of SARS-CoV-1 nucleoprotein
• IFN gamma signaling pathway
• Pilocytic astrocytoma
• 17072321: Genes encoding the NF-kB core signaling proteins.
• MAPK1 (ERK2) activation
• TNFR2 Signaling Pathway
• CLEC7A/inflammasome pathway
• Regulation of PTEN mRNA translation
• 16140955: Genes up-regulated synergistically in NB4 cells (acute promyelocytic leukemia, APL) by tretinoin and NSC682994 [PubChem=444795;388304].

3.1.2.2 Clinical C7

In contrast, the C7 set is relatively focused on the immune system, and is therefore likely to be of direct interest.

t_clinical_gsva_c7 <- simple_gsva(
    t_clinical,
    signatures = broad_c7,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c7")
t_clinical_gsva_c7_sig <- get_sig_gsva_categories(
    t_clinical_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_clinical_gsva_c7.xlsx")

t_clinical_gsva_c7 <- simple_gsva(
    t_clinical, signature_category = "c7")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_clinical_gsva_c7_sig <- get_sig_gsva_categories(
    t_clinical_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_clinical_gsva_c7.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: failure_vs_cure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: cure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: failure.  Adjust = BH

## The factor cure has 67 rows.

## The factor failure has 56 rows.

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

With the above in mind, here are a few of the most highly scored papers/sets when scoring the failed samples:

• 21357945: Genes positively correlated with titer response index in peripheral blood mononuclear cell in Caucasian male adults (18-40) (high responders) after exposure to Fluarix/Fluvirin , time point 1D and 3DY. Comment: Signature predictive of titer response index (TRI). Day 1 and day 3 values averaged.
• 27764254: Genes negatively correlated with T cell responses (long term) in peripheral blood mononuclear cell in seniors (50-75) after exposure to Zostavax , time point 1D. Comment: (B) Network of genes informative of long-term responses.
• 29535712: Genes down-regulated in blood 1d vs 0hr in adults (18-45) after exposure to CN54gp140 adjuvanted with GLA-AF , time point 1D , administered i.m.
• 15789058: Genes down-regulated in comparison of naive B cells versus day 0 monocytes.
• 15789058: Genes down-regulated in comparison of naive CD8 T cells versus day 0 monocytes.
• 15789058: Genes down-regulated in comparison of dendritic cells (DC) versus monocytes.
• 15789058: Genes up-regulated in comparison of monocytes cultured for 0 days versus those cultured for 7 days.
• 21743478: Genes up-regulated in comparison of monocytes versus plasmacytoid dendritic cells (pDC).
• 27764254: Genes positively correlated with expansion of VZV specific T cells (0d to peak) in peripheral blood mononuclear cell in seniors (50-75) after exposure to Zostavax , time point 1D
• 15789058: Genes down-regulated in comparison of naive CD4 [GeneID=920] T cells versus day 0 monocytes.
• 24808365: Genes up-regulated in dendritic cells infected by Leishmania major: 2h versus 24h.
• 18275831: Genes down-regulated in comparison of systemic lupus erythematosus B cells versus systemic lupus erythromatosus myeloid cells.
• 21743478: Genes up-regulated in comparison of monocytes from influenza vaccinee at day 7 post-vaccination versus plasmacytoid dendritic cells (mDC) at day 7 post-vaccination.
• 15789058: Genes down-regulated in comparison of naive CD4 [GeneID=920] CD8 T cells versus monocytes cultured for 0 days.
• 21743478: Genes down-regulated in peripheral blood mononuclear cell 7d vs 0d in adults (18-50) after exposure to FluMist , time point 7D. Comment: Molecular signature induced by LAIV vaccination. (a) Interferon (IFN)-related genes differentially expressed after LAIV vaccination
• 27764254: Genes positively correlated with contraction of VZV specific T cells (peak to 28d) in peripheral blood mononuclear cell in seniors (50-75) after exposure to Zostavax , time point 1D
• 25596819: Genes down-regulated in peripheral blood mononuclear cell 2d vs 0d in seniors (70+) (nonresponder) after exposure to Inactivated influenza vaccine , time point 2D
• 21743478: Genes up-regulated in comparison of monocytes from influenza vaccinee at day 7 post-vaccination versus myeloid dendritic cells at day 7 post-vaccination.
• 24808365: Genes up-regulated in dendritic cells: untreated versus 24h after infection of Leishmania major.
• 21743478: Genes down-regulated in comparison of B cells versus monocytes.
• 17105821: Genes down-regulated in comparison of peripheral blood mononuclear cells (PBMC) from healthy donors versus PBMC from patients with acute S. aureus infection.
• 21743478: Genes down-regulated in comparison of B cells from influenza vaccinee at day 7 versus monocytes from influenza vaccinee at day 7.
• 18275831: Genes down-regulated in comparison of systemic lupus erythematosus CD4 [GeneID=920] T cells versus systemic lupus erythematosus myeloid cells.
• 15789058: Genes down-regulated in comparison of naive B cells versus unstimulated neutrophils.
• 19047440: Genes up-regulated in peripheral blood mononuclear cell 28d vs 0d in unknown after exposure to YF-Vax/Stamaril , time point 28D
• 21093321: Genes up-regulated in bone marrow-derived macrophages with STAT6 [GeneID=6778] knockout treated with rosiglitazone [PubChem=77999]: control versus IL4 [GeneID=3565].
• 21743478: Genes up-regulated in peripheral blood mononuclear cell 3d vs 0d in adults (18-50) after exposure to FluMist , time point 3D. Comment: Molecular signature induced by LAIV vaccination. (a) Interferon (IFN)-related genes differentially expressed after LAIV vaccination
• 18292579: Genes down-regulated in comparison of monocytes treated with anti-TREM1 [GeneID=54210] versus monocytes treated with control IgG.
• 21743478: Genes down-regulated in comparison of plasmacytoid dendritic cells (DC) versus myeloid DCs.
• 22986450: Genes up-regulated in comparison of control polymorphonuclear leukocytes (PMN) at 12 h versus PMN treated with F. tularensis vaccine at 24 h.
• 28099485: Genes up-regulated in B cell 1d vs 0d in adults (18-49) after exposure to inactivated monovalent influenza A/Indonesia/05/2005 H5N1 split-virus vaccine , time point 1D , administered i.m.
• 18292579: Genes down-regulated in comparison of monocytes treated with anti-TREM1 [GeneID=54210] versus untreated monocytes.
• 17105821: Genes up-regulated in comparison of peripheral blood mononuclear cells (PBMC) from patients with acute influenza infection versus PBMC from patients with acute E. coli infection.
• 21093321: Genes up-regulated in bone marrow-derived macrophages with STAT6 [GeneID=6778] knockout treated with IL4 [GeneID=3565]: control versus rosiglitazone [PubChem=77999].
• 21743478: Genes up-regulated in comparison of monocytes versus myeloid dendritic cells (mDC).
• 21743478: Genes up-regulated in peripheral blood mononuclear cell 7d vs 0d in adults (18-50) after exposure to FluMist , time point 7D. Comment: Molecular signature induced by LAIV vaccination. (a) Interferon (IFN)-related genes differentially expressed after LAIV vaccination
• 17595242: Genes up-regulated in comparison of peripheral blood mononuclear cells (PBMC) from healthy donors versus PBMCs from patients with type 1 diabetes at 4 month after the diagnosis.
• 16474395: Genes up-regulated in comparison of neutrophils versus central memory CD4 [GeneID=920] T cells.
• 22986450: Genes up-regulated in comparison of control polymorphonuclear leukocytes (PMN) at 0 h versus PMN treated with F. tularensis vaccine at 24 h.
• 21743478: Genes down-regulated in comparison of plasmacytoid dendritic cells (DC) from influenza vaccinee at day 7 post-vaccination versus myeloid DCs at day 7 post-vaccination.

3.1.2.3 Clinical H

The H set is both much larger and smaller, it is comprised of just 50 sets, but they have many more genes associated with them.

t_clinical_gsva_h <- simple_gsva(
    t_clinical,
    signatures = broad_h,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category="h")
t_clinical_gsva_h_sig <- get_sig_gsva_categories(
    t_clinical_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_clinical_gsva_h.xlsx")

t_clinical_gsva_h <- simple_gsva(
    t_clinical, signature_category = "h")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_clinical_gsva_h_sig <- get_sig_gsva_categories(
    t_clinical_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_clinical_gsva_h.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: failure_vs_cure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: cure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: failure.  Adjust = BH

## The factor cure has 67 rows.

## The factor failure has 56 rows.

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

• Genes up-regulated in response to alpha interferon proteins.
• Genes involved in protein secretion pathway.
• Genes up-regulated in response to IFNG [GeneID=3458].
• A subgroup of genes regulated by MYC - version 1 (v1).
• Genes up-regulated by activation of the PI3K/AKT/mTOR pathway.
• Genes up-regulated during unfolded protein response, a cellular stress response related to the endoplasmic reticulum.
• Genes up-regulated in response to TGFB1 [GeneID=7040].
• Genes up-regulated by reactive oxigen species (ROS).
• Genes important for mitotic spindle assembly.
• Genes regulated by NF-kB in response to TNF [GeneID=7124].
• Genes up-regulated through activation of mTORC1 complex.
• Genes encoding proteins involved in oxidative phosphorylation.
• Genes defining response to androgens.
• Genes mediating programmed cell death (apoptosis) by activation of caspases.
• Genes up-regulated by IL6 [GeneID=3569] via STAT3 [GeneID=6774], e.g., during acute phase response.
• Genes involved in p53 pathways and networks.

3.1.2.4 Tumaco Biopsies c2

t_biopsy_gsva_c2 <- simple_gsva(
    t_biopsies,
    signatures = broad_c2,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category="c2")
t_biopsy_gsva_c2_sig <- get_sig_gsva_categories(
    t_biopsy_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_biopsy_gsva_c2.xlsx")

t_biopsy_gsva_c2 <- simple_gsva(
    t_biopsies, signature_category = "c2")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_biopsy_gsva_c2_sig <- get_sig_gsva_categories(
    t_biopsy_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_biopsy_gsva_c2.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 9 rows.

## The factor tumaco_failure has 5 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

• 17440165: Proteins with reduced expression in mulignant glioma cell line (A172) which bears loss of heterozygosity (LOH) in the 1p region.
• 15710396: Genes down-regulated during transition from L0 (non-tumor, not infected with HCV) to L1 (non-tumor, infected with HCV) in the development of hepatocellular carcinoma.
• Macrophage markers
• 18285459: Proteins significantly induced by oxidative stress (hydrogen peroxide [PubChem=784] in 786-O cells (renal clear cell carcinoma, RCC) expressing VHL [GeneID=7428].
• IL22 Soluble Receptor Signaling Pathway
• Maturation of SARS-CoV-1 spike protein
• 16424014: Genes from the 12p region that were up-regulated in choriocarcinoma cells compared to normal testis.
• Uptake and function of diphtheria toxin
• Lck and Fyn tyrosine kinases in initiation of TCR Activation
• Type I hemidesmosome assembly
• 12618007: Genes down-regulated after heat shock in peripheral lympocytes from old donors, compared to those from the young ones.
• Stat3 Signaling Pathway
• Nef and signal transduction
• T Helper Cell Surface Molecules
• LRR FLII-interacting protein 1 (LRRFIP1) activates type I IFN production
• Endosomal/Vacuolar pathway
• T Cytotoxic Cell Surface Molecules
• Folding of actin by CCT/TriC
• CTL mediated immune response against target cells
• PD-1 signaling
• DEx/H-box helicases activate type I IFN and inflammatory cytokines production
• 16140955: Genes up-regulated synergistically in NB4 cells (acute promyelocytic leukemia, APL) by tretinoin and NSC682994 [PubChem=444795;388304].
• 15710396: Genes down-regulated during transition from G2 (moderately differentiated tumor, infected with HCV) to G3 (poorly differentiated tumor, infected with HCV) in the development of hepatocellular carcinoma.
• Interleukin-6 signaling
• 18757430: Immune response genes up-regulated in zenograft tumors formed by SNU-601 cells (gastric cancer) made to express LRRC3B [GeneID=116135].
• TCA cycle nutrient use and invasiveness of ovarian cancer
• 18245477: Genes in the green cluster of protein kinases distinguishing between luminal A and basal breast cancer subtypes.
• 20308323: Proteins identified by mass spectrometry in complexes containing ALKBH8 [GeneID=91801].
• Cycling of Ran in nucleocytoplasmic transport
• Activation of the mRNA upon binding of the cap-binding complex and eIFs, and subsequent binding to 43S
• 18451147: Genes down-regulated in HSA/c and KYSE140 cells (esophageal squamous cell carcinoma, ESCC) after knockdown of PTTG1 [GeneID=9232] by RNAi.
• 17187432: Down-regulated genes in hepatocellular carcinoma (HCC) subclass G5, defined by unsupervised clustering.
• 15608688: Cluster 4: genes with similar expression profiles across follicular thyroid carcinoma (FTC) samples.
• Basic Mechanisms of SUMOylation
• IFN gamma signaling pathway
• 15531917: Selected down-regulated genes distinguishing between Wilms tumors of different histological types: anaplastic vs favorable histology.
• 15608688: Cluster 3: genes with similar expression profiles across follicular thyrorid carcinoma (FTC) samples; genes in this cluster correlated well with the presence of PAX8-PPARG [GeneID=7849;5468] fusion protein.
• Internal Ribosome entry pathway
• 11773596: Housekeeping genes identified as expressed across 19 normal tissues.
• 16909099: Genes down-regulated in adult T-cell leukemia (ATL), chronic vs acute clinical condition.
• 17486082: The ‘TEB profile genes’: down-regulated during pubertal mammary gland development specifically in the TEB (terminal end bud) structures.
• CLEC7A/inflammasome pathway
• Nef Mediated CD4 Down-regulation
• 15608688: Genes down-regulated in follicular thyroid carcinoma (FTC) samples that bear PAX8-PPARG [GeneID=7849;5468] fusion protein.
• SUMO is conjugated to E1 (UBA2:SAE1)
• 18724378: Genes up-regulated in MCF7-ADR cell line (breast cancer) resistant to docetaxel [PubChem=148124].
• MAPK1 (ERK2) activation
• 16760442: Down-regulated genes constituting the molecular signature of Burkitt ’s lymphoma.
• Antigen Processing and Presentation
• 15295046: Genes distinguishing asparaginase resistant and sensitive B-lineage ALL; here - genes up-regulated in the drug resistant samples.

3.1.2.5 Tumaco Biopsies c7

t_biopsy_gsva_c7 <- simple_gsva(
    t_biopsies,
    signatures = broad_c7,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c7")
t_biopsy_gsva_c7_sig <- get_sig_gsva_categories(
    t_biopsy_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_biopsy_gsva_c7.xlsx")

t_biopsy_gsva_c7 <- simple_gsva(
    t_biopsies, signature_category = "c7")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_biopsy_gsva_c7_sig <- get_sig_gsva_categories(
    t_biopsy_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_biopsy_gsva_c7.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 9 rows.

## The factor tumaco_failure has 5 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

• 28099485: Genes up-regulated in B cell 1d vs 0d in adults (18-49) after exposure to inactivated monovalent influenza A/Indonesia/05/2005 H5N1 split-virus vaccine , time point 1D , administered i.m.
• 24912498: Genes up-regulated in peripheral blood mononuclear cell 28d vs 7d in infants (4-6m) (BCG-primed) after exposure to Modified Vaccinia Ankara (MVA) virus vaccine vector , time point 28D
• 24912498: Genes up-regulated in peripheral blood mononuclear cell vaccinated vs candin placebo in infants (4-6m) (BCG-primed) after exposure to Modified Vaccinia Ankara (MVA) virus vaccine vector , time point 28D
• 27764254Genes negatively correlated with T cell responses (long term) in peripheral blood mononuclear cell in seniors (50-75) after exposure to Zostavax , time point 1D. Comment: (B) Network of genes informative of long-term responses.
• 15879137Genes up-regulated in polymorphonuclear leukocytes (24h): control versus infection by A. phagocytophilum.
• 29868000: Genes positively correlated with high anti-HBs concentration at week 30 in blood in young/old adults (20-40)/(60-84) (primary vaccination) after exposure to Twinrix , time point 1D. Comment: Correlation between pre-immunization expression levels of single genes (log2-transformed) and anti-HBs concentrations (log10-transformed) at week 30 post-primary vaccination
• 21357945: Genes positively correlated with titer response index in peripheral blood mononuclear cell in Caucasian male adults (18-40) (high responders) after exposure to Fluarix/Fluvirin , time point 1D and 3DY. Comment: Signature predictive of titer response index (TRI). Day 1 and day 3 values averaged.
• 15879137: Genes down-regulated in polymorphonuclear leukocytes (9h): control versus infection by A. phagocytophilum.
• 28193898: Genes up-regulated in peripheral blood mononuclear cell immunized with ARR vs immunized by RRR in unknown (primary immunization with recombinant adenovirus 35 (Ad35)) after exposure to P. falciparum RTS,S/AS01 , time point 1D
• 15789058: Genes down-regulated in comparison of neutrophils versusl monocytes.
• 17651872: Genes down-regulated in peripheral blood mononuclear cell post-vaccination vs pre-vaccination in adults (18-40) after exposure to YF-Vax or APSV Wetvax (identical responses) , time point anyD. Comment: Significantly Modulated Genes Common to Vaccinia and Yellow Fever Vaccination
• 26726811: Genes down-regulated in T cell 7d vs 1d in adults (18-64) after exposure to Pandemrix (A/California/7/09 (H1N1)) , time point 7D. Comment: - roughly 60/40 female:male ratio, over 70% were Causasian
• 15789058: Genes down-regulated in comparison of naive B cells versus day 0 monocytes.
• 29535712: Genes down-regulated in blood 3d and 7d vs 0hr in adults (18-45) (high IgM responders) after exposure to CN54gp140 adjuvanted with GLA-AF , time point 3D, 7D combined (identical signatures) , administered i.m.
• 17651872: Genes down-regulated in peripheral blood mononuclear cell in adults (18-40) after exposure to YF-Vax , time point anyD
• 24495909: Genes positively correlated with H3N2 VN titer in blood in children (0.5-14y) after exposure to FluMist , time point 7D. Comment: ~80% of cohort were white, ~50/50 Female:male
• 17349694: Genes up-regulated in peripheral blood mononuclear cell (18 to 336)h vs 0h in adults (22-54) after exposure to F. tularensis vaccine LVS , time point 18 to 336H. Comment: Pattern 3, sustained-up. These approx 9 of 42 genes in pattern linked to immune function.
• 28099485: Genes up-regulated in T cell 1d vs 0d in adults (18-49) after exposure to inactivated monovalent influenza A/Indonesia/05/2005 H5N1 split-virus vaccine , time point 1D , administered i.m.
• 21743478: Genes up-regulated in peripheral blood mononuclear cell 7d vs 0d in adults (18-50) after exposure to FluMist , time point 7D. Comment: Molecular signature induced by LAIV vaccination. (a) Interferon (IFN)-related genes differentially expressed after LAIV vaccination
• 21743478: Genes up-regulated in peripheral blood mononuclear cell 3d vs 0d in adults (18-50) after exposure to FluMist , time point 3D. Comment: Molecular signature induced by LAIV vaccination. (a) Interferon (IFN)-related genes differentially expressed after LAIV vaccination
• 26726811: Genes up-regulated in T cell 1d vs 0d in adults (18-64) after exposure to Pandemrix (A/California/7/09 (H1N1)) , time point 1D. Comment: - roughly 60/40 female:male ratio, over 70% were Causasian
• 27870591: Genes up-regulated in blood vaccinated vs control in adults (23-48) after exposure to Live attenuated vaccine TC-83 , time point 7D
• 26755593: Genes up-regulated in peripheral blood mononuclear cell 1d postboost vs 0d pre-imm in children (14-27m) (MF59-adjuvanted) after exposure to Fluad , time point 1D. Comment: ATIV
• 23420886: Genes down-regulated in thymic T reg: CD24 high [GeneID=100133941] versus CD24 int [GeneID=100133941].
• 20643338: Genes up-regulated in ex vivo follicular dendritic cells from peripheral lymph nodes: naïve versus immunized mice.
• 23420886: Genes down-regulated in T reg: peripheral lymph nodes versus thymic CD24 int [GeneID=100133941].
• 23844129: Genes up-regulated in peripheral blood mononuclear cell low responders vs high responders in adults (18-55) after exposure to Modified Vaccinia Ankara (MVA) virus vaccine vector , time point 2D. Comment: Enriched for GO terms associated with regulation of T-cell activation and co-stimulation signal (DAVID, fdr<0.05), from 176 DE genes
• 23878721: Genes positively correlated with antibody response in blood in adults (18-40) after exposure to Sanofi Pasteur, SA, Inactivated influenza vaccine , time point 1D
• 26755593: Genes up-regulated in peripheral blood mononuclear cell 1d postboost vs 0d pre-imm in children (14-27m) (MF59-adjuvanted) after exposure to Fluad , time point 1D. Comment: (C) Genes in BTM M40; (D) Genes in BTM M53
• 17105821: Genes up-regulated in comparison of peripheral blood mononuclear cells (PBMC) from patients with acute influenza infection versus PBMC from patients with acute E. coli infection.
• 21743478: Genes down-regulated in peripheral blood mononuclear cell 7d vs 0d in adults (18-50) after exposure to FluMist , time point 7D. Comment: Molecular signature induced by LAIV vaccination. (a) Interferon (IFN)-related genes differentially expressed after LAIV vaccination
• 15789058: Genes down-regulated in comparison of naive CD8 T cells versus day 0 monocytes.
• 17595242: Genes up-regulated in comparison of peripheral blood mononuclear cells (PBMC) from healthy donors versus PBMCs from patients with type 1 diabetes at 4 month after the diagnosis.
• 24336226: Genes negatively correlated with antibody response in peripheral blood mononuclear cell in adults (18-45) (anti-DT antibody-correlation profile) after exposure to Menactra , time point 3D
• 24808365: Genes up-regulated in dendritic cells infected by Leishmania major: 2h versus 24h.
• 21743478: Genes down-regulated in comparison of B cells versus myeloid dendritic cells (mDC).
• 21743478: Genes down-regulated in comparison of plasmacytoid dendritic cells (DC) versus myeloid DCs.
• 26755593: Genes positively correlated with HAI response in peripheral blood mononuclear cell in children (14-27m) (MF59-adjuvanted and non-adjuvanted) after exposure to Fluad/Imuvac , time point 1D. Comment: Genes in BTM M75
• 22617845: Genes down-regulated in peripheral blood mononuclear cell 24h vs 0h in adults (18-45) (non-responders (previously immunized)) after exposure to Live attenuated vaccine TC-83 , time point 24H. Comment: initial exposure 2-10 months before PBMCs drawn. significant genes chosen for membership in canonical pathways
• 27870591: Genes up-regulated in blood vaccinated vs control in adults (23-48) after exposure to Live attenuated vaccine TC-83 , time point 2D
• 29535712: Genes down-regulated in blood 1d vs 0hr in adults (18-45) after exposure to CN54gp140 adjuvanted with GLA-AF , time point 1D , administered i.m.
• 22722857: Genes down-regulated in CD4 [GeneID=920] over-expressing: FOXP3 [GeneID=50943] and PPARg1 form of PPARG [GeneID=5468] versus FOXP3 [GeneID=50943].
• 19029902: Genes up-regulated in peripheral blood mononuclear cell 3d vs 0d in adults (18-45) after exposure to YF-17D vaccine , time point 3D
• 19029902: Genes up-regulated in peripheral blood mononuclear cell 7d vs 0d in adults (18-45) after exposure to YF-17D vaccine , time point 7D
• 22722857: Genes down-regulated in CD4 [GeneID=920] T cells treated with pioglitazone [PubChem=4829] and over-expressing: FOXP3 [GeneID=50943] and PPARg1 isoform of PPARG [GeneID=5468] versus FOXP3 [GeneID=50943] and PPARg2 form of PPARG [GeneID=5468].
• 21743478: Genes down-regulated in comparison of plasmacytoid dendritic cells (DC) from influenza vaccinee at day 7 post-vaccination versus myeloid DCs at day 7 post-vaccination.
• 24912498: Genes up-regulated in blood vaccinated vs candin placebo in infants (4-6m) (BCG-primed) after exposure to Modified Vaccinia Ankara (MVA) virus vaccine vector , time point 1D
• 21636294: Genes down-regulated in BCL6 [GeneID=604] high follicular helper T cells (Tfh) versus all Tfh.
• 15789058: Genes down-regulated in comparison of naive CD4 [GeneID=920] T cells versus day 0 monocytes.
• 17204652: Genes down-regulated in cells from Flt3L Melanom injected mice: splenic DEC205+ dendritic cells versus CD8 T cells.

3.1.2.6 Tumaco biopsies H

t_biopsy_gsva_h <- simple_gsva(
    t_biopsies,
    signatures = broad_h,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "h")
t_biopsy_gsva_h_sig <- get_sig_gsva_categories(
    t_biopsy_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_biopsy_gsva_h.xlsx")

t_biopsy_gsva_h <- simple_gsva(
    t_biopsies, signature_category = "h")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_biopsy_gsva_h_sig <- get_sig_gsva_categories(
    t_biopsy_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_biopsy_gsva_h.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 9 rows.

## The factor tumaco_failure has 5 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.7 Tumaco Eosinophils C7

t_eosinophil_gsva_c7 <- simple_gsva(
    t_eosinophils,
    signatures = broad_c7,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c7")
t_eosinophil_gsva_c7_sig <- get_sig_gsva_categories(
    t_eosinophil_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_eosinophil_gsva_c7.xlsx")

t_eosinophil_gsva_c7 <- simple_gsva(
    t_eosinophils, signature_category = "c7")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_eosinophil_gsva_c7_sig <- get_sig_gsva_categories(
    t_eosinophil_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_eosinophil_gsva_c7.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 17 rows.

## The factor tumaco_failure has 9 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.8 Tumaco Eosinophils C2

t_eosinophil_gsva_c2 <- simple_gsva(
    t_eosinophils,
    signatures = broad_c2,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c2")
t_eosinophil_gsva_c2_sig <- get_sig_gsva_categories(
    t_eosinophil_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_eosinophil_gsva_c2.xlsx")

t_eosinophil_gsva_c2 <- simple_gsva(
    t_eosinophils, signature_category = "c2")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_eosinophil_gsva_c2_sig <- get_sig_gsva_categories(
    t_eosinophil_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_eosinophil_gsva_c2.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 17 rows.

## The factor tumaco_failure has 9 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.9 Tumaco Eosinophils H

t_eosinophil_gsva_h <- simple_gsva(
    t_eosinophils,
    signatures = broad_h,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "h")
t_eosinophil_gsva_h_sig <- get_sig_gsva_categories(
    t_eosinophil_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_eosinophil_gsva_h.xlsx")

t_eosinophil_gsva_h <- simple_gsva(
    t_eosinophils, signature_category = "h")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_eosinophil_gsva_h_sig <- get_sig_gsva_categories(
    t_eosinophil_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_eosinophil_gsva_h.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 17 rows.

## The factor tumaco_failure has 9 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.10 Tumaco Monocytes C7

t_monocyte_gsva_c7 <- simple_gsva(
    t_monocytes,
    signatures = broad_c7,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c7")
t_monocyte_gsva_c7_sig <- get_sig_gsva_categories(
    t_monocyte_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_monocyte_gsva_c7.xlsx")

t_monocyte_gsva_c7 <- simple_gsva(
    t_monocytes, signature_category = "c7")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_monocyte_gsva_c7_sig <- get_sig_gsva_categories(
    t_monocyte_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_monocyte_gsva_c7.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 21 rows.

## The factor tumaco_failure has 21 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.11 Tumaco Monocytes C2

t_monocyte_gsva_c2 <- simple_gsva(
    t_monocytes,
    signatures = broad_c2,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c2")
t_monocyte_gsva_c2_sig <- get_sig_gsva_categories(
    t_monocyte_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_monocyte_gsva_c2.xlsx")

t_monocyte_gsva_c2 <- simple_gsva(
    t_monocytes, signature_category = "c2")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_monocyte_gsva_c2_sig <- get_sig_gsva_categories(
    t_monocyte_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_monocyte_gsva_c2.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 21 rows.

## The factor tumaco_failure has 21 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.12 Tumaco Monocytes H

t_monocyte_gsva_h <- simple_gsva(
    t_monocytes,
    signatures = broad_h,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "h")
t_monocyte_gsva_h_sig <- get_sig_gsva_categories(
    t_monocyte_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_monocyte_gsva_h.xlsx")

t_monocyte_gsva_h <- simple_gsva(
    t_monocytes, signature_category = "h")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_monocyte_gsva_h_sig <- get_sig_gsva_categories(
    t_monocyte_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_monocyte_gsva_h.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 21 rows.

## The factor tumaco_failure has 21 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.13 Tumaco Neutrophils c7

t_neutrophil_gsva_c7 <- simple_gsva(
    t_neutrophils,
    signatures = broad_c7,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c7")
t_neutrophil_gsva_c7_sig <- get_sig_gsva_categories(
    t_neutrophil_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_neutrophil_gsva_c7.xlsx")

t_neutrophil_gsva_c7 <- simple_gsva(
    t_neutrophils, signature_category = "c7")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_neutrophil_gsva_c7_sig <- get_sig_gsva_categories(
    t_neutrophil_gsva_c7,
    excel = "analyses/4_tumaco/GSVA/t_neutrophil_gsva_c7.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 20 rows.

## The factor tumaco_failure has 21 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.14 Tumaco Neutrophils c2

t_neutrophil_gsva_c2 <- simple_gsva(
    t_neutrophils,
    signatures = broad_c2,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "c2")
t_neutrophil_gsva_c2_sig <- get_sig_gsva_categories(
    t_neutrophil_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_neutrophil_gsva_c2.xlsx")

t_neutrophil_gsva_c2 <- simple_gsva(
    t_neutrophils, signature_category = "c2")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_neutrophil_gsva_c2_sig <- get_sig_gsva_categories(
    t_neutrophil_gsva_c2,
    excel = "analyses/4_tumaco/GSVA/t_neutrophil_gsva_c2.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 20 rows.

## The factor tumaco_failure has 21 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

3.1.2.15 Tumaco Neutrophils H

t_neutrophil_gsva_h <- simple_gsva(
    t_neutrophils,
    signatures = broad_h,
    msig_xml = "reference/msigdb/msigdb_v7.5.1.xml",
    signature_category = "h")
t_neutrophil_gsva_h_sig <- get_sig_gsva_categories(
    t_neutrophil_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_neutrophil_gsva_h.xlsx")

t_neutrophil_gsva_h <- simple_gsva(
    t_neutrophils, signature_category = "h")

## Converting the rownames() of the expressionset to ENTREZID.

## 574 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 19952 entries.

## After conversion, the expressionset has 19337 entries.

## Adding descriptions and IDs to the gene set annotations.

t_neutrophil_gsva_h_sig <- get_sig_gsva_categories(
    t_neutrophil_gsva_h,
    excel = "analyses/4_tumaco/GSVA/t_neutrophil_gsva_h.xlsx")

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Assuming this data is similar to a micro array and not performign voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: tumacofailure_vs_tumacocure.  Adjust = BH

## Limma step 6/6: 1/2: Creating table: tumacocure.  Adjust = BH

## Limma step 6/6: 2/2: Creating table: tumacofailure.  Adjust = BH

## The factor tumaco_cure has 20 rows.

## The factor tumaco_failure has 21 rows.

## Testing each factor against the others.

## Scoring tumaco_cure against everything else.

## Scoring tumaco_failure against everything else.

Bibliography

“GSVAdata.” n.d. Bioconductor. http://bioconductor.org/packages/GSVAdata/. Accessed September 16, 2024.

Hänzelmann, Sonja, Robert Castelo, and Justin Guinney. 2013. “GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data.” BMC Bioinformatics 14 (1): 7. https://doi.org/10.1186/1471-2105-14-7.

Liberzon, Arthur, Aravind Subramanian, Reid Pinchback, Helga Thorvaldsdóttir, Pablo Tamayo, and Jill P. Mesirov. 2011. “Molecular Signatures Database (MSigDB) 3.0.” Bioinformatics 27 (12): 1739–40. https://doi.org/10.1093/bioinformatics/btr260.

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” bioRxiv. https://doi.org/10.1101/002832.

McCarthy, Davis J., Yunshun Chen, and Gordon K. Smyth. 2012. “Differential Expression Analysis of Multifactor RNA-Seq Experiments with Respect to Biological Variation.” Nucleic Acids Research 40 (10): 4288–97. https://doi.org/10.1093/nar/gks042.