J. Taroni 2018

The NARES study (Grayson, et al. Arthritis Rheumatol. 2015) included three sets of patients with ANCA-associated vasculitis: those with active nasal disease, prior nasal disease, and those with no history of nasal disease. There is also a composite comparator group.

Here, we’ll test for PLIER LVs that are differentially expressed between groups. We’ll focus on the multiPLIER LVs (e.g., the recount2 PLIER model), as we have demonstrated the validity of this approach in prior analyses.

Functions and directory set up

`%>%` <- dplyr::`%>%` 
source(file.path("util", "test_LV_differences.R"))
# plot and result directory setup for this notebook
plot.dir <- file.path("plots", "18")
dir.create(plot.dir, recursive = TRUE, showWarnings = FALSE)
results.dir <- file.path("results", "18")
dir.create(results.dir, recursive = TRUE, showWarnings = FALSE)

Read in data

Clinical data

This clinical data includes information like measures of disease severity (BVAS) and disease duration. We’re most interested in the disease labels (Classification).

nares.demo <- readr::read_tsv(file.path("data", "sample_info", 
                                        "NARES_demographic_data.tsv")) %>%
  dplyr::mutate(Classification = dplyr::recode(Classification, 
                                               C1 = "Healthy",
                                               C2 = "Sarcoidosis",
                                               C3 = "Allergic Rhinitis",
                                               C4 = "EGPA",
                                               V1 = "GPA (no prior nasal disease)",
                                               V2 = "GPA with prior nasal disease",
                                               V3 = "GPA with active nasal disease"))
Parsed with column specification:
cols(
  .default = col_character(),
  Batch = col_integer(),
  Age = col_integer(),
  Smoking_pkyrs = col_integer(),
  Steroids_daily_pred_mg = col_double(),
  Steroids_Cat = col_integer(),
  GC_or_Immune = col_integer(),
  Immune_or_Nasal = col_integer(),
  Any_Immune = col_integer()
)
See spec(...) for full column specifications.

recount2 LVs

recount.b <- readRDS(file.path("results", "13", "NARES_recount2_B.RDS"))

NARES differential latent variable expression (DLVE)

We’re going to compare the seven classification groups, which correspond to different disease labels (and were recoded above).

nares.df.list <- LVTestWrapper(b.matrix = recount.b,
                                sample.info.df = nares.demo,
                                phenotype.col = "Classification",
                                file.lead = "NARES_recount2_model", 
                                plot.dir = plot.dir,
                                results.dir = results.dir)
Column `Sample` joining factor and character vector, coercing into character vector

Example LV

nares.limma.df <- nares.df.list$limma
nares.limma.df %>% 
  dplyr::filter(adj.P.Val < 0.05) %>% 
  dplyr::arrange(adj.P.Val)
# read in recount2 summary data.frame
summary.file <- file.path("results", "08", "recount2_PLIER_summary.tsv")
recount.summary <- readr::read_tsv(summary.file)
Parsed with column specification:
cols(
  pathway = col_character(),
  `LV index` = col_integer(),
  AUC = col_double(),
  `p-value` = col_double(),
  FDR = col_double()
)

Let’s plot LV96, as it looks biologically interesting because of it’s putative association with natural killer (NK) cells. We’d expect myeloid cell, specifically granulocyte, signatures to be differentially expressed in this dataset, but not necessarily lymphocyte signatures. First, does the name agree with the recount2 PLIER summary?

recount.summary %>%
  dplyr::filter(`LV index` == 96)

Yes, SVM NK cells activated is significantly associated with LV96.

recount.b.df <- nares.df.list$b.df
recount.b.df %>%
  dplyr::filter(LV == "96,SVM NK cells activated") %>%
  ggplot2::ggplot(ggplot2::aes(x = Classification, y = Value)) +
  ggplot2::geom_boxplot() +
  ggplot2::geom_jitter(position = ggplot2::position_jitter(0.15)) + 
  ggplot2::theme_bw() +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1),
                 plot.title = ggplot2::element_text(hjust = 0.5, face = "bold")) +
  ggplot2::ggtitle("LV 96 SVM NK cells activated")

LS0tCnRpdGxlOiAiTkFSRVMgZGlmZmVyZW50aWFsIGV4cHJlc3Npb24iCm91dHB1dDogICAKICBodG1sX25vdGVib29rOiAKICAgIHRvYzogdHJ1ZQogICAgdG9jX2Zsb2F0OiB0cnVlCi0tLQoKKipKLiBUYXJvbmkgMjAxOCoqCgpUaGUgTkFSRVMgc3R1ZHkgKFtHcmF5c29uLCBldCBhbC4gX0FydGhyaXRpcyBSaGV1bWF0b2xfLiAyMDE1XShodHRwczovL2R4LmRvaS5vcmcvMTAuMTAwMi9hcnQuMzkxODUpKSAKaW5jbHVkZWQgdGhyZWUgc2V0cyBvZiBwYXRpZW50cyB3aXRoIEFOQ0EtYXNzb2NpYXRlZCB2YXNjdWxpdGlzOiB0aG9zZSB3aXRoIAphY3RpdmUgbmFzYWwgZGlzZWFzZSwgcHJpb3IgbmFzYWwgZGlzZWFzZSwgYW5kIHRob3NlIHdpdGggbm8gaGlzdG9yeSBvZiBuYXNhbCAKZGlzZWFzZS4KVGhlcmUgaXMgYWxzbyBhIGNvbXBvc2l0ZSBjb21wYXJhdG9yIGdyb3VwLgoKSGVyZSwgd2UnbGwgdGVzdCBmb3IgUExJRVIgTFZzIHRoYXQgYXJlIGRpZmZlcmVudGlhbGx5IGV4cHJlc3NlZCBiZXR3ZWVuIGdyb3Vwcy4KV2UnbGwgZm9jdXMgb24gdGhlIG11bHRpUExJRVIgTFZzIChlLmcuLCB0aGUgcmVjb3VudDIgUExJRVIgbW9kZWwpLCBhcyAKd2UgaGF2ZSBkZW1vbnN0cmF0ZWQgdGhlIHZhbGlkaXR5IG9mIHRoaXMgYXBwcm9hY2ggaW4gcHJpb3IgYW5hbHlzZXMuCgojIyBGdW5jdGlvbnMgYW5kIGRpcmVjdG9yeSBzZXQgdXAKCmBgYHtyfQpgJT4lYCA8LSBkcGx5cjo6YCU+JWAgCnNvdXJjZShmaWxlLnBhdGgoInV0aWwiLCAidGVzdF9MVl9kaWZmZXJlbmNlcy5SIikpCmBgYAoKYGBge3J9CiMgcGxvdCBhbmQgcmVzdWx0IGRpcmVjdG9yeSBzZXR1cCBmb3IgdGhpcyBub3RlYm9vawpwbG90LmRpciA8LSBmaWxlLnBhdGgoInBsb3RzIiwgIjE4IikKZGlyLmNyZWF0ZShwbG90LmRpciwgcmVjdXJzaXZlID0gVFJVRSwgc2hvd1dhcm5pbmdzID0gRkFMU0UpCnJlc3VsdHMuZGlyIDwtIGZpbGUucGF0aCgicmVzdWx0cyIsICIxOCIpCmRpci5jcmVhdGUocmVzdWx0cy5kaXIsIHJlY3Vyc2l2ZSA9IFRSVUUsIHNob3dXYXJuaW5ncyA9IEZBTFNFKQpgYGAKCiMjIFJlYWQgaW4gZGF0YQoKIyMjIENsaW5pY2FsIGRhdGEKClRoaXMgY2xpbmljYWwgZGF0YSBpbmNsdWRlcyBpbmZvcm1hdGlvbiBsaWtlIG1lYXN1cmVzIG9mIGRpc2Vhc2Ugc2V2ZXJpdHkgCihbYEJWQVNgXShodHRwczovL2RvaS5vcmcvMTAuMTEzNi9hcmQuMjAwOC4xMDEyNzkpKSBhbmQgZGlzZWFzZSBkdXJhdGlvbi4KV2UncmUgbW9zdCBpbnRlcmVzdGVkIGluIHRoZSBkaXNlYXNlIGxhYmVscyAoYENsYXNzaWZpY2F0aW9uYCkuCgpgYGB7cn0KbmFyZXMuZGVtbyA8LSByZWFkcjo6cmVhZF90c3YoZmlsZS5wYXRoKCJkYXRhIiwgInNhbXBsZV9pbmZvIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAiTkFSRVNfZGVtb2dyYXBoaWNfZGF0YS50c3YiKSkgJT4lCiAgZHBseXI6Om11dGF0ZShDbGFzc2lmaWNhdGlvbiA9IGRwbHlyOjpyZWNvZGUoQ2xhc3NpZmljYXRpb24sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIEMxID0gIkhlYWx0aHkiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIEMyID0gIlNhcmNvaWRvc2lzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBDMyA9ICJBbGxlcmdpYyBSaGluaXRpcyIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgQzQgPSAiRUdQQSIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgVjEgPSAiR1BBIChubyBwcmlvciBuYXNhbCBkaXNlYXNlKSIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgVjIgPSAiR1BBIHdpdGggcHJpb3IgbmFzYWwgZGlzZWFzZSIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgVjMgPSAiR1BBIHdpdGggYWN0aXZlIG5hc2FsIGRpc2Vhc2UiKSkKYGBgCgojIyMgcmVjb3VudDIgTFZzCgpgYGB7cn0KcmVjb3VudC5iIDwtIHJlYWRSRFMoZmlsZS5wYXRoKCJyZXN1bHRzIiwgIjEzIiwgIk5BUkVTX3JlY291bnQyX0IuUkRTIikpCmBgYAoKIyMgTkFSRVMgZGlmZmVyZW50aWFsIGxhdGVudCB2YXJpYWJsZSBleHByZXNzaW9uIChETFZFKQoKV2UncmUgZ29pbmcgdG8gY29tcGFyZSB0aGUgc2V2ZW4gY2xhc3NpZmljYXRpb24gZ3JvdXBzLCB3aGljaCBjb3JyZXNwb25kIHRvCmRpZmZlcmVudCBkaXNlYXNlIGxhYmVscyAoYW5kIHdlcmUgcmVjb2RlZCBhYm92ZSkuCgpgYGB7cn0KbmFyZXMuZGYubGlzdCA8LSBMVlRlc3RXcmFwcGVyKGIubWF0cml4ID0gcmVjb3VudC5iLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc2FtcGxlLmluZm8uZGYgPSBuYXJlcy5kZW1vLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcGhlbm90eXBlLmNvbCA9ICJDbGFzc2lmaWNhdGlvbiIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmaWxlLmxlYWQgPSAiTkFSRVNfcmVjb3VudDJfbW9kZWwiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHBsb3QuZGlyID0gcGxvdC5kaXIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICByZXN1bHRzLmRpciA9IHJlc3VsdHMuZGlyKQpgYGAKCiMjIyBFeGFtcGxlIExWCgpgYGB7cn0KbmFyZXMubGltbWEuZGYgPC0gbmFyZXMuZGYubGlzdCRsaW1tYQpuYXJlcy5saW1tYS5kZiAlPiUgCiAgZHBseXI6OmZpbHRlcihhZGouUC5WYWwgPCAwLjA1KSAlPiUgCiAgZHBseXI6OmFycmFuZ2UoYWRqLlAuVmFsKQpgYGAKCmBgYHtyfQojIHJlYWQgaW4gcmVjb3VudDIgc3VtbWFyeSBkYXRhLmZyYW1lCnN1bW1hcnkuZmlsZSA8LSBmaWxlLnBhdGgoInJlc3VsdHMiLCAiMDgiLCAicmVjb3VudDJfUExJRVJfc3VtbWFyeS50c3YiKQpyZWNvdW50LnN1bW1hcnkgPC0gcmVhZHI6OnJlYWRfdHN2KHN1bW1hcnkuZmlsZSkKYGBgCgpMZXQncyBwbG90IGBMVjk2YCwgYXMgaXQgbG9va3MgYmlvbG9naWNhbGx5IGludGVyZXN0aW5nIGJlY2F1c2Ugb2YgaXQncyAKcHV0YXRpdmUgYXNzb2NpYXRpb24gd2l0aCBuYXR1cmFsIGtpbGxlciAoTkspIGNlbGxzLiAKV2UnZCBleHBlY3QgbXllbG9pZCBjZWxsLCBzcGVjaWZpY2FsbHkgZ3JhbnVsb2N5dGUsIHNpZ25hdHVyZXMgdG8gYmUgCmRpZmZlcmVudGlhbGx5IGV4cHJlc3NlZCBpbiB0aGlzIGRhdGFzZXQsIGJ1dCBub3QgbmVjZXNzYXJpbHkgbHltcGhvY3l0ZSAKc2lnbmF0dXJlcy4KRmlyc3QsIGRvZXMgdGhlIF9uYW1lXyBhZ3JlZSB3aXRoIHRoZSByZWNvdW50MiBQTElFUiBzdW1tYXJ5PwoKYGBge3J9CnJlY291bnQuc3VtbWFyeSAlPiUKICBkcGx5cjo6ZmlsdGVyKGBMViBpbmRleGAgPT0gOTYpCmBgYAoKWWVzLCBgU1ZNIE5LIGNlbGxzIGFjdGl2YXRlZGAgaXMgc2lnbmlmaWNhbnRseSBhc3NvY2lhdGVkIHdpdGggYExWOTZgLgoKYGBge3J9CnJlY291bnQuYi5kZiA8LSBuYXJlcy5kZi5saXN0JGIuZGYKcmVjb3VudC5iLmRmICU+JQogIGRwbHlyOjpmaWx0ZXIoTFYgPT0gIjk2LFNWTSBOSyBjZWxscyBhY3RpdmF0ZWQiKSAlPiUKICBnZ3Bsb3QyOjpnZ3Bsb3QoZ2dwbG90Mjo6YWVzKHggPSBDbGFzc2lmaWNhdGlvbiwgeSA9IFZhbHVlKSkgKwogIGdncGxvdDI6Omdlb21fYm94cGxvdCgpICsKICBnZ3Bsb3QyOjpnZW9tX2ppdHRlcihwb3NpdGlvbiA9IGdncGxvdDI6OnBvc2l0aW9uX2ppdHRlcigwLjE1KSkgKyAKICBnZ3Bsb3QyOjp0aGVtZV9idygpICsKICBnZ3Bsb3QyOjp0aGVtZShheGlzLnRleHQueCA9IGdncGxvdDI6OmVsZW1lbnRfdGV4dChhbmdsZSA9IDQ1LCBoanVzdCA9IDEpLAogICAgICAgICAgICAgICAgIHBsb3QudGl0bGUgPSBnZ3Bsb3QyOjplbGVtZW50X3RleHQoaGp1c3QgPSAwLjUsIGZhY2UgPSAiYm9sZCIpKSArCiAgZ2dwbG90Mjo6Z2d0aXRsZSgiTFYgOTYgU1ZNIE5LIGNlbGxzIGFjdGl2YXRlZCIpCmBgYAoK