L. panamensis 202207: 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

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 5633 low-count genes (14290 remaining).

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

## Warning in deseq_lrt(tc_clinical_filt, transform = "vst", interaction = FALSE, :
## The clinic 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

## -- 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 6272 genes.

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

## Joining, by = "merge"
## Joining, by = "merge"

summary(lrt_visit[["favorite_genes"]])

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

t_clinical_filt <- normalize_expt(t_clinical, filter=TRUE)

## Removing 5774 low-count genes (14149 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 5439 genes.

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

## Joining, by = "merge"
## Joining, by = "merge"

summary(lrt_visit[["favorite_genes"]])

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

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

## Deleting the file excel/lrt_clinical_visit-v202207.xlsx before writing the tables.

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, by = "merge"Joining, by = "merge"

lrt_monocyte_visit$cluster_data$plot

lrt_monocyte_visit_v2 <- deseq_lrt(t_visitcf_monocyte, transform = "vst",
                                   interaction = FALSE, minc = 1,
                                   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 > 1
## Joining, by = "merge"Joining, by = "merge"

2.2 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 552 genes.

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

## Joining, by = "merge"
## Joining, by = "merge"

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::glue("excel/lrt_clinical_celltype-v{ver}.xlsx"))

## Deleting the file excel/lrt_clinical_celltype-v202207.xlsx before writing the tables.

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

3 GSVA

3.0.1 Clinical samples

hs_celltype_gsva_c2 <- simple_gsva(tc_valid)

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

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

## Before conversion, the expressionset has 19923 entries.

## After conversion, the expressionset has 19549 entries.

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

## Starting limma pairwise comparison.

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

## Choosing the non-intercept containing 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)

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

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

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

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

hs_celltype_gsva_c2_sig$subset_plot

hs_celltype_gsva_c2_sig$score_plot

broad_c7 <- load_gmt_signatures(signatures="reference/msigdb/c7.all.v7.2.entrez.gmt",
                                signature_category="c7")
hs_celltype_gsva_c7 <- simple_gsva(tc_valid,
                                   signatures="reference/msigdb/c7.all.v7.2.entrez.gmt",
                                   signature_category="c7",
                                   msig_xml="reference/msigdb_v7.2.xml",
                                   cores=10)

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

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

## Before conversion, the expressionset has 19923 entries.

## After conversion, the expressionset has 19549 entries.

## Adding annotations from reference/msigdb_v7.2.xml.

hs_celltype_gsva_c7_sig <- get_sig_gsva_categories(
    hs_celltype_gsva_c7,
    excel="excel/individual_celltypes_gsva_c7.xlsx")

## Starting limma pairwise comparison.

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

## Choosing the non-intercept containing 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)

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

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

## Testing each factor against the others.

## Scoring cure against everything else.

## Scoring failure against everything else.

3.0.2 Print some plots of the GSVA outputs

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

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

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

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

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

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

tmp <- loadme(filename=savefile)