PLIER Utils proof-of-concept

J. Taroni 2018

In this notebook, I aim to demonstrate the validity of implementation for a subset of PLIER custom functions. Specifically, checking the operations for:

• Training a PLIER model on a new dataset PLIERNewData()
• Applying a previously computed PLIER to a new dataset to get the LV x sample matrix (B) GetNewDataB()
• Reconstruction of input gene expression data with a PLIER model GetReconstructedExprs() and the evaluation function GetReconstructionCorrelation()

Here, I use the NARES dataset for convenience due to its relatively small size (n = 77).

Load libraries and custom functions

library(AnnotationDbi)

Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport,
    clusterMap, parApply, parCapply, parLapply, parLapplyLB, parRapply, parSapply,
    parSapplyLB

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, cbind, colMeans, colnames, colSums,
    do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect,
    is.unsorted, lapply, lengths, Map, mapply, match, mget, order, paste, pmax,
    pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rowMeans, rownames,
    rowSums, sapply, setdiff, sort, table, tapply, union, unique, unsplit, which,
    which.max, which.min

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite
    Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: IRanges
Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

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

    expand.grid

library(PLIER)

Loading required package: RColorBrewer
Loading required package: gplots

Attaching package: ‘gplots’

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

    space

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

    space

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

    lowess

Loading required package: pheatmap
Loading required package: glmnet
Loading required package: Matrix

Attaching package: ‘Matrix’

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

    expand

Loading required package: foreach
Loaded glmnet 2.0-13

Loading required package: knitr
Loading required package: rsvd
Loading required package: qvalue

source(file.path("util", "plier_util.R"))

# plots directory specifically for this notebook
dir.create(file.path("plots", "01"), recursive = TRUE,
           showWarnings = FALSE)

NARES expression data

# Read in PCL
nares.data <- readr::read_tsv(file.path("data", "expression_data", 
                                         "NARES_SCANfast_ComBat.pcl"),
                              progress = FALSE)

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

Building a PLIER model requires HGNC symbol annotation, as this is what is included in the prior information (e.g., pathways, genesets) that is included in the package.

symbol.obj <- org.Hs.eg.db::org.Hs.egSYMBOL

mapped.genes <- AnnotationDbi::mappedkeys(symbol.obj)
symbol.list <- as.list(symbol.obj[mapped.genes])
symbol.df <- as.data.frame(cbind(names(symbol.list), unlist(symbol.list)))
colnames(symbol.df) <- c("EntrezID", "GeneSymbol")
# get gene column name to match to facilitate use with dplyr
colnames(nares.data)[1] <- "EntrezID"
# matching types
symbol.df$EntrezID <- as.integer(as.character(symbol.df$EntrezID))
# inner join
annot.nares.data <- dplyr::inner_join(symbol.df, nares.data, by = "EntrezID")
# only leave the matrix with gene symbol as rownames
exprs.mat <- dplyr::select(annot.nares.data, -EntrezID)
rownames(exprs.mat) <- exprs.mat$GeneSymbol
exprs.mat <- as.matrix(dplyr::select(exprs.mat, -GeneSymbol))

Train PLIER model

PLIERNewData is a wrapper function for applying PLIER to a new dataset. Expression data is row-normalized for use with PLIER. We use the following genesets that come with PLIER: bloodCellMarkersIRISDMAP, svmMarkers, and canonicalPathways. See also the PLIER vignette.

nares.plier <- PLIERNewData(exprs.mat)

Removing 1 pathways with too few genesComputing SVD
Done

[1] 47.1281
[1] "L2 is set to 47.1281023042499"
[1] "L1 is set to 23.5640511521249"

errorY (SVD based:best possible) = 0.116
iter1 errorY= 0.9617, Bdiff= 3358, Bkappa=69.2
iter2 errorY= 0.4158, Bdiff= 0.5971, Bkappa=41.7
iter3 errorY= 0.3548, Bdiff= 0.06627, Bkappa=26.01
iter4 errorY= 0.3462, Bdiff= 0.02215, Bkappa=18.06
iter5 errorY= 0.3439, Bdiff= 0.0142, Bkappa=17.25
iter6 errorY= 0.3424, Bdiff= 0.01032, Bkappa=17.89
iter7 errorY= 0.3407, Bdiff= 0.008008, Bkappa=19.53
iter8 errorY= 0.3388, Bdiff= 0.006558, Bkappa=21.84
iter9 errorY= 0.3368, Bdiff= 0.005559, Bkappa=24.73
iter10 errorY= 0.3349, Bdiff= 0.004753, Bkappa=27.64
iter11 errorY= 0.3333, Bdiff= 0.004037, Bkappa=29.76
iter12 errorY= 0.3319, Bdiff= 0.003411, Bkappa=32.67
iter13 errorY= 0.3309, Bdiff= 0.002883, Bkappa=34.5
iter14 errorY= 0.3301, Bdiff= 0.002455, Bkappa=38.06
iter15 errorY= 0.3295, Bdiff= 0.002122, Bkappa=39.39
iter16 errorY= 0.329, Bdiff= 0.001864, Bkappa=41.45
iter17 errorY= 0.3286, Bdiff= 0.001669, Bkappa=36.21
iter18 errorY= 0.3283, Bdiff= 0.001511, Bkappa=38.81
iter19 errorY= 0.3279, Bdiff= 0.001381, Bkappa=45.1
Updating L3, current fraction= 0, target=0.7
L3 is set to 0.004884 in 10 iterations
iter20 errorY= 0.3261, prior information ratio= 0.03, Bdiff= 0.001376, Bkappa=45.6;pos. col. U=24
iter21 errorY= 0.3251, prior information ratio= 0.04, Bdiff= 0.001292, Bkappa=47.26;pos. col. U=24
iter22 errorY= 0.3245, prior information ratio= 0.04, Bdiff= 0.001206, Bkappa=49.19;pos. col. U=24
iter23 errorY= 0.3239, prior information ratio= 0.04, Bdiff= 0.001109, Bkappa=51.49;pos. col. U=23
iter24 errorY= 0.3234, prior information ratio= 0.05, Bdiff= 0.001001, Bkappa=54.24;pos. col. U=23
iter25 errorY= 0.3229, prior information ratio= 0.05, Bdiff= 0.0008939, Bkappa=57.87;pos. col. U=23
iter26 errorY= 0.3225, prior information ratio= 0.05, Bdiff= 0.0007924, Bkappa=62.17;pos. col. U=23
iter27 errorY= 0.3222, prior information ratio= 0.05, Bdiff= 0.0007015, Bkappa=65.38;pos. col. U=22
iter28 errorY= 0.3219, prior information ratio= 0.05, Bdiff= 0.0006201, Bkappa=69.5;pos. col. U=23
iter29 errorY= 0.3216, prior information ratio= 0.05, Bdiff= 0.0005482, Bkappa=73.66;pos. col. U=23
iter30 errorY= 0.3214, prior information ratio= 0.05, Bdiff= 0.0004862, Bkappa=77.86;pos. col. U=24
iter31 errorY= 0.3212, prior information ratio= 0.04, Bdiff= 0.0004342, Bkappa=82.14;pos. col. U=24
iter32 errorY= 0.321, prior information ratio= 0.04, Bdiff= 0.0003932, Bkappa=86.77;pos. col. U=24
iter33 errorY= 0.3208, prior information ratio= 0.05, Bdiff= 0.0003647, Bkappa=91.5;pos. col. U=25
iter34 errorY= 0.3207, prior information ratio= 0.04, Bdiff= 0.0003424, Bkappa=96.35;pos. col. U=25
iter35 errorY= 0.3205, prior information ratio= 0.04, Bdiff= 0.0003255, Bkappa=102.2;pos. col. U=25
iter36 errorY= 0.3203, prior information ratio= 0.04, Bdiff= 0.000313, Bkappa=108.4;pos. col. U=25
iter37 errorY= 0.3201, prior information ratio= 0.04, Bdiff= 0.0003047, Bkappa=117.6;pos. col. U=25
iter38 errorY= 0.3199, prior information ratio= 0.04, Bdiff= 0.0002974, Bkappa=125.3;pos. col. U=25
iter39 errorY= 0.3198, prior information ratio= 0.04, Bdiff= 0.0002906, Bkappa=133.7;pos. col. U=25
Updating L3, current fraction= 0.7353, target=0.7
L3 is set to 0.00571 in 12 iterations
iter40 errorY= 0.3197, prior information ratio= 0.05, Bdiff= 0.0002853, Bkappa=142.3;pos. col. U=23
iter41 errorY= 0.3196, prior information ratio= 0.05, Bdiff= 0.0002782, Bkappa=150.6;pos. col. U=21
iter42 errorY= 0.3195, prior information ratio= 0.05, Bdiff= 0.0002711, Bkappa=159.2;pos. col. U=21
iter43 errorY= 0.3193, prior information ratio= 0.05, Bdiff= 0.0002642, Bkappa=167.7;pos. col. U=21
iter44 errorY= 0.3192, prior information ratio= 0.05, Bdiff= 0.0002559, Bkappa=174.7;pos. col. U=21
iter45 errorY= 0.3191, prior information ratio= 0.05, Bdiff= 0.000246, Bkappa=182.7;pos. col. U=21
iter46 errorY= 0.3189, prior information ratio= 0.04, Bdiff= 0.000235, Bkappa=190.5;pos. col. U=21
iter47 errorY= 0.3188, prior information ratio= 0.05, Bdiff= 0.0002231, Bkappa=201;pos. col. U=20
iter48 errorY= 0.3187, prior information ratio= 0.05, Bdiff= 0.0002111, Bkappa=207.8;pos. col. U=21
iter49 errorY= 0.3187, prior information ratio= 0.05, Bdiff= 0.0002003, Bkappa=214.1;pos. col. U=21
iter50 errorY= 0.3186, prior information ratio= 0.05, Bdiff= 0.0001911, Bkappa=219.8;pos. col. U=21
iter51 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 0.0001837, Bkappa=225.1;pos. col. U=21
iter52 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 0.0001764, Bkappa=229.6;pos. col. U=21
iter53 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 0.0001699, Bkappa=233.4;pos. col. U=21
iter54 errorY= 0.3185, prior information ratio= 0.04, Bdiff= 0.0001642, Bkappa=235.6;pos. col. U=21
iter55 errorY= 0.3184, prior information ratio= 0.05, Bdiff= 0.0001591, Bkappa=237.4;pos. col. U=20
iter56 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 0.0001543, Bkappa=237.8;pos. col. U=20
iter57 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 0.0001498, Bkappa=237;pos. col. U=20
iter58 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 0.0001454, Bkappa=234.9;pos. col. U=20
iter59 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 0.0001409, Bkappa=231.6;pos. col. U=20
Updating L3, current fraction= 0.5882, target=0.7
L3 is set to 0.004297 in 2 iterations
iter60 errorY= 0.3183, prior information ratio= 0.04, Bdiff= 0.0001358, Bkappa=226.7;pos. col. U=24
iter61 errorY= 0.3182, prior information ratio= 0.04, Bdiff= 0.0001292, Bkappa=220.7;pos. col. U=24
iter62 errorY= 0.3182, prior information ratio= 0.05, Bdiff= 0.0001236, Bkappa=214.4;pos. col. U=24
iter63 errorY= 0.3182, prior information ratio= 0.05, Bdiff= 0.0001184, Bkappa=206.7;pos. col. U=24
iter64 errorY= 0.3182, prior information ratio= 0.05, Bdiff= 0.0001133, Bkappa=198.4;pos. col. U=24
iter65 errorY= 0.3182, prior information ratio= 0.05, Bdiff= 0.0001085, Bkappa=189.7;pos. col. U=24
iter66 errorY= 0.3183, prior information ratio= 0.05, Bdiff= 0.0001041, Bkappa=180.9;pos. col. U=23
iter67 errorY= 0.3183, prior information ratio= 0.05, Bdiff= 9.994e-05, Bkappa=170.9;pos. col. U=23
iter68 errorY= 0.3183, prior information ratio= 0.05, Bdiff= 9.651e-05, Bkappa=160.9;pos. col. U=23
iter69 errorY= 0.3184, prior information ratio= 0.05, Bdiff= 9.415e-05, Bkappa=153.3;pos. col. U=23
iter70 errorY= 0.3184, prior information ratio= 0.05, Bdiff= 9.123e-05, Bkappa=143.7;pos. col. U=23
iter71 errorY= 0.3184, prior information ratio= 0.05, Bdiff= 8.832e-05, Bkappa=134.5;pos. col. U=23
iter72 errorY= 0.3184, prior information ratio= 0.05, Bdiff= 8.572e-05, Bkappa=126.2;pos. col. U=23
iter73 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 8.337e-05, Bkappa=117.6;pos. col. U=23
iter74 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 8.131e-05, Bkappa=109.4;pos. col. U=23
iter75 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.944e-05, Bkappa=101.7;pos. col. U=23
iter76 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.849e-05, Bkappa=94.4;pos. col. U=23
iter77 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.716e-05, Bkappa=87.53;pos. col. U=23
iter78 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.62e-05, Bkappa=80.36;pos. col. U=23
iter79 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.519e-05, Bkappa=74.61;pos. col. U=23
Updating L3, current fraction= 0.6765, target=0.7
L3 not changed
iter80 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.415e-05, Bkappa=70.78;pos. col. U=23
iter81 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.268e-05, Bkappa=66.16;pos. col. U=23
iter82 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.167e-05, Bkappa=61.99;pos. col. U=23
iter83 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.082e-05, Bkappa=56.16;pos. col. U=24
iter84 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 7.022e-05, Bkappa=52.68;pos. col. U=24
iter85 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 6.987e-05, Bkappa=49.57;pos. col. U=24
iter86 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 6.967e-05, Bkappa=46.41;pos. col. U=24
iter87 errorY= 0.3185, prior information ratio= 0.05, Bdiff= 6.956e-05, Bkappa=44.04;pos. col. U=23
iter88 errorY= 0.3184, prior information ratio= 0.06, Bdiff= 7.023e-05, Bkappa=42.75;pos. col. U=23
iter89 errorY= 0.3184, prior information ratio= 0.06, Bdiff= 7.058e-05, Bkappa=44.27;pos. col. U=23
iter90 errorY= 0.3184, prior information ratio= 0.06, Bdiff= 7.069e-05, Bkappa=46.13;pos. col. U=23
iter91 errorY= 0.3184, prior information ratio= 0.06, Bdiff= 7.064e-05, Bkappa=44.97;pos. col. U=23
iter92 errorY= 0.3183, prior information ratio= 0.06, Bdiff= 7.055e-05, Bkappa=41.8;pos. col. U=22
iter93 errorY= 0.3183, prior information ratio= 0.06, Bdiff= 7.017e-05, Bkappa=40.92;pos. col. U=22
iter94 errorY= 0.3183, prior information ratio= 0.06, Bdiff= 6.969e-05, Bkappa=40.26;pos. col. U=22
iter95 errorY= 0.3183, prior information ratio= 0.06, Bdiff= 6.908e-05, Bkappa=39.74;pos. col. U=22
iter96 errorY= 0.3182, prior information ratio= 0.06, Bdiff= 6.844e-05, Bkappa=39.31;pos. col. U=22
iter97 errorY= 0.3182, prior information ratio= 0.06, Bdiff= 6.767e-05, Bkappa=34.23;pos. col. U=22
iter98 errorY= 0.3182, prior information ratio= 0.06, Bdiff= 6.696e-05, Bkappa=34.05;pos. col. U=22
iter99 errorY= 0.3182, prior information ratio= 0.06, Bdiff= 6.637e-05, Bkappa=38.6;pos. col. U=22
Updating L3, current fraction= 0.6471, target=0.7
L3 is set to 0.004031 in 4 iterations
iter100 errorY= 0.3181, prior information ratio= 0.06, Bdiff= 6.529e-05, Bkappa=38.55;pos. col. U=24
iter101 errorY= 0.318, prior information ratio= 0.06, Bdiff= 6.683e-05, Bkappa=38.54;pos. col. U=24
iter102 errorY= 0.318, prior information ratio= 0.06, Bdiff= 6.59e-05, Bkappa=38.54;pos. col. U=25
iter103 errorY= 0.318, prior information ratio= 0.06, Bdiff= 6.518e-05, Bkappa=38.62;pos. col. U=25
iter104 errorY= 0.3179, prior information ratio= 0.05, Bdiff= 6.561e-05, Bkappa=38.73;pos. col. U=25
iter105 errorY= 0.3179, prior information ratio= 0.05, Bdiff= 6.535e-05, Bkappa=38.86;pos. col. U=25
iter106 errorY= 0.3179, prior information ratio= 0.05, Bdiff= 6.505e-05, Bkappa=37.64;pos. col. U=25
iter107 errorY= 0.3179, prior information ratio= 0.05, Bdiff= 6.484e-05, Bkappa=37.68;pos. col. U=25
iter108 errorY= 0.3179, prior information ratio= 0.05, Bdiff= 6.482e-05, Bkappa=34.67;pos. col. U=25
iter109 errorY= 0.3178, prior information ratio= 0.05, Bdiff= 6.51e-05, Bkappa=34.74;pos. col. U=26
iter110 errorY= 0.3178, prior information ratio= 0.05, Bdiff= 6.535e-05, Bkappa=32.4;pos. col. U=26
iter111 errorY= 0.3178, prior information ratio= 0.05, Bdiff= 6.569e-05, Bkappa=32.26;pos. col. U=26
iter112 errorY= 0.3178, prior information ratio= 0.06, Bdiff= 6.596e-05, Bkappa=32.11;pos. col. U=26
iter113 errorY= 0.3178, prior information ratio= 0.05, Bdiff= 6.622e-05, Bkappa=31.97;pos. col. U=26
iter114 errorY= 0.3177, prior information ratio= 0.05, Bdiff= 6.61e-05, Bkappa=31.39;pos. col. U=26
iter115 errorY= 0.3177, prior information ratio= 0.05, Bdiff= 6.636e-05, Bkappa=31.36;pos. col. U=26
iter116 errorY= 0.3177, prior information ratio= 0.06, Bdiff= 6.695e-05, Bkappa=31.33;pos. col. U=25
iter117 errorY= 0.3177, prior information ratio= 0.05, Bdiff= 6.69e-05, Bkappa=31.32;pos. col. U=25
iter118 errorY= 0.3176, prior information ratio= 0.06, Bdiff= 6.685e-05, Bkappa=31.31;pos. col. U=25
iter119 errorY= 0.3176, prior information ratio= 0.06, Bdiff= 6.697e-05, Bkappa=31.48;pos. col. U=23
Updating L3, current fraction= 0.6765, target=0.7
L3 not changed
iter120 errorY= 0.3176, prior information ratio= 0.06, Bdiff= 6.719e-05, Bkappa=33.2;pos. col. U=23
iter121 errorY= 0.3176, prior information ratio= 0.06, Bdiff= 6.723e-05, Bkappa=32.89;pos. col. U=23
iter122 errorY= 0.3176, prior information ratio= 0.06, Bdiff= 6.703e-05, Bkappa=32.98;pos. col. U=23
iter123 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 6.694e-05, Bkappa=31.78;pos. col. U=23
iter124 errorY= 0.3175, prior information ratio= 0.05, Bdiff= 6.672e-05, Bkappa=43.35;pos. col. U=23
iter125 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 6.617e-05, Bkappa=43.33;pos. col. U=24
iter126 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 6.532e-05, Bkappa=43.29;pos. col. U=24
iter127 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 6.432e-05, Bkappa=43.23;pos. col. U=23
iter128 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 6.327e-05, Bkappa=42.23;pos. col. U=23
iter129 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 6.205e-05, Bkappa=42.16;pos. col. U=23
iter130 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 6.084e-05, Bkappa=42.08;pos. col. U=23
iter131 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 5.968e-05, Bkappa=42.83;pos. col. U=23
iter132 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 5.849e-05, Bkappa=42.74;pos. col. U=23
iter133 errorY= 0.3175, prior information ratio= 0.06, Bdiff= 5.728e-05, Bkappa=42.62;pos. col. U=23
iter134 errorY= 0.3174, prior information ratio= 0.05, Bdiff= 6.302e-05, Bkappa=39.48;pos. col. U=23
iter135 errorY= 0.3174, prior information ratio= 0.05, Bdiff= 6.486e-05, Bkappa=39.35;pos. col. U=23
iter136 errorY= 0.3174, prior information ratio= 0.05, Bdiff= 6.52e-05, Bkappa=39.21;pos. col. U=23
iter137 errorY= 0.3174, prior information ratio= 0.05, Bdiff= 6.48e-05, Bkappa=39.06;pos. col. U=23
iter138 errorY= 0.3174, prior information ratio= 0.05, Bdiff= 6.393e-05, Bkappa=38.9;pos. col. U=23
iter139 errorY= 0.3174, prior information ratio= 0.05, Bdiff= 6.293e-05, Bkappa=38.3;pos. col. U=23
Updating L3, current fraction= 0.6765, target=0.7
L3 not changed
iter140 errorY= 0.3175, prior information ratio= 0.05, Bdiff= 6.207e-05, Bkappa=38.24;pos. col. U=23
iter141 errorY= 0.3175, prior information ratio= 0.05, Bdiff= 6.134e-05, Bkappa=38.19;pos. col. U=23
iter142 errorY= 0.3175, prior information ratio= 0.05, Bdiff= 6.265e-05, Bkappa=38.14;pos. col. U=23
iter143 errorY= 0.3175, prior information ratio= 0.05, Bdiff= 6.208e-05, Bkappa=43.24;pos. col. U=23
iter144 errorY= 0.3175, prior information ratio= 0.05, Bdiff= 6.165e-05, Bkappa=43.11;pos. col. U=23
iter145 errorY= 0.3175, prior information ratio= 0.05, Bdiff= 6.125e-05, Bkappa=42.99;pos. col. U=23
iter146 errorY= 0.3176, prior information ratio= 0.05, Bdiff= 6.105e-05, Bkappa=42.41;pos. col. U=23
iter147 errorY= 0.3176, prior information ratio= 0.05, Bdiff= 6.087e-05, Bkappa=42.31;pos. col. U=23
iter148 errorY= 0.3176, prior information ratio= 0.04, Bdiff= 6.139e-05, Bkappa=42.26;pos. col. U=23
iter149 errorY= 0.3176, prior information ratio= 0.04, Bdiff= 6.263e-05, Bkappa=42.3;pos. col. U=23
iter150 errorY= 0.3176, prior information ratio= 0.04, Bdiff= 6.237e-05, Bkappa=42.35;pos. col. U=23
iter151 errorY= 0.3177, prior information ratio= 0.04, Bdiff= 6.224e-05, Bkappa=42.41;pos. col. U=23
iter152 errorY= 0.3177, prior information ratio= 0.04, Bdiff= 6.322e-05, Bkappa=38.18;pos. col. U=23
iter153 errorY= 0.3177, prior information ratio= 0.04, Bdiff= 6.573e-05, Bkappa=38.27;pos. col. U=23
Bdiff is not decreasing
iter154 errorY= 0.3177, prior information ratio= 0.04, Bdiff= 6.611e-05, Bkappa=40.64;pos. col. U=23
Bdiff is not decreasing
iter155 errorY= 0.3177, prior information ratio= 0.04, Bdiff= 6.637e-05, Bkappa=41.07;pos. col. U=23
Bdiff is not decreasing
iter156 errorY= 0.3178, prior information ratio= 0.04, Bdiff= 6.662e-05, Bkappa=41.52;pos. col. U=23
Bdiff is not decreasing
iter157 errorY= 0.3178, prior information ratio= 0.04, Bdiff= 6.675e-05, Bkappa=41.98;pos. col. U=23
Bdiff is not decreasing
iter158 errorY= 0.3178, prior information ratio= 0.04, Bdiff= 6.668e-05, Bkappa=42.44;pos. col. U=23
Bdiff is not decreasing
converged at  iteration 158 Bdiff is not decreasing
There are 15  LVs with AUC>0.70

Apply PLIER model to “new” expression dataset

The GetNewDataB function will first row-normalize and reorder the “new” input gene expression data (exprs.mat), and then using a previously computed PLIER model (plier.model, specifically the gene loadings and the L2 constant), get the new data into the PLIER model LV space. Here, we supply the NARES data as the expression data and the PLIER model that has already been trained on the same gene expression matrix.

new.b.mat <- GetNewDataB(exprs.mat = exprs.mat,
                         plier.model = nares.plier)
# NARES B matrix from PLIERNewData
nares.b.mat <- nares.plier$B

We want to ensure that the two B matrices are the same.

all.equal(nares.b.mat, new.b.mat)

[1] TRUE

Reconstruction of gene expression data

We reconstruct gene expression data from the gene loadings and LVs.

nares.recon <- GetReconstructedExprs(z.matrix = as.matrix(nares.plier$Z),
                                     b.matrix = as.matrix(nares.plier$B))

Let’s evaluate the reconstruction. We’ll need the row-normalized NARES expression data (input) for comparison.

nares.rownorm <- PLIER::rowNorm(exprs.mat)
nares.rownorm <- nares.rownorm[rownames(nares.recon), ]

Spearman correlation between input, row-normalized expression data and the reconstructed data. If correlation between the input and the reconstructed data is high, that suggests that reconstruction is “successful.” Given the different constraints in PLIER, we would not expect to perfectly (rho = 1) reconstruct the input data. This particular evaluation will be most useful when we look at applying a trained PLIER model to a test dataset.

recon.correlation <- GetReconstructionCorrelation(true.mat = nares.rownorm,
                                                  recon.mat = nares.recon)

Plot density

# density plot
p <- ggplot2::ggplot(as.data.frame(recon.correlation), 
                ggplot2::aes(x = recon.correlation)) + 
      ggplot2::geom_density() +
      ggplot2::theme_bw() +
      ggplot2::labs(x = "Spearman Correlation",
                    title = "Input vs. PLIER reconstructed NARES data") +
      ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5, 
                                                        face = "bold"))
p

png.file <- file.path("plots", "01",
                      "NARES_reconstructed_data_correlation_density.png")
ggplot2::ggsave(filename = png.file, plot = p, width = 7, height = 5,
                units = "in")