Explore the RTX data

J. Taroni 2018

Functions and directory set up

# magrittr pipe
`%>%` <- dplyr::`%>%`
source(file.path("util", "plier_util.R"))
source(file.path("util", "test_LV_differences.R"))

# plot and result directory setup for this notebook
plot.dir <- file.path("plots", "24")
dir.create(plot.dir, recursive = TRUE, showWarnings = FALSE)
results.dir <- file.path("results", "24")
dir.create(results.dir, recursive = TRUE, showWarnings = FALSE)

Read in data

Specifically, we’ll look at RNA-seq data that’s been processed various ways (e.g., taking into account experimental design or not, etc.).

Covariate information

covariate.df <- readr::read_tsv(file.path("data", "rtx", 
                                          "RTX_full_covariates.tsv"))

Parsed with column specification:
cols(
  .default = col_character(),
  barcode = col_integer(),
  AGE = col_integer(),
  bcells = col_double(),
  HGB = col_double(),
  `Platelet Count` = col_double(),
  WBC = col_double(),
  Lymphs = col_double(),
  Neutrophils = col_double(),
  Eosinophils = col_double(),
  Tscore = col_integer()
)
See spec(...) for full column specifications.

Counts

# tximport gene-level data -- includes counts, TPM
gene.summary <- readRDS(file.path("data", "rtx", "tximport.RDS"))
# extract count matrix -- we used `countsFromAbundance = "no"` to generate this
count.matrix <- gene.summary$counts
# log2 transform plus some constant, here 1, to account for zeros
log.count <- log2(count.matrix + 1)

VST - blind to experimental design

vst.blind <- readr::read_tsv(file.path("data", "rtx", "VST_blind.pcl"))

Parsed with column specification:
cols(
  .default = col_double(),
  Gene = col_character()
)
See spec(...) for full column specifications.

VST - batch, timepoint, responder status

Note that the experimental design we gave the DESeq2::vst function may ultimately not be the one we’d like to use.

vst.design <- readr::read_tsv(file.path("data", "rtx", "VST_design.pcl"))

Parsed with column specification:
cols(
  .default = col_double(),
  Gene = col_character()
)
See spec(...) for full column specifications.

PCA on full expression matrix

PCAPlotWrapper <- function(exprs, plot.title) {
  # This is only intended to be used in this global environment where we have
  # covariate.df and the magrittr pipe already
  # takes the rituximab expression data (exprs) where rows are genes and columns 
  # are samples and makes a plot of PC1-2 with the title supplied as plot.title
  # argument
  
  # PCA
  pc.results <- prcomp(t(exprs))
  cum.var.exp <- cumsum(pc.results$sdev^2 / sum(pc.results$sdev^2))
  # PC1-2 in form suitable for ggplot2
  pc.df <- as.data.frame(cbind(rownames(pc.results$x), pc.results$x[, 1:2]))
  colnames(pc.df)[1] <- "Sample"
  
  # check that the barcodes are in agreement
  barcode.check <- all(covariate.df$barcode == substr(pc.df[[1]], start = 1, 
                                                      stop = 6))
  if (!barcode.check) {
    stop("Something went wrong, the sample barcodes are not in the same order!")
  }
  
  # now that we've ensured that the barcodes are in the same order, we'll 
  # add in the batch information
  pc.df$Batch <- covariate.df$procbatch
  
  # plot!
  pc.df %>% 
    dplyr::mutate(PC1 = as.numeric(as.character(PC1)),
                  PC2 = as.numeric(as.character(PC2))) %>%
    ggplot2::ggplot(ggplot2::aes(x = PC1, y = PC2, color = Batch, 
                                 shape = Batch)) +
    ggplot2::geom_point(size = 3, alpha = 0.75) +
    ggplot2::labs(x = paste("PC1 (cum. var. exp. =", round(cum.var.exp[1], 3), 
                            ")"), 
                  y =  paste("PC2 (cum. var. exp. =", round(cum.var.exp[2], 3), 
                            ")"),
                  title = plot.title) +
    ggplot2::theme_bw() +
    ggplot2::theme(text = ggplot2::element_text(size = 15)) +
    ggplot2::scale_color_manual(values = c("#000000", "#969696"))
  
}

# the count matrix itself does not have column names because of the way it was
# processed
colnames(log.count) <- colnames(vst.blind)[2:ncol(vst.blind)]
PCAPlotWrapper(exprs = log.count, plot.title = "log2(counts + 1)")

PCAPlotWrapper(exprs = as.matrix(vst.blind[, 2:ncol(vst.blind)]),
               plot.title = "VST (blind)")

PCAPlotWrapper(exprs = as.matrix(vst.design[, 2:ncol(vst.design)]),
               plot.title = "VST")

Gene identifier conversion

All of the gene expression data from the RTX data set uses Ensembl gene IDs. PLIER, on the other hand, uses gene symbols. So we’ll need to do conversion. Let’s continue with variance stabilizing transformed data blinded to experimental design.

mart <- biomaRt::useDataset("hsapiens_gene_ensembl", 
                            biomaRt::useMart("ensembl"))
gene.df <- biomaRt::getBM(filters = "ensembl_gene_id",
                          attributes = c("ensembl_gene_id", "hgnc_symbol"),
                          values = vst.blind$Gene, 
                          mart = mart)

Batch submitting query [=---------------------------------]   3% eta: 26s
Batch submitting query [=---------------------------------]   4% eta: 29s
Batch submitting query [==--------------------------------]   6% eta: 29s
Batch submitting query [==--------------------------------]   7% eta: 29s
Batch submitting query [===-------------------------------]   9% eta: 29s
Batch submitting query [===-------------------------------]  10% eta: 31s
Batch submitting query [====------------------------------]  11% eta: 36s
Batch submitting query [====------------------------------]  13% eta: 34s
Batch submitting query [=====-----------------------------]  14% eta: 34s
Batch submitting query [=====-----------------------------]  16% eta: 33s
Batch submitting query [======----------------------------]  17% eta: 32s
Batch submitting query [======----------------------------]  19% eta: 32s
Batch submitting query [=======---------------------------]  20% eta: 31s
Batch submitting query [=======---------------------------]  21% eta: 30s
Batch submitting query [========--------------------------]  23% eta: 29s
Batch submitting query [========--------------------------]  24% eta: 29s
Batch submitting query [=========-------------------------]  26% eta: 28s
Batch submitting query [=========-------------------------]  27% eta: 27s
Batch submitting query [==========------------------------]  29% eta: 27s
Batch submitting query [==========------------------------]  30% eta: 27s
Batch submitting query [===========-----------------------]  31% eta: 27s
Batch submitting query [===========-----------------------]  33% eta: 27s
Batch submitting query [============----------------------]  34% eta: 26s
Batch submitting query [============----------------------]  36% eta: 26s
Batch submitting query [=============---------------------]  37% eta: 25s
Batch submitting query [=============---------------------]  39% eta: 24s
Batch submitting query [==============--------------------]  40% eta: 24s
Batch submitting query [==============--------------------]  41% eta: 23s
Batch submitting query [===============-------------------]  43% eta: 22s
Batch submitting query [===============-------------------]  44% eta: 22s
Batch submitting query [================------------------]  46% eta: 21s
Batch submitting query [================------------------]  47% eta: 20s
Batch submitting query [=================-----------------]  49% eta: 22s
Batch submitting query [=================-----------------]  50% eta: 22s
Batch submitting query [=================-----------------]  51% eta: 21s
Batch submitting query [==================----------------]  53% eta: 20s
Batch submitting query [==================----------------]  54% eta: 20s
Batch submitting query [===================---------------]  56% eta: 19s
Batch submitting query [===================---------------]  57% eta: 18s
Batch submitting query [====================--------------]  59% eta: 18s
Batch submitting query [====================--------------]  60% eta: 17s
Batch submitting query [=====================-------------]  61% eta: 17s
Batch submitting query [=====================-------------]  63% eta: 17s
Batch submitting query [======================------------]  64% eta: 16s
Batch submitting query [======================------------]  66% eta: 16s
Batch submitting query [=======================-----------]  67% eta: 15s
Batch submitting query [=======================-----------]  69% eta: 14s
Batch submitting query [========================----------]  70% eta: 14s
Batch submitting query [========================----------]  71% eta: 13s
Batch submitting query [=========================---------]  73% eta: 13s
Batch submitting query [=========================---------]  74% eta: 13s
Batch submitting query [==========================--------]  76% eta: 12s
Batch submitting query [==========================--------]  77% eta: 11s
Batch submitting query [===========================-------]  79% eta: 10s
Batch submitting query [===========================-------]  80% eta: 10s
Batch submitting query [============================------]  81% eta:  9s
Batch submitting query [============================------]  83% eta:  8s
Batch submitting query [=============================-----]  84% eta:  8s
Batch submitting query [=============================-----]  86% eta:  7s
Batch submitting query [==============================----]  87% eta:  6s
Batch submitting query [==============================----]  89% eta:  6s
Batch submitting query [===============================---]  90% eta:  5s
Batch submitting query [===============================---]  91% eta:  4s
Batch submitting query [================================--]  93% eta:  4s
Batch submitting query [================================--]  94% eta:  3s
Batch submitting query [=================================-]  96% eta:  2s
Batch submitting query [=================================-]  97% eta:  1s
Batch submitting query [==================================]  99% eta:  1s
Batch submitting query [==================================] 100% eta:  0s
                                                                         

# filter to remove genes without a gene symbol
gene.df <- gene.df %>% dplyr::filter(complete.cases(.))

vst.blind.annot <- dplyr::inner_join(gene.df, vst.blind,
                                     by = c("ensembl_gene_id" = "Gene")) %>%
  dplyr::select(-ensembl_gene_id)
colnames(vst.blind.annot)[1] <- "Gene"
# aggregate duplicate gene symbols to gene mean
vst.agg <- PrepExpressionDF(vst.blind.annot)

Read in recount2 PLIER model

recount.plier <- readRDS(file.path("data", "recount2_PLIER_data", 
                                   "recount_PLIER_model.RDS"))

Now filter the VST data to just the genes included in the recount2 model and prep it for use with various PLIER helper functions.

vst.filt <- vst.agg %>% 
  dplyr::filter(Gene %in% rownames(recount.plier$Z)) %>%
  tibble::column_to_rownames(var = "Gene")

Setting row names on a tibble is deprecated.

# save to file
vst.file <- file.path("data", "rtx", "VST_blind_filtered.RDS")
saveRDS(vst.filt, vst.file)

# let's remove the other VST data.frames or matrices and extraneous biomart 
# objects
rm(vst.agg, vst.blind, vst.blind.annot, vst.design, mart, count.matrix)

PCAPlotWrapper(exprs = vst.filt,
               plot.title = "VST filtered to genes in model")

PLIER

We’ll give PLIER row-normalized (z-scored) TPM data, as this is closest to the form we had the recount2 data in (RPKM) and, given this use case, the slight difference is unlikely to effect the results much.

# get the gene-level TPM
rtx.tpm <- gene.summary$abundance
# we'll need to convert to HGNC gene symbol
rtx.tpm <- tibble::rownames_to_column(data.frame(rtx.tpm), var = "Gene")
rtx.tpm.annot <- dplyr::inner_join(gene.df, rtx.tpm,
                                    by = c("ensembl_gene_id" = "Gene")) %>%
  dplyr::select(-ensembl_gene_id)
colnames(rtx.tpm.annot) <- c("Gene", colnames(log.count))
# Aggregate duplicate gene identifiers to mean
rtx.agg <- PrepExpressionDF(rtx.tpm.annot)
# remove first row without identifier & write to file
rtx.agg <- rtx.agg %>%
  dplyr::filter(Gene != "")
readr::write_tsv(rtx.agg, file.path("data", "rtx", "tpm_aggregated.pcl"))
# get in expression matrix form
rtx.exprs <- as.matrix(rtx.agg %>% tibble::column_to_rownames(var = "Gene"))

Setting row names on a tibble is deprecated.

# remove intermediates
rm(rtx.tpm.annot, rtx.tpm, rtx.agg)

# remove genes that are all zero
remove.index <- which(apply(rtx.exprs, 1, function(x) all(x == 0)))
rtx.exprs <- rtx.exprs[-remove.index, ]

Dataset-specific model

Training a PLIER model on this dataset

rtx.plier <- PLIERNewData(exprs.mat = as.matrix(rtx.exprs))

Loading required package: PLIER
Loading required package: RColorBrewer
Loading required package: gplots

Attaching package: ‘gplots’

The following object is masked from ‘package:stats’:

    lowess

Loading required package: pheatmap
Loading required package: glmnet
Loading required package: Matrix
Loading required package: foreach
Loaded glmnet 2.0-13

Loading required package: knitr
Loading required package: rsvd
Loading required package: qvalue
Removing 1 pathways with too few genesComputing SVD
Done

[1] 58.36176
[1] "L2 is set to 58.361764073198"
[1] "L1 is set to 29.180882036599"

errorY (SVD based:best possible) = 0.1835
iter1 errorY= 0.9724, Bdiff= 5546, Bkappa=63.6
iter2 errorY= 0.4744, Bdiff= 0.6175, Bkappa=40.06
iter3 errorY= 0.3814, Bdiff= 0.09271, Bkappa=32.25
iter4 errorY= 0.3606, Bdiff= 0.03585, Bkappa=28.31
iter5 errorY= 0.3526, Bdiff= 0.01795, Bkappa=26.55
iter6 errorY= 0.3493, Bdiff= 0.009761, Bkappa=25.69
iter7 errorY= 0.347, Bdiff= 0.006178, Bkappa=25.2
iter8 errorY= 0.3447, Bdiff= 0.004487, Bkappa=25.13
iter9 errorY= 0.3425, Bdiff= 0.003575, Bkappa=24.3
iter10 errorY= 0.3405, Bdiff= 0.002988, Bkappa=24.48
iter11 errorY= 0.3388, Bdiff= 0.002519, Bkappa=24.44
iter12 errorY= 0.3374, Bdiff= 0.002119, Bkappa=22.51
iter13 errorY= 0.3363, Bdiff= 0.001784, Bkappa=22.17
iter14 errorY= 0.3354, Bdiff= 0.00151, Bkappa=24.51
iter15 errorY= 0.3347, Bdiff= 0.001287, Bkappa=23.86
iter16 errorY= 0.3341, Bdiff= 0.001118, Bkappa=23.51
iter17 errorY= 0.3335, Bdiff= 0.0009953, Bkappa=20.35
iter18 errorY= 0.3331, Bdiff= 0.0009098, Bkappa=19.44
iter19 errorY= 0.3327, Bdiff= 0.0008494, Bkappa=19.55
Updating L3, current fraction= 0, target=0.7
L3 is set to 0.01172 in 8 iterations
iter20 errorY= 0.3301, prior information ratio= 0.03, Bdiff= 0.0009455, Bkappa=19.41;pos. col. U=17
iter21 errorY= 0.3292, prior information ratio= 0.03, Bdiff= 0.0008496, Bkappa=19.37;pos. col. U=16
iter22 errorY= 0.3286, prior information ratio= 0.05, Bdiff= 0.0007274, Bkappa=19.29;pos. col. U=14
iter23 errorY= 0.3282, prior information ratio= 0.05, Bdiff= 0.0006296, Bkappa=20.19;pos. col. U=14
iter24 errorY= 0.3279, prior information ratio= 0.05, Bdiff= 0.0005488, Bkappa=20.44;pos. col. U=14
iter25 errorY= 0.3276, prior information ratio= 0.04, Bdiff= 0.0004819, Bkappa=20.63;pos. col. U=15
iter26 errorY= 0.3274, prior information ratio= 0.05, Bdiff= 0.0004247, Bkappa=20.81;pos. col. U=15
iter27 errorY= 0.3272, prior information ratio= 0.05, Bdiff= 0.0003731, Bkappa=20.99;pos. col. U=15
iter28 errorY= 0.327, prior information ratio= 0.04, Bdiff= 0.0003264, Bkappa=21.13;pos. col. U=16
iter29 errorY= 0.3268, prior information ratio= 0.04, Bdiff= 0.0002838, Bkappa=21.27;pos. col. U=16
iter30 errorY= 0.3267, prior information ratio= 0.04, Bdiff= 0.0002474, Bkappa=21.43;pos. col. U=16
iter31 errorY= 0.3265, prior information ratio= 0.05, Bdiff= 0.0002154, Bkappa=21.91;pos. col. U=16
iter32 errorY= 0.3264, prior information ratio= 0.04, Bdiff= 0.0001876, Bkappa=21.97;pos. col. U=17
iter33 errorY= 0.3263, prior information ratio= 0.04, Bdiff= 0.0001639, Bkappa=21.96;pos. col. U=17
iter34 errorY= 0.3262, prior information ratio= 0.04, Bdiff= 0.0001442, Bkappa=21.89;pos. col. U=17
iter35 errorY= 0.3261, prior information ratio= 0.04, Bdiff= 0.0001276, Bkappa=21.76;pos. col. U=17
iter36 errorY= 0.326, prior information ratio= 0.04, Bdiff= 0.0001132, Bkappa=21.57;pos. col. U=17
iter37 errorY= 0.326, prior information ratio= 0.04, Bdiff= 0.0001006, Bkappa=21.32;pos. col. U=17
iter38 errorY= 0.3259, prior information ratio= 0.04, Bdiff= 8.93e-05, Bkappa=21.02;pos. col. U=17
iter39 errorY= 0.3258, prior information ratio= 0.04, Bdiff= 7.94e-05, Bkappa=20.76;pos. col. U=17
Updating L3, current fraction= 0.7391, target=0.7
L3 not changed
iter40 errorY= 0.3258, prior information ratio= 0.04, Bdiff= 7.111e-05, Bkappa=20.54;pos. col. U=17
iter41 errorY= 0.3257, prior information ratio= 0.04, Bdiff= 6.371e-05, Bkappa=20.32;pos. col. U=17
iter42 errorY= 0.3256, prior information ratio= 0.04, Bdiff= 5.733e-05, Bkappa=20.1;pos. col. U=17
iter43 errorY= 0.3256, prior information ratio= 0.04, Bdiff= 5.18e-05, Bkappa=19.88;pos. col. U=17
iter44 errorY= 0.3255, prior information ratio= 0.04, Bdiff= 4.707e-05, Bkappa=19.66;pos. col. U=17
iter45 errorY= 0.3255, prior information ratio= 0.04, Bdiff= 4.307e-05, Bkappa=19.48;pos. col. U=17
iter46 errorY= 0.3254, prior information ratio= 0.04, Bdiff= 3.973e-05, Bkappa=19.31;pos. col. U=17
iter47 errorY= 0.3254, prior information ratio= 0.04, Bdiff= 3.69e-05, Bkappa=19.14;pos. col. U=17
iter48 errorY= 0.3253, prior information ratio= 0.04, Bdiff= 3.456e-05, Bkappa=18.98;pos. col. U=17
iter49 errorY= 0.3253, prior information ratio= 0.04, Bdiff= 3.255e-05, Bkappa=18.82;pos. col. U=17
iter50 errorY= 0.3252, prior information ratio= 0.04, Bdiff= 3.086e-05, Bkappa=18.5;pos. col. U=17
iter51 errorY= 0.3252, prior information ratio= 0.04, Bdiff= 2.947e-05, Bkappa=18.47;pos. col. U=17
iter52 errorY= 0.3251, prior information ratio= 0.04, Bdiff= 2.84e-05, Bkappa=18.43;pos. col. U=17
iter53 errorY= 0.3251, prior information ratio= 0.04, Bdiff= 2.762e-05, Bkappa=18.39;pos. col. U=17
iter54 errorY= 0.325, prior information ratio= 0.04, Bdiff= 2.704e-05, Bkappa=18.35;pos. col. U=17
iter55 errorY= 0.325, prior information ratio= 0.04, Bdiff= 2.667e-05, Bkappa=18.3;pos. col. U=17
iter56 errorY= 0.3249, prior information ratio= 0.04, Bdiff= 2.647e-05, Bkappa=18.26;pos. col. U=17
iter57 errorY= 0.3249, prior information ratio= 0.04, Bdiff= 2.642e-05, Bkappa=18.21;pos. col. U=17
iter58 errorY= 0.3249, prior information ratio= 0.04, Bdiff= 2.649e-05, Bkappa=18.17;pos. col. U=17
iter59 errorY= 0.3248, prior information ratio= 0.04, Bdiff= 2.671e-05, Bkappa=18.29;pos. col. U=17
Updating L3, current fraction= 0.7391, target=0.7
L3 not changed
iter60 errorY= 0.3248, prior information ratio= 0.04, Bdiff= 2.709e-05, Bkappa=18.26;pos. col. U=17
iter61 errorY= 0.3247, prior information ratio= 0.04, Bdiff= 2.763e-05, Bkappa=18.44;pos. col. U=17
iter62 errorY= 0.3247, prior information ratio= 0.04, Bdiff= 2.811e-05, Bkappa=18.43;pos. col. U=17
iter63 errorY= 0.3247, prior information ratio= 0.04, Bdiff= 2.852e-05, Bkappa=18.43;pos. col. U=17
iter64 errorY= 0.3246, prior information ratio= 0.04, Bdiff= 2.89e-05, Bkappa=18.42;pos. col. U=17
iter65 errorY= 0.3246, prior information ratio= 0.04, Bdiff= 2.914e-05, Bkappa=18.42;pos. col. U=16
iter66 errorY= 0.3245, prior information ratio= 0.04, Bdiff= 2.942e-05, Bkappa=18.42;pos. col. U=17
iter67 errorY= 0.3245, prior information ratio= 0.04, Bdiff= 2.96e-05, Bkappa=18.41;pos. col. U=17
iter68 errorY= 0.3245, prior information ratio= 0.04, Bdiff= 2.955e-05, Bkappa=18.41;pos. col. U=17
iter69 errorY= 0.3244, prior information ratio= 0.04, Bdiff= 2.934e-05, Bkappa=18.4;pos. col. U=17
iter70 errorY= 0.3244, prior information ratio= 0.04, Bdiff= 2.887e-05, Bkappa=18.39;pos. col. U=17
iter71 errorY= 0.3243, prior information ratio= 0.04, Bdiff= 2.825e-05, Bkappa=18.39;pos. col. U=17
iter72 errorY= 0.3243, prior information ratio= 0.04, Bdiff= 2.749e-05, Bkappa=18.38;pos. col. U=17
iter73 errorY= 0.3243, prior information ratio= 0.04, Bdiff= 2.634e-05, Bkappa=18.37;pos. col. U=17
iter74 errorY= 0.3242, prior information ratio= 0.04, Bdiff= 2.51e-05, Bkappa=18.36;pos. col. U=17
iter75 errorY= 0.3242, prior information ratio= 0.04, Bdiff= 2.378e-05, Bkappa=18.35;pos. col. U=17
iter76 errorY= 0.3242, prior information ratio= 0.04, Bdiff= 2.254e-05, Bkappa=18.37;pos. col. U=17
iter77 errorY= 0.3241, prior information ratio= 0.04, Bdiff= 2.139e-05, Bkappa=18.39;pos. col. U=17
iter78 errorY= 0.3241, prior information ratio= 0.04, Bdiff= 2.029e-05, Bkappa=18.4;pos. col. U=17
iter79 errorY= 0.3241, prior information ratio= 0.04, Bdiff= 1.922e-05, Bkappa=18.42;pos. col. U=17
Updating L3, current fraction= 0.7391, target=0.7
L3 not changed
iter80 errorY= 0.3241, prior information ratio= 0.04, Bdiff= 1.824e-05, Bkappa=18.44;pos. col. U=17
iter81 errorY= 0.324, prior information ratio= 0.04, Bdiff= 1.73e-05, Bkappa=18.46;pos. col. U=17
iter82 errorY= 0.324, prior information ratio= 0.04, Bdiff= 1.643e-05, Bkappa=18.48;pos. col. U=17
iter83 errorY= 0.324, prior information ratio= 0.04, Bdiff= 1.566e-05, Bkappa=18.5;pos. col. U=17
iter84 errorY= 0.324, prior information ratio= 0.04, Bdiff= 1.5e-05, Bkappa=18.53;pos. col. U=17
iter85 errorY= 0.324, prior information ratio= 0.04, Bdiff= 1.44e-05, Bkappa=18.3;pos. col. U=17
iter86 errorY= 0.324, prior information ratio= 0.04, Bdiff= 1.387e-05, Bkappa=18.32;pos. col. U=17
iter87 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.339e-05, Bkappa=18.34;pos. col. U=17
iter88 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.296e-05, Bkappa=18.37;pos. col. U=17
iter89 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.258e-05, Bkappa=18.39;pos. col. U=17
iter90 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.224e-05, Bkappa=18.41;pos. col. U=17
iter91 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.192e-05, Bkappa=18.44;pos. col. U=17
iter92 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.162e-05, Bkappa=18.46;pos. col. U=17
iter93 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.133e-05, Bkappa=18.51;pos. col. U=17
iter94 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.106e-05, Bkappa=18.58;pos. col. U=17
iter95 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.079e-05, Bkappa=18.65;pos. col. U=17
iter96 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.048e-05, Bkappa=18.72;pos. col. U=17
iter97 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 1.019e-05, Bkappa=18.78;pos. col. U=17
iter98 errorY= 0.3238, prior information ratio= 0.04, Bdiff= 9.922e-06, Bkappa=18.85;pos. col. U=17
converged at  iteration 98
There are 15  LVs with AUC>0.70

saveRDS(rtx.plier, file.path("data", "rtx", "RTX_PLIER_model.RDS"))

PLIER::plotU(rtx.plier, auc.cutoff = 0.75, fontsize_row = 7,
             fontsize_col = 10)

There’s no neutrophil signature with AUC > 0.75. That’s interesting because this is a signal that we would expect in this whole blood data.

pdf(file.path(plot.dir, "RTX_model_U_plot.pdf"))
PLIER::plotU(rtx.plier, auc.cutoff = 0.75, fontsize_row = 7,
             fontsize_col = 10)
dev.off()

null device 
          1 

PCAPlotWrapper(exprs = rtx.plier$B, 
               plot.title = "Dataset-specific PLIER B, All LVs")

recount2 model

recount.b <- GetNewDataB(exprs.mat = as.matrix(rtx.exprs),
                         plier.model = recount.plier)
saveRDS(recount.b, file.path("data", "rtx", "RTX_recount2_B.RDS"))

PCAPlotWrapper(exprs = recount.b, 
               plot.title = "recount2 PLIER B, All LVs")

Final plotting

We’ll plot the following three plots, as these are the expected into any kind of supervised machine learning model (predicting response):

• VST expression data (only genes that overlap with recount2 PLIER model)
• Dataset-specific PLIER B, all LVs
• recount2 PLIER B, all LVs

exprs.plot <- PCAPlotWrapper(exprs = vst.filt,
                             plot.title = "VST filtered to genes in model")
rtx.plot <- PCAPlotWrapper(exprs = rtx.plier$B, 
                           plot.title = "Dataset-specific PLIER B, All LVs")
recount.plot <- PCAPlotWrapper(exprs = recount.b, 
                               plot.title = "recount2 PLIER B, All LVs")

pdf(file.path(plot.dir, "RTX_expected_ML_input_PCA.pdf"), height = 14, 
    width = 7)
cowplot::plot_grid(exprs.plot, rtx.plot, recount.plot, align = "h", ncol = 1)
dev.off()

null device 
          1