Train a PLIER model on the NARES data

J. Taroni 2018

The NARES study was published in the following article:

Grayson PC, Steiling K, Platt M, et al. Defining the Nasal Transcriptome in Granulomatosis with Polyangiitis. Arthritis & Rheumatology, 2015. doi: 10.1002/art.39185.

We processed this data in the greenelab/rheum-data-plier repo.

This dataset (n = 76) contains nasal brushing samples from patients with granulomatosis with polyangiitis (GPA), a form of ANCA-associated vasculitis. Patients with GPA either had active, a prior history of, or no history of nasal disease. A composite comparator group of patients with sarcoidosis, patients with allergic rhinitis and healthy controls is also included.

Here we’ll train a PLIER model and explore:

• The sparsity of the U matrix
• How many LVs are recovered? How many pathways?
• What are the pathways that are represented?

Libraries

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

NARES data

nares.data <- readr::read_tsv(file.path("data", "expression_data", 
                                         "NARES_SCANfast_ComBat.pcl"))

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

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")
symbol.file <- 
  file.path("data", "expression_data", 
            "NARES_SCANfast_ComBat_with_GeneSymbol.pcl")
readr::write_delim(annot.nares.data, path = symbol.file, delim = "\t")
# 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))
rm(list = setdiff(ls(), c("exprs.mat")))

Functions and directory set up

Organized this way for convenience w.r.t. removing objects from the environment after gene identifier conversion.

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

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

PLIER model

Training

# train model and save to file
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

plier.file <- file.path(results.dir, "NARES_PLIER_model.RDS")
saveRDS(nares.plier, file = plier.file)
# save summary data.frame to results
summary.file <- file.path(results.dir, "NARES_summary.tsv")
readr::write_tsv(nares.plier$summary, path = summary.file)

Number of LVs

nrow(nares.plier$B)

[1] 34

The model learns 34 latent variables. From the trace during training, we know that 15 of them have AUC>0.70.

Examine U and pathway coverage

As noted in the SLE WB PLIER notebook:

U is the prior information coefficient matrix; it tells us how the prior information in the form of pathways/gene sets relates to LVs.

PLIER::plotU(nares.plier, fontsize_row = 7)

To summarize some of these results:

• A latent variable associated with extracellular matrix formation is captured. LV14 is associated with REACTOME_EXTRACELLULAR_MATRIX_ORGANIZATION, REACTOME_COLLAGEN_FORMATION, and PID_INTEGRIN1_PATHWAY.
• There’s a neutrophil-associated LV, LV3.
• Interferon signaling is probably captured by LV25, but two immune cell types, M1 macrophages and CD8 T cells, are also associated with this LV. This is likely due to their production of or activation by IFN and we saw a bit of this overlap with cell types in other models (SLE WB, recount2).

# save U plot to file
pdf(file.path(plot.dir, "NARES_U_plot.pdf"), height = 5, width = 7)
PLIER::plotU(nares.plier, fontsize_row = 8, fontsize_col = 8)
dev.off()

null device 
          1 

What proportion of the pathways input into the PLIER model are significantly associated with at least one LV?

coverage.results <- GetPathwayCoverage(nares.plier)  # FDR < 0.05 by default
coverage.results$pathway

[1] 0.09394904

Less than 10% of the pathways we are interested in are significantly associated with a latent variable.

How sparse is U?

Take into account all pathways

u.sparsity.all <- CalculateUSparsity(nares.plier,
                                     significant.only = FALSE)
ggplot2::ggplot(as.data.frame(u.sparsity.all),
                ggplot2::aes(x = u.sparsity.all)) +
  ggplot2::geom_density(fill = "blue", alpha = 0.5) +
  ggplot2::theme_bw() +
  ggplot2::labs(x = "proportion of positive entries in U",
                title = "All LVs") +
  ggplot2::theme(text = ggplot2::element_text(size = 15))

summary(u.sparsity.all)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000000 0.000000 0.004777 0.005761 0.011943 0.015924 

plot.file <- file.path(plot.dir, "NARES_U_sparsity_all.png")
ggplot2::ggsave(plot.file, plot = ggplot2::last_plot())

Saving 7 x 7 in image

What proportion of entries in the U matrix for each LV are significantly associated with that LV?

u.sparsity.sig <- CalculateUSparsity(nares.plier,
                                     significant.only = TRUE)
ggplot2::ggplot(as.data.frame(u.sparsity.sig),
                ggplot2::aes(x = u.sparsity.sig)) +
  ggplot2::geom_density(fill = "blue", alpha = 0.5) +
  ggplot2::theme_bw() +
  ggplot2::labs(x = "proportion of positive entries in U",
                title = "Pathway-associated LVs")  +
  ggplot2::theme(text = ggplot2::element_text(size = 15))

summary(u.sparsity.sig)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000000 0.000000 0.001592 0.003185 0.004777 0.012739 

plot.file <- file.path(plot.dir, "NARES_U_sparsity_significant.png")
ggplot2::ggsave(plot.file, plot = ggplot2::last_plot())

Saving 7 x 7 in image