Guide to the effect of felzartamab anti-CD38 treatment on the molecular phenotype of antibody-mediated rejection in kidney transplant biopsies.
Patrick T Gauthier
2024-09-08
0.1 Foreword
This guide is meant to assist the interpretation of analyses carried out within “Diebold and Gauthier et al. Effect of felzartamab anti-CD38 treatment on the molecular phenotype of antibody-mediated rejection in kidney transplant biopsies. under review. Nature Medicine”. This manuscript is currently under peer-review and has not yet been accepted for publication. A pre-print of the manuscript is available for review via Research Sqaure (https://www.researchsquare.com/article/rs-4870641/v1). The manuscript describes a comprehensive transcriptomics assessment of the effect of felzartamab on gene expression in kidney allograft biopsies taken from patients with antibody-mediated rejection (ABMR). The work described herein was part of a randomized controlled trial (ClinicalTrials.gov; NCT05021484). Biospies were taken at baseline, as well as 24 and 52 weeks post-baseline. The efficacy of felzartamab in resolving antibody-mediated rejection was compared to the progression of molecular activity in biopsies collected from patients receiving a placebo. This guide is a companion to the GitHub repository (https://github.com/TSI-PTG/CD38-effect-of-treatment).
1 The effect of treatment on molecular scores
1.1 INTRODUCTION
Transcriptomics
assessment of biopsies from kidney allografts enables the interrogation
of a multitude of endpoints beyond individual gene expression. For
example, multiple gene expression-based machine learning classifiers
have been trained to predict the probability of histology lesion scores
indicative of antibody-mediated rejection (ABMR), T cell-mediated
rejection (TCMR), and overall ABMR and TCMR diagnoses. Predicted
probabilities from these classifiers can be used as a continuous
variable to assess the molecular activity of ABMR and TCMR. There are
also numerous gene lists that have been annotated in the transplantation
space to describe molecular profiles of rejection, cellular activity,
and responses to injury (https://www.ualberta.ca/en/medicine/institutes-centres-groups/atagc/research/gene-lists.html).
These gene lists have been translated into pathogenesis-based transcript
set (PBT) scores, and have been widely applied to depict gradients of
molecular activity related to allograft rejection and injury.
Supplementary Table 3. Molecular scores (PBTs, classifiers) | ||
|---|---|---|
Category | Abbreviation | Description |
TCMR-related | TCMR-RAT | TCMR-associated (1) |
QCAT | Cytotoxic T cell associated (2) | |
TCB | T-cell burden related (3) | |
i>1Prob | Top 20 transcripts in classifier of probability of histologic kidney i-lesion score > 1 (4) | |
t>1Prob | Top 20 transcripts in classifier of probability of histologic kidney t-lesion score > 1 (4) | |
TCMRProb | Top 20 transcripts in classifier of probability of histologic kidney TCMR diagnosis (4) | |
ABMR-related | NKB | Natural killer (NK) cell burden (5) (only valid if there is no TCMR) |
DSAST | Kidney derived DSA selective transcripts (3) | |
ABMRProb | Top 20 transcripts in classifier of probability of histologic kidney ABMR diagnosis (4) | |
g>0Prob | Top 20 transcripts in classifier of probability of histologic g-lesion score > 0 (4) | |
ptc>0Prob | Top 20 transcripts in classifier of probability of histologic ptc-lesion score > 0 (4) | |
New ABMR-related PBTs | AAG | ABMR-associated activity genes |
IIAAG | IFNG-inducible ABMR-associated activity genes | |
NKAAG | NK cell expressed ABMR-associated activity genes | |
AEG | ABMR-associated endothelial cell genes | |
Macrophage-related | AMAT1 | Alternatively activated macrophage transcripts (6) |
QCMAT | Constitutive macrophage transcripts (7) | |
Recent injury-related | IRITD3 | Injury-repair induced, day 3 (IRITD3) (8) |
IRITD5 | Injury-repair induced, day 5 (IRITD3) (8) | |
IRRAT30 | Injury-repair associated (IRRAT30) (9) | |
Atrophy-fibrosis-related | ci>1Prob | Fibrosis classifier (10) |
ct>1Prob | Atrophy classifier (10) | |
Normal parenchyma | KT1 | Normal kidney transcripts (except solute carriers) (11) |
KT2 | Normal kidney solute carrier transcripts (11) | |
A https://www.ualberta.ca/medicine/institutes-centres-groups/atagc/research/gene-lists. | ||
Reference List | ||
1. Halloran, P.F., Venner, J.M. & Famulski, K.S. Comprehensive Analysis of Transcript Changes Associated With Allograft Rejection: Combining Universal and Selective Features. Am J Transplant 17, 1754-1769 (2017). | ||
2. Hidalgo, L.G., et al. The transcriptome of human cytotoxic T cells: measuring the burden of CTL-associated transcripts in human kidney transplants. Am J Transplant 8, 637-646 (2008). | ||
3. Hidalgo, L.G., et al. NK cell transcripts and NK cells in kidney biopsies from patients with donor-specific antibodies: evidence for NK cell involvement in antibody-mediated rejection. Am J Transplant 10, 1812-1822 (2010). | ||
4. Reeve, J., et al. Assessing rejection-related disease in kidney transplant biopsies based on archetypal analysis of molecular phenotypes. Jci Insight 2, e94197 (2017). | ||
5. Hidalgo, L.G., et al. Interpreting NK cell transcripts versus T cell transcripts in renal transplant biopsies. Am J Transplant 12, 1180-1191 (2012). | ||
6. Famulski, K.S., Sis, B., Billesberger, L. & Halloran, P.F. Interferon-gamma and donor MHC class I control alternative macrophage activation and activin expression in rejecting kidney allografts: a shift in the Th1-Th2 paradigm. Am J Transplant 8, 547-556 (2008). | ||
7. Famulski, K.S., et al. Defining the canonical form of T-cell-mediated rejection in human kidney transplants. Am J Transplant 10, 810-820 (2010). | ||
8. Famulski, K.S., et al. Transcriptome analysis reveals heterogeneity in the injury response of kidney transplants. Am J Transplant 7, 2483-2495 (2007). | ||
9. Famulski, K.S., et al. Molecular phenotypes of acute kidney injury in kidney transplants. J Am Soc Nephrol 23, 948-958 (2012). | ||
10. Halloran, P.F., et al. Molecular phenotype of kidney transplant indication biopsies with inflammation in scarred areas. Am J Transplant 19, 1356-1370 (2019). | ||
11. Einecke, G., Broderick, G., Sis, B. & Halloran, P.F. Early loss of renal transcripts in kidney allografts: relationship to the development of histologic lesions and alloimmune effector mechanisms. Am J Transplant 7, 1121-1130 (2007). | ||
Abbreviations: pathogenesis-based transcript set (PBT); donor-derived cell-free DNA (dd-cfDNA); antibody-mediated rejection (ABMR); T cell mediated rejection (TCMR) | ||
A primary
analysis in the present work assessed how felzartamab therapy affected
classifier and PBT scores in patients with ABMR. The study design
included two arms, placebo and felzartamab. Patients from both arms were
biopsied prior to receiving placebo or felzartamab, and followed up with
two additional biopsies at 24 and 52 weeks. The goal of the present work
was to determine if felzartamab suppressed the molecular activity of
ABMR. This effect, the effect of treatment, was defined as the
interactive effect between treatment (placebo and felzartamab) and
follow up (baseline, week 24, and week 52). This design was chosen to
account for the confounding effects of time between placebo and
felzartamab. While it is statistically straightforward to measure
interactive effects using parametric models, the assumptions of
normality and heterogeneity of variance were not met across all
molecular scores, and thus, for consistency we simply chose to carry out
non-parametric tests for all scores. This presented a challenge as
standard non-parametric rank-based tests do no allow for testing
interactions among independent variables. Thus, we availed to
aligned-rank transformation, a technique that rank-orders independent
variables while maintaining rankings across independent variables. An
analysis of variance (ANOVA) can then be carried out on the aligned-rank
transformed data to allow for a fully factorial non-parametric ANOVA.
This is powerful in this context because it allows us to test for the
interactive effect of treatment (i.e., placebo vs felzartamab) and
follow-up (i.e., change in scores from biopsy to biopsy). Moreover,
aligned-rank transformation also allows for us to specify random effects
terms in the ANOVA to handle repeated measurements (i.e., biopsies)
across patients.
References:
• Halloran, P.F., Famulski, K.S. & Reeve, J. Molecular assessment of disease states in kidney transplant biopsy samples. Nat Rev Nephrol 12, 534-548 (2016).
• Wobbrock J, F.L., Gergle D, Higgins J The Aligned Rank Transform for Nonparametric Factorial Analyses Using Only ANOVA Procedures. In Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ’11) 143–146. doi:10.1145/1978942.1978963, https://depts.washington.edu/acelab/proj/art/.(2011).
• Kay M, E.L., Higgins J, Wobbrock J. ARTool: Aligned Rank Transform for Nonparametric Factorial ANOVAs. doi:10.5281/zenodo.594511, R package version 0.11.1, https://github.com/mjskay/ARTool. (2021).
• Elkin, L.A., Kay, M., Higgins, J.J. & Wobbrock, J.O. An Aligned Rank Transform Procedure for Multifactor Contrast Tests. in The 34th Annual ACM Symposium on User Interface Software and Technology 754–768 (Association for Computing Machinery, Virtual Event, USA, 2021).
1.2 DATA WRANGLING
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# Load necessary libraries
library(tidyverse)
library(broom)
library(rstatix)
library(flextable)
library(officer)
library(emmeans)
library(multcomp)
library(ARTool)
library(rcompanion)
library(Biobase)
library(ggpubr)
library(patchwork)
library(ggprism)
# Custom operators
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Suppress dplyr reframe info
options(dplyr.reframe.inform = FALSE)
# Load data
load("C:/R/CD38-effect-of-treatment/data/cd38.RData")
load("C:/R/CD38-effect-of-treatment/results/flextable_artANOVA.RData")
# load plot data
load("C:/R/CD38-effect-of-treatment/results/plots_artANOVA.RData")
# Seed for reproducibility
set.seed(42)
Now we
define the categories of molecular scores for easier organization
downstream.
# Feature categories
vars_cfDNA <- c("cfDNA_cpml")
vars_abmr <- c("ABMRpm", "ggt0", "ptcgt0", "DSAST", "NKB")
vars_tcmr <- c("TCMRt", "tgt1", "igt1", "QCAT", "TCB")
vars_injury <- c("IRRAT30", "IRITD3", "IRITD5", "cigt1", "ctgt1")
vars_parenchyma <- c("KT1", "KT2")
vars_macrophage <- c("AMAT1", "QCMAT")
# DEFINE VARIABLES TO ASSESS ####
vars <- c(vars_cfDNA, vars_abmr, vars_tcmr, vars_macrophage, vars_injury, vars_parenchyma)
Now we
process the data so we can iteratively carry out the analyses on each
score.
# DEFINE THE DATA ####
data <- set %>%
pData() %>%
dplyr::select(Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, all_of(vars)) %>%
dplyr::filter(Patient %nin% c(15, 18)) %>%
left_join(., summarise(., sample_pairs = n(), .by = Felzartamab_Followup), by = "Felzartamab_Followup") %>%
relocate(sample_pairs, .after = Felzartamab_Followup) %>%
arrange(Felzartamab, Patient, Group)
# WRANGLE THE PHENOTYPE DATA ####
data_00 <- data %>%
expand_grid(category = c("cfDNA", "ABMR", "TCMR", "macrophage", "injury", "parenchyma")) %>%
nest(.by = category) %>%
mutate(
features = map(
category,
function(category) {
if (category == "cfDNA") {
vars_cfDNA
} else if (category == "ABMR") {
vars_abmr
} else if (category == "TCMR") {
vars_tcmr
} else if (category == "macrophage") {
vars_macrophage
} else if (category == "injury") {
vars_injury
} else if (category == "parenchyma") {
vars_parenchyma
}
}
),
data = pmap(
list(features, data),
function(features, data) {
data %>%
dplyr::select(
Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, sample_pairs,
all_of(features)
)
}
)
)
# WRANGLE THE DATA FOR UNIVARIATE TESTS ####
data_01 <- data_00 %>%
mutate(
data_univariate = pmap(
list(data, features),
function(data, features) {
data %>%
pivot_longer(
cols = all_of(features),
names_to = "variable",
values_to = "value"
) %>%
nest(.by = variable) %>%
mutate(
variable = variable %>%
factor(
levels = vars
),
annotation = case_when(
variable %in% vars_cfDNA ~ "cfDNA",
variable %in% vars_tcmr ~ "TCMR-related",
variable %in% vars_abmr ~ "ABMR-related",
variable %in% vars_macrophage ~ "macrophage-related",
variable %in% vars_injury ~ "injury-related",
variable %in% vars_parenchyma ~ "parenchyma-related",
TRUE ~ " "
) %>%
factor(
levels = c(
"cfDNA",
"ABMR-related",
"TCMR-related",
"macrophage-related",
"injury-related",
"parenchyma-related",
"archetypes",
"rejection PC",
"injury PC"
)
),
score = case_when(
variable == "cfDNA_cpml" ~ "Donor-derived cell-free DNA (dd-cfDNA, cp/mL)",
variable == "TCMRt" ~ "TCMR classifier (TCMRProb)",
variable == "TCB" ~ "T cell burden (TCB)",
variable == "tgt1" ~ "Tubulitis classifier (t>1Prob)",
variable == "igt1" ~ "Interstitial infiltrate classifier (i>1Prob)",
variable == "QCAT" ~ "Cytotoxic T cell-associated (QCAT)",
variable == "ABMRpm" ~ "ABMR classifier (ABMRProb)",
variable == "DSAST" ~ "DSA-selective transcripts (DSAST)",
variable == "NKB" ~ "NK cell burden (NKB)",
variable == "ggt0" ~ "Glomerulitis classifier (g>0Prob)",
variable == "ptcgt0" ~ "Peritubular capillaritis classifier (ptc>0Prob)",
variable == "AMAT1" ~ "Alternatively activated macrophage (AMAT1)",
variable == "QCMAT" ~ "Constitutive macrophage (QCMAT)",
variable == "IRITD3" ~ "Injury-repair induced, day 3 (IRITD3)",
variable == "IRITD5" ~ "Injury-repair induced, day 5 (IRITD5)",
variable == "IRRAT30" ~ "Injury-repair associated (IRRAT30)",
variable == "cigt1" ~ "Fibrosis classifier (ci>1Prob)",
variable == "ctgt1" ~ "Atrophy classifier (ct>1Prob)",
variable == "KT1" ~ "Kidney parenchymal (KT1)",
variable == "KT2" ~ "Kidney parenchymal - no solute carriers (KT2)"
), .before = 1
) %>%
arrange(annotation, variable)
}
)
) %>%
dplyr::select(category, data_univariate) %>%
unnest(data_univariate) %>%
mutate(n_cat = score %>% unique() %>% length(), .by = category, .after = variable)
This is
what the data look like after we have organized them. Each row of the
dataframe contains the data for each score, the category it belongs to,
its category annotation, its long-format name, and the name of the
variable in the data.
data_01 %>% print(n="all")## # A tibble: 20 x 6
## category annotation score variable
## <chr> <fct> <chr> <fct>
## 1 cfDNA cfDNA Dono~ cfDNA_c~
## 2 ABMR ABMR-related ABMR~ ABMRpm
## 3 ABMR ABMR-related Glom~ ggt0
## 4 ABMR ABMR-related Peri~ ptcgt0
## 5 ABMR ABMR-related DSA-~ DSAST
## 6 ABMR ABMR-related NK c~ NKB
## 7 TCMR TCMR-related TCMR~ TCMRt
## 8 TCMR TCMR-related Tubu~ tgt1
## 9 TCMR TCMR-related Inte~ igt1
## 10 TCMR TCMR-related Cyto~ QCAT
## 11 TCMR TCMR-related T ce~ TCB
## 12 macrophage macrophage-~ Alte~ AMAT1
## 13 macrophage macrophage-~ Cons~ QCMAT
## 14 injury injury-rela~ Inju~ IRRAT30
## 15 injury injury-rela~ Inju~ IRITD3
## 16 injury injury-rela~ Inju~ IRITD5
## 17 injury injury-rela~ Fibr~ cigt1
## 18 injury injury-rela~ Atro~ ctgt1
## 19 parenchyma parenchyma-~ Kidn~ KT1
## 20 parenchyma parenchyma-~ Kidn~ KT2
## # i 2 more variables: n_cat <int>,
## # data <list>1.3 MEDIAN SCORES AND MEDIAN EFFECT OF TREATMENT
At this
stage we can summarize the median values for each score by treatment and
follow-up.
(Note that the median scores are used for general interpretation of
treatment effects because aligned-ranks, the data being assessed for
inference regarding the effects of treatment, are not easily
interpretable.)
# UNIVARIATE MEDIANS ####
data_02 <- data_01 %>%
mutate(
medians = map(
data,
function(data) {
data %>%
reframe(
median = value %>% median(),
IQR = value %>% IQR(),
sample_pairs = n(),
.by = c(Followup, Felzartamab, Felzartamab_Followup)
) %>%
arrange(Felzartamab_Followup)
}
),
medians_delta = map(
medians,
function(medians) {
medians %>%
reframe(
median_delta = combn(median, 2, diff) %>% as.numeric(),
IQR_delta = combn(IQR, 2, mean) %>% as.numeric(),
sample_pairs = sample_pairs,
.by = c(Felzartamab)
) %>%
mutate(
Followup_pairwise = rep(c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52"), 2) %>%
factor(levels = c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52")),
.before = sample_pairs
) %>%
mutate(
median_delta_delta = combn(median_delta, 2, diff) %>% as.numeric(),
IQR_delta_delta = combn(IQR_delta, 2, mean) %>% as.numeric(),
.by = c(Followup_pairwise),
.before = sample_pairs
)
}
)
)
Let’s
take a look at what the median summaries look like for the ABMRprob
score.
data_02 %>%
dplyr::filter(variable == "ABMRpm") %>%
pull(medians) %>%
pluck(1) %>%
flextable::flextable()Followup | Felzartamab | Felzartamab_Followup | median | IQR | sample_pairs |
|---|---|---|---|---|---|
Baseline | Placebo | Baseline_Placebo | 0.6695 | 0.35525 | 10 |
Week24 | Placebo | Week24_Placebo | 0.7740 | 0.39025 | 10 |
Week52 | Placebo | Week52_Placebo | 0.5185 | 0.49825 | 10 |
Baseline | Felzartamab | Baseline_Felzartamab | 0.7360 | 0.30775 | 10 |
Week24 | Felzartamab | Week24_Felzartamab | 0.1675 | 0.35700 | 10 |
Week52 | Felzartamab | Week52_Felzartamab | 0.7025 | 0.53875 | 10 |
Now we
can look at the change in median ABMRprob scores in the placebo and
felzartamab arms, as well as the difference in those medians (i.e., ΔΔ
median - this is our summarized treatment effect).
data_02 %>%
dplyr::filter(variable == "ABMRpm") %>%
pull(medians_delta) %>%
pluck(1) %>%
flextable::flextable()Felzartamab | median_delta | IQR_delta | Followup_pairwise | median_delta_delta | IQR_delta_delta | sample_pairs |
|---|---|---|---|---|---|---|
Placebo | 0.1045 | 0.372750 | Baseline - Week24 | -0.6730 | 0.3525625 | 10 |
Placebo | -0.1510 | 0.426750 | Baseline - Week52 | 0.1175 | 0.4250000 | 10 |
Placebo | -0.2555 | 0.444250 | Week24 - Week52 | 0.7905 | 0.4460625 | 10 |
Felzartamab | -0.5685 | 0.332375 | Baseline - Week24 | -0.6730 | 0.3525625 | 10 |
Felzartamab | -0.0335 | 0.423250 | Baseline - Week52 | 0.1175 | 0.4250000 | 10 |
Felzartamab | 0.5350 | 0.447875 | Week24 - Week52 | 0.7905 | 0.4460625 | 10 |
1.4 FIT THE MODELS
Now we
carry out the analysis. This includes testing the specific contrasts for
the interaction of treatment and time (i.e., our effect of
treatment).
# UNIVARIATE NONPARAMETRIC TESTS ####
artANOVA <- data_02 %>%
mutate(
art = map(data, art, formula = value ~ Followup * Felzartamab + (1 | Patient)),
art_aov = map(art, anova, type = "II"),
art_aov_tidy = map(art_aov, tidy) %>% suppressWarnings(),
art_aov_contrast = map(
art,
art.con,
formula = "Followup:Felzartamab", adjust = "fdr", method = "pairwise", interaction = TRUE, response = "art"
),
art_aov_contrast_tidy = map(art_aov_contrast, tidy),
art_aov_contrast_cld = map(
art_aov_contrast_tidy,
function(art_aov_contrast_tidy) {
art_aov_contrast_tidy %>%
as.data.frame() %>%
cldList(adj.p.value ~ Followup:Felzartamab, data = .)
}
)
)1.5 SUMMARIZE RESULTS
The
artANOVA effects are observed to reject or confirm a treatment effect on
the interactive effect for each score (i.e., to support the
interpretation of post-hoc tests for effect of
treatment).
# CREATE FLEXTABLE SUMMARY OF ART MODELS ####
title_art <- paste("Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients")
artANOVA %>%
dplyr::select(annotation, n_cat, score, art_aov) %>%
unnest(everything()) %>%
dplyr::rename(p.value = `Pr(>F)`) %>%
mutate(FDR = p.value %>% p.adjust(method = "fdr"), .by = annotation) %>%
dplyr::select(annotation, score, Term, `F`, p.value, FDR) %>%
flextable::flextable() %>%
flextable::add_header_row(values = rep(title_art, ncol_keys(.))) %>%
flextable::merge_h(part = "header") %>%
flextable::merge_v(j = 1:2) %>%
flextable::fontsize(size = 8, part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::bg(i = ~ FDR < 0.05 & Term == "Followup:Felzartamab", j = 3:6, bg = "#fbff00") %>%
flextable::colformat_double(j = 2:3, digits = 2) %>%
flextable::colformat_double(j = 4:ncol_keys(.), digits = 3) %>%
flextable::border_remove() %>%
flextable::bold(part = "header") %>%
flextable::padding(padding = 0, part = "all") %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::autofit()Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients | |||||
|---|---|---|---|---|---|
annotation | score | Term | F | p.value | FDR |
cfDNA | Donor-derived cell-free DNA (dd-cfDNA, cp/mL) | Followup | 4.418 | 0.019 | 0.029 |
Felzartamab | 1.598 | 0.222 | 0.222 | ||
Followup:Felzartamab | 10.396 | 0.000 | 0.001 | ||
ABMR-related | ABMR classifier (ABMRProb) | Followup | 6.659 | 0.003 | 0.007 |
Felzartamab | 1.620 | 0.219 | 0.235 | ||
Followup:Felzartamab | 9.605 | 0.000 | 0.001 | ||
Glomerulitis classifier (g>0Prob) | Followup | 11.768 | 0.000 | 0.001 | |
Felzartamab | 1.511 | 0.235 | 0.235 | ||
Followup:Felzartamab | 10.743 | 0.000 | 0.001 | ||
Peritubular capillaritis classifier (ptc>0Prob) | Followup | 14.188 | 0.000 | 0.000 | |
Felzartamab | 3.178 | 0.092 | 0.125 | ||
Followup:Felzartamab | 7.749 | 0.002 | 0.004 | ||
DSA-selective transcripts (DSAST) | Followup | 10.508 | 0.000 | 0.001 | |
Felzartamab | 1.747 | 0.203 | 0.234 | ||
Followup:Felzartamab | 4.179 | 0.023 | 0.035 | ||
NK cell burden (NKB) | Followup | 6.512 | 0.004 | 0.007 | |
Felzartamab | 1.956 | 0.179 | 0.224 | ||
Followup:Felzartamab | 5.114 | 0.011 | 0.018 | ||
TCMR-related | TCMR classifier (TCMRProb) | Followup | 3.555 | 0.039 | 0.165 |
Felzartamab | 1.417 | 0.249 | 0.468 | ||
Followup:Felzartamab | 1.124 | 0.336 | 0.504 | ||
Tubulitis classifier (t>1Prob) | Followup | 1.016 | 0.372 | 0.506 | |
Felzartamab | 1.853 | 0.190 | 0.408 | ||
Followup:Felzartamab | 0.065 | 0.937 | 0.968 | ||
Interstitial infiltrate classifier (i>1Prob) | Followup | 3.414 | 0.044 | 0.165 | |
Felzartamab | 0.997 | 0.331 | 0.504 | ||
Followup:Felzartamab | 0.032 | 0.968 | 0.968 | ||
Cytotoxic T cell-associated (QCAT) | Followup | 5.693 | 0.007 | 0.107 | |
Felzartamab | 4.118 | 0.057 | 0.172 | ||
Followup:Felzartamab | 0.927 | 0.405 | 0.506 | ||
T cell burden (TCB) | Followup | 4.360 | 0.020 | 0.151 | |
Felzartamab | 2.868 | 0.108 | 0.269 | ||
Followup:Felzartamab | 0.301 | 0.742 | 0.856 | ||
macrophage-related | Alternatively activated macrophage (AMAT1) | Followup | 1.549 | 0.226 | 0.304 |
Felzartamab | 1.047 | 0.320 | 0.320 | ||
Followup:Felzartamab | 2.406 | 0.105 | 0.209 | ||
Constitutive macrophage (QCMAT) | Followup | 3.378 | 0.045 | 0.209 | |
Felzartamab | 3.100 | 0.095 | 0.209 | ||
Followup:Felzartamab | 1.426 | 0.253 | 0.304 | ||
injury-related | Injury-repair associated (IRRAT30) | Followup | 0.055 | 0.946 | 0.946 |
Felzartamab | 3.842 | 0.066 | 0.223 | ||
Followup:Felzartamab | 1.915 | 0.162 | 0.347 | ||
Injury-repair induced, day 3 (IRITD3) | Followup | 0.174 | 0.841 | 0.946 | |
Felzartamab | 3.268 | 0.087 | 0.223 | ||
Followup:Felzartamab | 3.125 | 0.056 | 0.223 | ||
Injury-repair induced, day 5 (IRITD5) | Followup | 0.073 | 0.930 | 0.946 | |
Felzartamab | 1.088 | 0.311 | 0.583 | ||
Followup:Felzartamab | 2.586 | 0.089 | 0.223 | ||
Fibrosis classifier (ci>1Prob) | Followup | 0.360 | 0.700 | 0.876 | |
Felzartamab | 3.482 | 0.078 | 0.223 | ||
Followup:Felzartamab | 0.938 | 0.401 | 0.668 | ||
Atrophy classifier (ct>1Prob) | Followup | 0.568 | 0.572 | 0.780 | |
Felzartamab | 3.297 | 0.086 | 0.223 | ||
Followup:Felzartamab | 0.602 | 0.553 | 0.780 | ||
parenchyma-related | Kidney parenchymal (KT1) | Followup | 0.968 | 0.389 | 0.518 |
Felzartamab | 1.628 | 0.218 | 0.518 | ||
Followup:Felzartamab | 0.605 | 0.552 | 0.552 | ||
Kidney parenchymal - no solute carriers (KT2) | Followup | 0.859 | 0.432 | 0.518 | |
Felzartamab | 0.911 | 0.352 | 0.518 | ||
Followup:Felzartamab | 1.038 | 0.365 | 0.518 | ||
1.6 SUMMARIZE POST-HOC TESTS
Here we
present the effect of treatment (ΔΔ) for each score, along with the
effect of time (Δ) for placebo and felzartmab arms.
flextable_pairwiseTable 1. Effect of felzartamab on molecular scores | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Annotation | Score | Baseline - Week 24 | Week 24 - Week 52 | Baseline - Week 52 | |||||||||
Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | ||
cfDNA | Donor-derived cell-free DNA (dd-cfDNA, cp/mL) | 0.50 (76.12) | -55.00 (30) | -55.50 (53.06) | 2e-04 | -3.50 (69.12) | 14.50 (25.5) | 18.00 (47.31) | 0.027 | -3.00 (51) | -40.50 (45.75) | -37.50 (48.38) | 0.045 |
ABMR-related | ABMR classifier (ABMRProb) | 0.10 (0.37) | -0.57 (0.33) | -0.67 (0.35) | 1e-03 | -0.26 (0.44) | 0.54 (0.45) | 0.79 (0.45) | 1e-03 | -0.15 (0.43) | -0.03 (0.42) | 0.12 (0.43) | 0.709 |
Glomerulitis classifier (g>0Prob) | -0.02 (0.23) | -0.31 (0.35) | -0.29 (0.29) | 3e-03 | -0.19 (0.34) | 0.31 (0.41) | 0.50 (0.37) | 2e-04 | -0.22 (0.3) | 0.00 (0.38) | 0.21 (0.34) | 0.277 | |
Peritubular capillaritis classifier (ptc>0Prob) | -0.12 (0.25) | -0.55 (0.27) | -0.43 (0.26) | 3e-03 | -0.10 (0.37) | 0.37 (0.31) | 0.47 (0.34) | 3e-03 | -0.22 (0.3) | -0.18 (0.35) | 0.04 (0.33) | 0.977 | |
DSA-selective transcripts (DSAST) | -0.09 (0.22) | -0.37 (0.24) | -0.28 (0.23) | 0.026 | -0.16 (0.32) | 0.16 (0.33) | 0.32 (0.33) | 0.026 | -0.24 (0.33) | -0.20 (0.25) | 0.04 (0.29) | 0.988 | |
NK cell burden (NKB) | 0.02 (0.22) | -0.78 (0.49) | -0.80 (0.36) | 0.014 | -0.19 (0.4) | 0.68 (0.59) | 0.87 (0.49) | 0.014 | -0.17 (0.46) | -0.10 (0.62) | 0.07 (0.54) | 0.972 | |
TCMR-related | TCMR classifier (TCMRProb) | -0.01 (0.02) | 0.00 (0) | 0.01 (0.01) | 0.358 | 0.00 (0.02) | 0.00 (0) | 0.00 (0.01) | 0.856 | -0.01 (0.02) | 0.00 (0) | 0.01 (0.01) | 0.358 |
Tubulitis classifier (t>1Prob) | -0.01 (0.06) | -0.02 (0.04) | -0.01 (0.05) | 1.000 | 0.00 (0.04) | 0.01 (0.02) | 0.01 (0.03) | 1.000 | -0.01 (0.05) | -0.01 (0.03) | -0.01 (0.04) | 1.000 | |
Interstitial infiltrate classifier (i>1Prob) | -0.01 (0.04) | -0.02 (0.03) | -0.01 (0.04) | 0.964 | 0.01 (0.04) | 0.01 (0.04) | 0.00 (0.04) | 0.964 | 0.00 (0.05) | -0.01 (0.06) | -0.01 (0.05) | 0.964 | |
Cytotoxic T cell-associated (QCAT) | -0.20 (0.41) | -0.70 (0.46) | -0.50 (0.43) | 0.425 | 0.10 (0.33) | 0.54 (0.57) | 0.44 (0.45) | 0.425 | -0.10 (0.46) | -0.16 (0.41) | -0.06 (0.44) | 0.871 | |
T cell burden (TCB) | -0.44 (0.71) | -0.79 (0.56) | -0.36 (0.64) | 0.825 | 0.28 (0.57) | 0.34 (0.6) | 0.06 (0.58) | 0.905 | -0.16 (0.88) | -0.45 (0.43) | -0.29 (0.65) | 0.825 | |
macrophage-related | Alternatively activated macrophage (AMAT1) | 0.08 (0.2) | -0.27 (0.45) | -0.35 (0.33) | 0.128 | 0.09 (0.23) | 0.11 (0.47) | 0.02 (0.35) | 0.812 | 0.18 (0.19) | -0.16 (0.32) | -0.33 (0.25) | 0.128 |
Constitutive macrophage (QCMAT) | 0.00 (0.15) | -0.24 (0.27) | -0.24 (0.21) | 0.340 | 0.01 (0.12) | 0.09 (0.27) | 0.08 (0.2) | 0.685 | 0.01 (0.11) | -0.15 (0.24) | -0.16 (0.18) | 0.348 | |
injury-related | Injury-repair associated (IRRAT30) | 0.18 (0.52) | 0.05 (0.4) | -0.13 (0.46) | 0.353 | 0.05 (0.66) | -0.11 (0.72) | -0.16 (0.69) | 0.353 | 0.23 (0.6) | -0.06 (0.48) | -0.29 (0.54) | 0.175 |
Injury-repair induced, day 3 (IRITD3) | 0.03 (0.19) | 0.04 (0.12) | 0.01 (0.15) | 0.945 | 0.03 (0.18) | -0.18 (0.18) | -0.20 (0.18) | 0.060 | 0.06 (0.19) | -0.13 (0.14) | -0.19 (0.16) | 0.060 | |
Injury-repair induced, day 5 (IRITD5) | -0.01 (0.17) | 0.04 (0.14) | 0.05 (0.15) | 0.932 | 0.07 (0.2) | -0.17 (0.21) | -0.24 (0.21) | 0.093 | 0.06 (0.18) | -0.13 (0.17) | -0.18 (0.18) | 0.093 | |
Fibrosis classifier (ci>1Prob) | -0.02 (0.29) | 0.13 (0.32) | 0.15 (0.3) | 0.434 | 0.06 (0.25) | -0.08 (0.53) | -0.14 (0.39) | 0.844 | 0.04 (0.25) | 0.05 (0.37) | 0.01 (0.31) | 0.434 | |
Atrophy classifier (ct>1Prob) | 0.01 (0.33) | 0.12 (0.31) | 0.12 (0.32) | 0.587 | 0.09 (0.27) | -0.06 (0.44) | -0.15 (0.36) | 0.587 | 0.10 (0.32) | 0.07 (0.41) | -0.03 (0.36) | 0.587 | |
parenchyma-related | Kidney parenchymal (KT1) | -0.05 (0.22) | -0.15 (0.23) | -0.09 (0.22) | 0.731 | -0.01 (0.11) | 0.10 (0.22) | 0.11 (0.17) | 0.705 | -0.06 (0.25) | -0.04 (0.18) | 0.02 (0.22) | 0.705 |
Kidney parenchymal - no solute carriers (KT2) | -0.09 (0.39) | -0.16 (0.37) | -0.07 (0.38) | 0.577 | -0.09 (0.2) | 0.14 (0.38) | 0.23 (0.29) | 0.484 | -0.18 (0.42) | -0.02 (0.31) | 0.16 (0.37) | 0.577 | |
Grey shading denotes ANOVA interactive effect FDR < 0.05 | |||||||||||||
1.7 PLOT THE RESULTS
# EXTRACT LEGEND FOR PLOTS ####
panel_legend <- plots_artANOVA %>%
dplyr::filter(variable == "AMAT1") %>%
pull(plot_violin) %>%
ggpubr::get_legend() %>%
ggpubr::as_ggplot() +
theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))
# MOLECULAR ABMR PANELS ####
panel_violin_abmr <- plots_artANOVA %>%
dplyr::filter(category %in% c("ABMR")) %>%
pull(plot_violin) %>%
wrap_plots(nrow = 1, ncol = 5) &
theme(
axis.text = element_text(size = 10, colour = "black"),
legend.position = "none",
plot.background = element_rect(fill = "grey95", colour = " white")
)
panel_violin_abmr <- panel_violin_abmr %>%
ggarrange(
labels = "B",
font.label = list(size = 25, face = "bold"),
legend = "none"
) %>%
ggpubr::annotate_figure(
top = text_grob("Effect of felzartamab on molecular ABMR activity scores",
face = "bold.italic",
size = 25,
hjust = 1.27
)
)
# MOLECULAR TCMR PANELS ####
panel_violin_tcmr <- plots_artANOVA %>%
dplyr::filter(category %in% c("TCMR")) %>%
pull(plot_violin) %>%
wrap_plots(nrow = 1, ncol = 5) &
theme(
axis.text = element_text(size = 10, colour = "black"),
legend.position = "none",
)
panel_violin_tcmr <- panel_violin_tcmr %>%
ggarrange(
labels = "C",
font.label = list(size = 25, face = "bold"),
legend = "none"
) %>%
ggpubr::annotate_figure(
top = text_grob("Effect of felzartamab on molecular TCMR activity scores",
face = "bold.italic",
size = 25,
hjust = 1.27
)
)
# COMBINED VIOLIN PANELS ####
panels_violin <- ggarrange(
panel_violin_abmr,
panel_violin_tcmr,
nrow = 2,
heights = c(1, 1)
)
ggarrange(
panel_legend,
panels_violin,
nrow = 2,
heights = c(0.125, 1)
) 
2 The effect of treatment on slopes of injury scores
2.1 INTRODUCTION
While
ABMR-related scores changed monotonically over time, dipping at week 24
then rising again by week 52, initial observations of injury scores over
time alluded to a different relationship. Injury scores appeared to
gradual increase overtime in the placebo arm and a gradual decrease
overtime in the felzartamab arm. However, the use of artANOVA, a
non-parametric rank-based, test failed to detect such a relationship.
This led us to ask if estimating the slopes of injury curves over time
would be a more informative interrogation of these data. We applied a
linear mixed-model to answer this question.
References:
• Vonesh, E., et al. Mixed-effects models for slope-based endpoints in clinical trials of chronic kidney disease. Stat Med 38, 4218-4239 (2019).
2.2 DATA WRANGLING
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# CRAN packages
library(tidyverse)
library(lme4)
library(lmerTest)
library(performance)
library(ggeffects)
library(flextable)
library(ggpubr)
library(patchwork)
library(ggprism)
# Bioconductor libraries
library(Biobase) # BiocManager::install("Biobase")
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
source("C:/R/CD38-effect-of-treatment/code/functions/get_slope_function.r")
# load reference set
load("C:/R/CD38-effect-of-treatment/data/cd38.RData")
# load refression plots
load("C:/R/CD38-effect-of-treatment/results/plots_regression.RData")
# load violin plots
load("C:/R/CD38-effect-of-treatment/results/plots_artANOVA.RData")
# Seed for reproducibility
set.seed(42)
# DEFINE THE DATA ####
data <- set %>%
pData() %>%
dplyr::select(Center, Patient, Felzartamab, Group, IRRAT30, IRITD3, IRITD5) %>%
dplyr::filter(Patient %nin% c(15, 18)) %>%
mutate(
time = case_when(
Group == "Index" ~ 0,
Group == "FU1" ~ 0.46,
Group == "FU2" ~ 0.99
),
linetype = "solid" %>% factor()
)2.3 FIT THE MODELS
Now we
fit the linear mixed-models.
# FIT MODELS ####
fit_IRRAT30 <- lmer(IRRAT30 ~ Felzartamab * time + (time | Patient), data)
fit_IRITD3 <- lmer(IRITD3 ~ Felzartamab * time + (time | Patient), data)
fit_IRITD5 <- lmer(IRITD5 ~ Felzartamab * time + (time | Patient), data)
We check
the validatity of the model assumptions for IRRAT30.
# test model assumptions
performance::check_model(fit_IRRAT30)
performance::check_normality(fit_IRRAT30)## OK: residuals appear as normally distributed (p = 0.605).
We check
the validatity of the model assumptions for IRITD3.
# test model assumptions
performance::check_model(fit_IRITD3)
performance::check_normality(fit_IRITD3)## OK: residuals appear as normally distributed (p = 0.070).
We check
the validatity of the model assumptions for IRITD5.
# test model assumptions
performance::check_model(fit_IRITD5)
performance::check_normality(fit_IRITD5)## OK: residuals appear as normally distributed (p = 0.090).2.4 SUMMARIZE THE RESULTS
The next
step is to extract and summarize the model results.
# EXTRACT THE SLOPES FROM THE MODEL FITS ####
# IRRAT30
# Create list object into which to place model results
slopes_IRRAT30 <- list()
# Acute slopes
slopes_IRRAT30[[1]] <- get_slope(
.model_obj = fit_IRRAT30, .name = "Placebo",
.contrasts = c(0, 0, 1, 0) #
)
slopes_IRRAT30[[2]] <- get_slope(
.model_obj = fit_IRRAT30, .name = "Felzartamab",
.contrasts = c(0, 0, 1, 1)
)
slopes_IRRAT30[[3]] <- get_slope(
.model_obj = fit_IRRAT30, .name = "Intergroup Difference",
.contrasts = c(0, 0, 0, 1)
)
model_IRRAT30 <- slopes_IRRAT30 %>%
bind_rows() %>%
mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
dplyr::select(name, result, p_value, p_adjusted)
# IRITD3
# Create list object into which to place model results
slopes_IRITD3 <- list()
# Acute slopes
slopes_IRITD3[[1]] <- get_slope(
.model_obj = fit_IRITD3, .name = "Placebo",
.contrasts = c(0, 0, 1, 0) #
)
slopes_IRITD3[[2]] <- get_slope(
.model_obj = fit_IRITD3, .name = "Felzartamab",
.contrasts = c(0, 0, 1, 1)
)
slopes_IRITD3[[3]] <- get_slope(
.model_obj = fit_IRITD3, .name = "Intergroup Difference",
.contrasts = c(0, 0, 0, 1)
)
model_IRITD3 <- slopes_IRITD3 %>%
bind_rows() %>%
mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
dplyr::select(name, result, p_value, p_adjusted)
# IRITD5
# Create list object into which to place model results
slopes_IRITD5 <- list()
# Acute slopes
slopes_IRITD5[[1]] <- get_slope(
.model_obj = fit_IRITD5, .name = "Placebo",
.contrasts = c(0, 0, 1, 0) #
)
slopes_IRITD5[[2]] <- get_slope(
.model_obj = fit_IRITD5, .name = "Felzartamab",
.contrasts = c(0, 0, 1, 1)
)
slopes_IRITD5[[3]] <- get_slope(
.model_obj = fit_IRITD5, .name = "Intergroup Difference",
.contrasts = c(0, 0, 0, 1)
)
model_IRITD5 <- slopes_IRITD5 %>%
bind_rows() %>%
mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
dplyr::select(name, result, p_value, p_adjusted)
IRRAT30
size_font <- 10
model_IRRAT30 %>%
dplyr::select(-p_value) %>%
flextable::qflextable() %>%
flextable::set_table_properties(layout = "autofit") %>%
flextable::set_header_labels(
name = "Group", result = "IRRAT30 Slope (95%CI)",
"p_adjusted" = "FDR"
) %>%
flextable::bold(i = NULL, part = "header") %>%
flextable::fontsize(size = size_font, part = "all")Group | IRRAT30 Slope (95%CI) | FDR |
|---|---|---|
Placebo | 0.32 (-0.04 to 0.67) | 0.1133 |
Felzartamab | -0.29 (-0.64 to 0.07) | 0.1133 |
Intergroup Difference | -0.60 (-1.10 to -0.10) | 0.0546 |
IRITD3
model_IRITD3 %>%
dplyr::select(-p_value) %>%
flextable::qflextable() %>%
flextable::set_table_properties(layout = "autofit") %>%
flextable::set_header_labels(
name = "Group", result = "IRITD3 Slope (95%CI)",
"p_adjusted" = "FDR"
) %>%
flextable::bold(i = NULL, part = "header") %>%
flextable::fontsize(size = size_font, part = "all")Group | IRITD3 Slope (95%CI) | FDR |
|---|---|---|
Placebo | 0.06 (-0.06 to 0.18) | 0.3453 |
Felzartamab | -0.14 (-0.26 to -0.02) | 0.0408 |
Intergroup Difference | -0.20 (-0.37 to -0.02) | 0.0408 |
IRITD5
model_IRITD5 %>%
dplyr::select(-p_value) %>%
flextable::qflextable() %>%
flextable::set_table_properties(layout = "autofit") %>%
flextable::set_header_labels(
name = "Group", result = "IRITD5 Slope (95%CI)",
"p_adjusted" = "FDR"
) %>%
flextable::bold(i = NULL, part = "header") %>%
flextable::fontsize(size = size_font, part = "all")Group | IRITD5 Slope (95%CI) | FDR |
|---|---|---|
Placebo | 0.09 (-0.04 to 0.22) | 0.15720 |
Felzartamab | -0.13 (-0.26 to -0.01) | 0.05895 |
Intergroup Difference | -0.23 (-0.41 to -0.05) | 0.04200 |
2.5 PLOT THE RESULTS
# DEFINE LEGEND FOR VIOLIN PLOTS ####
panel_legend_violin <- plots_artANOVA %>%
dplyr::filter(variable == "AMAT1") %>%
pull(plot_violin) %>%
ggpubr::get_legend() %>%
ggpubr::as_ggplot() +
theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))
# DEFINE LEGEND FOR REGRESSION PLOTS ####
panel_legend_regression <- plots_regression %>%
dplyr::filter(variable == "IRRAT30") %>%
pull(plot) %>%
ggpubr::get_legend() %>%
ggpubr::as_ggplot()
# DEFIN JOINT LEGEND FOR REGRESSION AND VIOLIN PLOTS ####
panel_legend_all_injury <- panel_legend_violin + panel_legend_regression
# DEFINE INJURY VIOLIN PLOTS ####
plots_violin <- plots_artANOVA %>%
dplyr::filter(category %in% c("injury")) %>%
pull(plot_violin)
# MAKE PANEL OF SLOPE AND VIOLIN PLOTS ####
plot_panels_all_injury <- patchwork::wrap_plots(
plots_violin[[1]], plots_violin[[2]], plots_violin[[3]],
plots_regression$plot[[1]], plots_regression$plot[[2]], plots_regression$plot[[3]],
plots_regression$grob_table[[1]], plots_regression$grob_table[[2]], plots_regression$grob_table[[3]],
nrow = 3,
ncol = 3,
guides = "collect",
heights = c(1, 1, 0.5)
) +
patchwork::plot_annotation(tag_levels = list(c(LETTERS[1:6], rep("", 3)))) &
theme(
legend.position = "none",
# plot.title = element_blank(),
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
plot.tag = element_text(size = 20, face = "bold", vjust = 1)
)
ggarrange(
panel_legend_all_injury,
plot_panels_all_injury,
nrow = 2,
heights = c(0.15, 1.5)
)
3 The effect of treatment on gene expression
3.1 INTRODUCTION
In addition to
molecular scores, we had the opportunity to carry out genome-wide
assessement of the effect of felzartamab. We use the Linear Models for
Microarray Data (limma) pipeline for applying empirical Bayes mediated
t-tests on individual genes. Specific contrasts to assess the
interaction between treatment and time are specified to estimate the
effect of treatment on each gene.
References:
• Ritchie, M.E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).
3.2 DATA WRANGLING
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse)
library(flextable)
library(officer)
# Bioconductor libraries
library(Biobase)
library(limma)
library(genefilter)
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Load data
load("C:/R/CD38-effect-of-treatment/data/cd38.RData")
load("C:/R/CD38-effect-of-treatment/data/affymap219.RData")
# Seed for reproducibility
set.seed(42)
# DEFINE THE SET ####
set00 <- set[, set$Patient %nin% c(15, 18)]
• We
filter the gene expression data by to remove genes with minimal variance
across the population (IQR filtering).
# IQR FILTER THE DATA ####
f1 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1)
if (!exists("selected")) {
selected <- genefilter(set00, ff)
}
set01 <- set00[selected, ]
• We then
retain a single probeset for each gene based on whichever has the
highest mean expression in the population.
# KEEP UNIQUE GENES (keep probe with highest mean expression) ####
mean_exprs_by_probe <- set01 %>%
exprs() %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
mutate(mean_exprs = set01 %>%
exprs() %>% rowMeans(), .after = Symb)
genes <- mean_exprs_by_probe %>%
group_by(Symb) %>%
dplyr::slice_max(mean_exprs) %>%
dplyr::filter(Symb != "") %>%
distinct(Symb, .keep_all = TRUE) %>%
pull(AffyID)
set <- set01[featureNames(set01) %in% genes, ]3.3 REMOVE GENES THAT DIFFER AT BASELINE BETWEEN PLACEBO AND FELZARTAMAB ARMS
The
felzartamab trial was a phase II randomized control trial (RCT). While
RCTs are the ‘gold standard’ for accounting for confounders while
establishing causality (i.e., the effect of treatment), they are not
perfect, particularly when the number of study participants is low.
Thus, we pre-screened gene expression data to remove any genes that
differed between placebo and felzartamab arms at baseline. This was done
to reduce the potential for interpreting the effect of treatment on
genes with dissimilar baseline profiles.
# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set$Felzartamab_Followup %>% droplevels()
# BLOCK DESIGN week24 - baseline ####
design_block01 <- model.matrix(~ 0 + Felzartamab_Followup)
contrast_block_01 <- makeContrasts(
"x = (Felzartamab_FollowupBaseline_Felzartamab) -(Felzartamab_FollowupBaseline_Placebo)",
levels = design_block01
)
# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_block_1 <- limma::lmFit(set, design_block01)
cfit_block_1 <- limma::contrasts.fit(fit_block_1, contrast_block_01)
ebayes_block_1 <- limma::eBayes(cfit_block_1)
tab_block_1 <- limma::topTable(ebayes_block_1, adjust = "fdr", sort.by = "p", number = "all")
# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means_baseline <- fit_block_1 %>%
avearrays() %>%
data.frame() %>%
rownames_to_column("AffyID") %>%
tibble() %>%
mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
dplyr::select(AffyID, contains("Baseline"))
The top
20 genes that differ at baseline (sorted by p-value) between placebo and
felzartamab arms are summarized below.
(Note that all FDR are > 0.05. This is a consequence of differential expression across a large selection of genes (5,267 in this case) with small samples sizes (i.e., n = 10 by treatment group). Nonetheless, genes with p < 0.05 will be excluded from subsequent analyses).
# FORMAT TOPTABLES ####
baseline_deg <- tab_block_1 %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID") %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means_baseline, by = "AffyID") %>%
dplyr::select(
AffyID, Symb,
Baseline_Placebo, Baseline_Felzartamab,
t, logFC, P.Value, adj.P.Val,
) %>%
dplyr::rename(p = P.Value, FDR = adj.P.Val)
baseline_deg %>%
dplyr::slice_min(p, n = 20) %>%
mutate_at(vars(t, logFC, FDR), ~ formatC(.x, format = "f", digits = 2)) %>%
mutate_at(vars(p), ~ formatC(.x, format = "f", digits = 5)) %>%
flextable::flextable()AffyID | Symb | Baseline_Placebo | Baseline_Felzartamab | t | logFC | p | FDR |
|---|---|---|---|---|---|---|---|
11727331_at | SPAG4 | 90 | 50 | -4.23 | -0.85 | 0.00008 | 0.33 |
11763725_a_at | HNRNPU | 366 | 665 | 4.10 | 0.86 | 0.00013 | 0.33 |
11718286_a_at | UBL3 | 74 | 142 | 3.66 | 0.94 | 0.00053 | 0.45 |
11732590_x_at | TRIM5 | 42 | 63 | 3.65 | 0.59 | 0.00055 | 0.45 |
11725671_a_at | IPO7 | 42 | 66 | 3.54 | 0.64 | 0.00079 | 0.45 |
11741431_a_at | KCNJ16 | 203 | 311 | 3.49 | 0.61 | 0.00091 | 0.45 |
11726915_at | SLCO4C1 | 87 | 176 | 3.48 | 1.02 | 0.00096 | 0.45 |
11723109_a_at | GPBP1 | 153 | 263 | 3.46 | 0.78 | 0.00099 | 0.45 |
11739744_x_at | MXI1 | 123 | 191 | 3.45 | 0.64 | 0.00102 | 0.45 |
11720890_a_at | AGL | 39 | 63 | 3.40 | 0.70 | 0.00121 | 0.45 |
11746431_a_at | FBXO11 | 28 | 45 | 3.38 | 0.68 | 0.00127 | 0.45 |
11751841_a_at | TNFRSF17 | 56 | 20 | -3.36 | -1.51 | 0.00134 | 0.45 |
11744367_a_at | ZMYND11 | 30 | 49 | 3.33 | 0.68 | 0.00150 | 0.45 |
11748857_a_at | AK3 | 371 | 578 | 3.32 | 0.64 | 0.00155 | 0.45 |
11727307_x_at | ZNF320 | 46 | 66 | 3.31 | 0.51 | 0.00160 | 0.45 |
11747648_a_at | TCAIM | 96 | 152 | 3.29 | 0.66 | 0.00168 | 0.45 |
11728966_s_at | ZNF28 | 29 | 67 | 3.27 | 1.21 | 0.00181 | 0.45 |
11723387_a_at | PTAR1 | 149 | 261 | 3.25 | 0.81 | 0.00187 | 0.45 |
11739771_a_at | OXR1 | 31 | 51 | 3.24 | 0.73 | 0.00194 | 0.45 |
11730038_at | GNG13 | 44 | 28 | -3.24 | -0.63 | 0.00196 | 0.45 |
3.4 FIT THE MODELS
We now
carry out differential expression analysis to estimate the treatment
effect on individual genes.
# DEFINE GENES SIMILAR AT BASELINE ####
genes_baseline <- baseline_deg %>%
dplyr::filter(p > 0.05) %>%
pull(AffyID)
# WRANGLE THE MEAN EXPRESSION DATA ####
mean_exprs_by_probe <- set01 %>%
exprs() %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
mutate(
mean_exprs = set01 %>%
exprs() %>%
rowMeans(),
.after = Symb
)
genes <- mean_exprs_by_probe %>%
group_by(Symb) %>%
dplyr::slice_max(mean_exprs) %>%
dplyr::filter(Symb != "", AffyID %in% genes_baseline) %>%
distinct(Symb, .keep_all = TRUE) %>%
pull(AffyID)
set02 <- set01[featureNames(set01) %in% genes, ]
We need
to define the model structure for the contrasts we intend to make. Here
we will define 3 major contrast blocks; 1) effect of time in placebo
group, 2) effect of time in felzartamab group, and 3) the interactive
effect of time and treatment. Blocks 1 and 2 are assessed purely for
descriptive reasons as we are primarily concerned with block 3. We do
this for each time window (baseline to week 24, week 24 to week 52, and
baseline to week 52). We also include a molecular estimate for % cortex
as a covariate to adjust for difference in biopsy
composition.
# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set02$Felzartamab_Followup %>% droplevels()
cortex <- set02$Cortexprob
# DESIGN ####
design <- model.matrix(~ 0 + Felzartamab_Followup + cortex)
# CONTRAST DESIGN week24 - baseline ####
contrast_interaction_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
contrast_felzartamab_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
levels = design
)
contrast_placebo_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
# CONTRAST DESIGN week52 - week24 ####
contrast_interaction_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
levels = design
)
contrast_felzartamab_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2",
levels = design
)
contrast_placebo_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
levels = design
)
# CONTRAST DESIGN week52 - baseline####
contrast_interaction_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
contrast_felzartamab_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
levels = design
)
contrast_placebo_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
We can
check the specified contrast matrix for interaction terms. This is for
the baseline to week 24 window.
contrast_interaction_01 %>% as_tibble(rownames = "level")## # A tibble: 7 x 2
## level x = (Felzartamab_Fo~1
## <chr> <dbl>
## 1 Felzartamab_Fol~ 0.5
## 2 Felzartamab_Fol~ -0.5
## 3 Felzartamab_Fol~ 0
## 4 Felzartamab_Fol~ -0.5
## 5 Felzartamab_Fol~ 0.5
## 6 Felzartamab_Fol~ 0
## 7 cortex 0
## # i abbreviated name:
## # 1: `x = (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2`
and this
is for the baseline to week 52 window.
contrast_interaction_03 %>% as_tibble(rownames = "level")## # A tibble: 7 x 2
## level x = (Felzartamab_Fo~1
## <chr> <dbl>
## 1 Felzartamab_Fol~ 0.5
## 2 Felzartamab_Fol~ 0
## 3 Felzartamab_Fol~ -0.5
## 4 Felzartamab_Fol~ -0.5
## 5 Felzartamab_Fol~ 0
## 6 Felzartamab_Fol~ 0.5
## 7 cortex 0
## # i abbreviated name:
## # 1: `x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2`
We now
fit the models we’ve specified
# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_interaction_1 <- limma::lmFit(set02, design)
cfit_interaction_1 <- limma::contrasts.fit(fit_interaction_1, contrast_interaction_01)
ebayes_interaction_1 <- limma::eBayes(cfit_interaction_1)
tab_interaction_1 <- limma::topTable(ebayes_interaction_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_1$stdev.unscaled * sqrt(ebayes_interaction_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_1 <- limma::lmFit(set02, design)
cfit_felzartamab_1 <- limma::contrasts.fit(fit_felzartamab_1, contrast_felzartamab_01)
ebayes_felzartamab_1 <- limma::eBayes(cfit_felzartamab_1)
tab_felzartamab_1 <- limma::topTable(ebayes_felzartamab_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_1$stdev.unscaled * sqrt(ebayes_felzartamab_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_1 <- limma::lmFit(set02, design)
cfit_placebo_1 <- limma::contrasts.fit(fit_placebo_1, contrast_placebo_01)
ebayes_placebo_1 <- limma::eBayes(cfit_placebo_1)
tab_placebo_1 <- limma::topTable(ebayes_placebo_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_1$stdev.unscaled * sqrt(ebayes_placebo_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
# FIT BLOCK week52 - week24 LIMMA MODEL ####
fit_interaction_2 <- limma::lmFit(set02, design)
cfit_interaction_2 <- limma::contrasts.fit(fit_interaction_2, contrast_interaction_02)
ebayes_interaction_2 <- limma::eBayes(cfit_interaction_2)
tab_interaction_2 <- limma::topTable(ebayes_interaction_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_2$stdev.unscaled * sqrt(ebayes_interaction_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_2 <- limma::lmFit(set02, design)
cfit_felzartamab_2 <- limma::contrasts.fit(fit_felzartamab_2, contrast_felzartamab_02)
ebayes_felzartamab_2 <- limma::eBayes(cfit_felzartamab_2)
tab_felzartamab_2 <- limma::topTable(ebayes_felzartamab_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_2$stdev.unscaled * sqrt(ebayes_felzartamab_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_2 <- limma::lmFit(set02, design)
cfit_placebo_2 <- limma::contrasts.fit(fit_placebo_2, contrast_placebo_02)
ebayes_placebo_2 <- limma::eBayes(cfit_placebo_2)
tab_placebo_2 <- limma::topTable(ebayes_placebo_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_2$stdev.unscaled * sqrt(ebayes_placebo_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
# FIT BLOCK week52 - baseline LIMMA MODEL ####
fit_interaction_3 <- limma::lmFit(set02, design)
cfit_interaction_3 <- limma::contrasts.fit(fit_interaction_3, contrast_interaction_03)
ebayes_interaction_3 <- limma::eBayes(cfit_interaction_3)
tab_interaction_3 <- limma::topTable(ebayes_interaction_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_3$stdev.unscaled * sqrt(ebayes_interaction_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_3 <- limma::lmFit(set02, design)
cfit_felzartamab_3 <- limma::contrasts.fit(fit_felzartamab_3, contrast_felzartamab_03)
ebayes_felzartamab_3 <- limma::eBayes(cfit_felzartamab_3)
tab_felzartamab_3 <- limma::topTable(ebayes_felzartamab_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_3$stdev.unscaled * sqrt(ebayes_felzartamab_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_3 <- limma::lmFit(set02, design)
cfit_placebo_3 <- limma::contrasts.fit(fit_placebo_3, contrast_placebo_03)
ebayes_placebo_3 <- limma::eBayes(cfit_placebo_3)
tab_placebo_3 <- limma::topTable(ebayes_placebo_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_3$stdev.unscaled * sqrt(ebayes_placebo_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")3.5 SUMMARIZE THE RESULTS
The
results are now organized for tabulation
# MERGE BLOCKS ####
tab_block_1 <- tab_interaction_1 %>%
left_join(
tab_felzartamab_1 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_1 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
tab_block_2 <- tab_interaction_2 %>%
left_join(
tab_felzartamab_2 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_2 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
tab_block_3 <- tab_interaction_3 %>%
left_join(
tab_felzartamab_3 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_3 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means <- fit_interaction_1 %>%
avearrays() %>%
as_tibble(rownames = "AffyID") %>%
mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
dplyr::select(-contains("Week12"), -any_of(c("cortex")))
# FORMAT TOPTABLES ####
table_interaction_1 <- tab_block_1 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
table_interaction_2 <- tab_block_2 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
table_interaction_3 <- tab_block_3 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
deg <- tibble(
design = c(
"Baseline_vs_Week24",
"Week24_vs_Week52",
"Baseline_vs_Week52"
),
table = list(
table_interaction_1,
table_interaction_2,
table_interaction_3
)
)# FORMAT TABLES TO MAKE FLEXTABLES ####
table_deg <- deg %>%
expand_grid(direction = c("all", "increased", "decreased")) %>%
relocate(direction, .after = design) %>%
mutate(
gene_tables = pmap(
list(direction, table),
function(direction, table) {
if (direction == "increased") {
table <- table %>%
dplyr::filter(logFC > 0) %>%
arrange(p)
} else if (direction == "decreased") {
table <- table %>%
dplyr::filter(logFC < 0) %>%
arrange(p)
}
table %>%
dplyr::select(-t, -se, -FDR, -contains("AffyID"), -contains("MMDx")) %>%
dplyr::slice(1:20) %>%
mutate(
Gene = Gene %>% str_remove("///.*"),
plogFC = plogFC %>% round(2),
flogFC = flogFC %>% round(2),
logFC = logFC %>% round(2),
p = case_when(
p < 0.0001 ~ p %>% formatC(digits = 0, format = "e"),
TRUE ~ p %>% formatC(digits = 4, format = "f")
)
)
}
)
)
# GLOBAL PARAMETERS FOR FLEXTABLES ####
header1 <- c(
"Gene\nsymbol", "Gene", "PBT",
rep("Differential expression", 4),
rep("Mean expression by group", 6)
)
header2 <- c(
"Gene\nsymbol", "Gene", "PBT",
"\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
"\u394\u394\nlogFC", "\u394\u394\nP",
rep("Placebo", 3), rep("Felzartamab", 3)
)
header3 <- c(
"Gene\nsymbol", "Gene", "PBT",
"\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
"\u394\u394\nlogFC", "\u394\u394\nP",
"Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)", "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)"
)
# MAKE FORMATTED FLEXTABLES ####
flextables <- table_deg %>%
mutate(
flextable = pmap(
list(design, gene_tables),
function(design, gene_tables) {
if (design == "Baseline_vs_Week24") {
title <- paste("Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
} else if (design == "Week24_vs_Week52") {
title <- paste("Table i. Top 20 differentially expressed genes between week24 and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
} else if (design == "Baseline_vs_Week52") {
title <- paste("Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
}
gene_tables %>%
arrange(logFC) %>%
flextable::flextable() %>%
flextable::delete_part("header") %>%
flextable::add_header_row(top = TRUE, values = header3) %>%
flextable::add_header_row(top = TRUE, values = header2) %>%
flextable::add_header_row(top = TRUE, values = header1) %>%
flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
flextable::merge_v(part = "header") %>%
flextable::merge_h(part = "header") %>%
flextable::border_remove() %>%
flextable::border(part = "header", border = officer::fp_border()) %>%
flextable::border(part = "body", border = officer::fp_border()) %>%
flextable::align(align = "center") %>%
flextable::align(align = "center", part = "header") %>%
flextable::font(fontname = "Arial", part = "all") %>%
flextable::fontsize(size = 7, part = "all") %>%
flextable::fontsize(size = 7, part = "footer") %>%
flextable::fontsize(i = 1, size = 12, part = "header") %>%
flextable::bold(part = "header") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::padding(padding = 0, part = "all") %>%
flextable::autofit()
}
)
)
Genes
suppressed by felzartamab in week 24 compared to baseline
(Note
that all ΔΔFDR are > 0.05. This is a consequence of differential
expression across a large selection of genes (5,267 in this case) with
small samples sizes (i.e., n = 10 by treatment group). Thus, treatment
effects on individual genes should be interpreted with caution. Instead,
geneset enrichment analysis (GSEA) will be used to assess bulk
trends).
flextables %>%
dplyr::filter(design == "Baseline_vs_Week24", direction == "decreased") %>%
pull(flextable) %>%
pluck(1)Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value) | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Gene | Gene | PBT | Differential expression | Mean expression by group | ||||||||
Δ | Δ | ΔΔ | ΔΔ | Placebo | Felzartamab | |||||||
Baseline | Week24 | Week52 | Baseline | Week24 | Week52 | |||||||
CXCL11 | chemokine (C-X-C motif) ligand 11 | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.16 | -1.01 | -0.85 | 0.0176 | 494 | 395 | 372 | 617 | 153 | 317 |
CXCL9 | chemokine (C-X-C motif) ligand 9 | ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT | -0.14 | -0.94 | -0.79 | 0.0290 | 1,344 | 1,102 | 1,384 | 1,684 | 459 | 1,036 |
CXCL10 | chemokine (C-X-C motif) ligand 10 | ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT | -0.21 | -0.98 | -0.77 | 0.0228 | 844 | 632 | 720 | 1,063 | 274 | 615 |
IDO1 | indoleamine 2,3-dioxygenase 1 | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.09 | -0.77 | -0.68 | 0.0103 | 338 | 298 | 296 | 409 | 140 | 251 |
LYZ | lysozyme | IRITD5,QCMAT | -0.01 | -0.67 | -0.66 | 0.0386 | 1,676 | 1,658 | 1,884 | 2,244 | 884 | 1,184 |
FCGR3A | Fc fragment of IgG, low affinity IIIa, receptor (CD16a) | ABMR-RAT,GRIT2,IRRAT950,RAT,Rej-RAT | 0.06 | -0.58 | -0.64 | 0.0167 | 550 | 602 | 613 | 571 | 256 | 352 |
GZMB | granzyme B | ABMR-RAT,QCAT,RAT,Rej-RAT | -0.01 | -0.62 | -0.61 | 0.0019 | 149 | 147 | 130 | 143 | 61 | 114 |
GBP1 | guanylate binding protein 1, interferon-inducible | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.03 | -0.61 | -0.59 | 0.0096 | 678 | 655 | 643 | 909 | 389 | 556 |
GNLY | granulysin | ABMR-RAT,DSAST,QCAT,RAT,Rej-RAT | -0.18 | -0.75 | -0.57 | 0.0285 | 192 | 149 | 134 | 165 | 58 | 100 |
SH2D1B | SH2 domain containing 1B | ABMR-RAT,DSAST,NKB,RAT | 0.07 | -0.44 | -0.51 | 0.0010 | 38 | 42 | 34 | 40 | 22 | 33 |
FCGR1A | Fc fragment of IgG, high affinity Ia, receptor (CD64) | GRIT2,IRRAT950,RAT,TCMR-RAT | 0.09 | -0.41 | -0.50 | 0.0304 | 131 | 148 | 161 | 152 | 86 | 86 |
KLRD1 | killer cell lectin-like receptor subfamily D, member 1 | ABMR-RAT,RAT,Rej-RAT | -0.06 | -0.49 | -0.43 | 0.0295 | 110 | 101 | 96 | 126 | 63 | 102 |
EPB41L3 | erythrocyte membrane protein band 4.1-like 3 | 0.22 | -0.20 | -0.42 | 0.0388 | 209 | 286 | 258 | 339 | 258 | 283 | |
LOC339059 | uncharacterized LOC339059 | 0.26 | -0.15 | -0.41 | 0.0203 | 51 | 73 | 55 | 66 | 53 | 59 | |
APOL1 | apolipoprotein L1 | ABMR-RAT,GRIT3,RAT,Rej-RAT | 0.07 | -0.32 | -0.39 | 0.0205 | 626 | 687 | 744 | 754 | 483 | 636 |
LYPD5 | LY6/PLAUR domain containing 5 | ABMR-RAT,RAT,Rej-RAT | 0.04 | -0.33 | -0.37 | 0.0255 | 27 | 28 | 27 | 30 | 19 | 26 |
TRAV12-2 | T cell receptor alpha variable 12-2 | 0.25 | -0.10 | -0.35 | 0.0173 | 39 | 55 | 39 | 49 | 43 | 43 | |
IL18RAP | interleukin 18 receptor accessory protein | 0.02 | -0.32 | -0.34 | 0.0233 | 38 | 39 | 41 | 42 | 27 | 33 | |
PRF1 | perforin 1 (pore forming protein) | ABMR-RAT,QCAT,RAT,Rej-RAT | -0.02 | -0.35 | -0.33 | 0.0382 | 152 | 148 | 137 | 145 | 89 | 123 |
PMM2 | phosphomannomutase 2 | 0.20 | -0.07 | -0.27 | 0.0380 | 38 | 50 | 41 | 46 | 42 | 49 | |
Genes
suppressed by felzartamab in week 52 compared to baseline
(Note
that all ΔΔFDR are > 0.05. This is a consequence of differential
expression across a large selection of genes (5,267 in this case) with
small samples sizes (i.e., n = 10 by treatment group). Thus, treatment
effects on individual genes should be interpreted with caution. Instead,
geneset enrichment analysis (GSEA) will be used to assess bulk
trends).
flextables %>%
dplyr::filter(design == "Baseline_vs_Week52", direction == "decreased") %>%
pull(flextable) %>%
pluck(1)Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value) | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Gene | Gene | PBT | Differential expression | Mean expression by group | ||||||||
Δ | Δ | ΔΔ | ΔΔ | Placebo | Felzartamab | |||||||
Baseline | Week24 | Week52 | Baseline | Week24 | Week52 | |||||||
HNRNPK | heterogeneous nuclear ribonucleoprotein K | 0.28 | -0.38 | -0.66 | 0.0055 | 818 | 1,044 | 1,203 | 1,096 | 1,059 | 649 | |
COL1A1 | collagen, type I, alpha 1 | COLA1356,FICOL,IRITD5,LivGST_UP | 0.38 | -0.26 | -0.64 | 0.0075 | 310 | 333 | 526 | 329 | 320 | 230 |
NFIX | nuclear factor I/X (CCAAT-binding transcription factor) | CT1 | 0.42 | -0.15 | -0.57 | 0.0012 | 162 | 179 | 289 | 221 | 167 | 180 |
RHOA | ras homolog family member A | IRRAT950 | 0.28 | -0.29 | -0.57 | 0.0049 | 691 | 803 | 1,017 | 923 | 809 | 620 |
COL1A2 | collagen, type I, alpha 2 | FICOL,IRITD5,LivGST_UP | 0.26 | -0.31 | -0.57 | 0.0118 | 1,238 | 1,152 | 1,769 | 1,534 | 1,194 | 997 |
ITGB3 | integrin beta 3 | IRRAT30,IRRAT950 | 0.34 | -0.21 | -0.55 | 0.0100 | 105 | 137 | 169 | 110 | 103 | 82 |
LOC100129518 | uncharacterized LOC100129518 | GRIT3,IRRAT950 | 0.25 | -0.22 | -0.47 | 0.0050 | 429 | 546 | 605 | 527 | 461 | 390 |
ACTB | actin, beta | IRRAT950 | 0.15 | -0.31 | -0.47 | 0.0077 | 6,451 | 7,352 | 7,993 | 7,274 | 5,662 | 4,720 |
SPARC | secreted protein, acidic, cysteine-rich (osteonectin) | IRITD3 | 0.16 | -0.30 | -0.46 | 0.0100 | 2,411 | 2,532 | 3,011 | 2,343 | 2,094 | 1,550 |
DNAJA4 | DnaJ (Hsp40) homolog, subfamily A, member 4 | IRRAT950 | 0.27 | -0.14 | -0.42 | 0.0095 | 49 | 68 | 72 | 68 | 61 | 56 |
SLFN5 | schlafen family member 5 | 0.19 | -0.20 | -0.40 | 0.0050 | 39 | 44 | 51 | 49 | 38 | 37 | |
CYR61 | cysteine-rich, angiogenic inducer, 61 | IRITD3,IRRAT950 | 0.14 | -0.26 | -0.40 | 0.0074 | 271 | 268 | 330 | 338 | 234 | 236 |
SDCBP | syndecan binding protein | 0.24 | -0.14 | -0.38 | 0.0120 | 1,135 | 1,404 | 1,582 | 1,478 | 1,358 | 1,224 | |
SBSPON | somatomedin B and thrombospondin type 1 domain containing | 0.31 | -0.08 | -0.38 | 0.0123 | 217 | 250 | 332 | 188 | 249 | 169 | |
PAFAH1B2 | platelet-activating factor acetylhydrolase 1b, catalytic subunit 2 (30kDa) | 0.23 | -0.13 | -0.36 | 0.0099 | 309 | 372 | 423 | 386 | 396 | 323 | |
TMEM2 | transmembrane protein 2 | 0.10 | -0.26 | -0.36 | 0.0118 | 270 | 307 | 311 | 318 | 259 | 222 | |
CALM2 | calmodulin 2 (phosphorylase kinase, delta) | CT1,IRITD3 | 0.12 | -0.25 | -0.36 | 0.0138 | 1,304 | 1,427 | 1,538 | 1,623 | 1,342 | 1,155 |
WNK1 | WNK lysine deficient protein kinase 1 | 0.26 | -0.08 | -0.34 | 0.0104 | 703 | 884 | 1,012 | 894 | 996 | 799 | |
CEBPB | CCAAT/enhancer binding protein (C/EBP), beta | IRITD3 | 0.11 | -0.22 | -0.33 | 0.0123 | 419 | 452 | 488 | 417 | 333 | 307 |
STAT6 | signal transducer and activator of transcription 6, interleukin-4 induced | 0.09 | -0.20 | -0.30 | 0.0144 | 671 | 710 | 765 | 746 | 616 | 563 | |
4 Geneset enrichment analysis to identify pathways affected by felzartamab treatment
4.1 INTRODUCTION
We applied
geneset enrichment analysis (GSEA) to identify pathways affected by
felzartamab treatment. Multiple libraries were screened to encompass a
broad spectrum of biological pathways, including Gene Ontology (GO),
Disease Ontology Semantic and Enrichment analysis (DOSE), Kyoto
Encyclopedia of Genes and Genomes (KEGG), Reactome, Molecular Signatures
Database (MsigDB), and WikiPathways (wiki) libraries. An additional
de-novo injury geneset was also assesed. This injury geneset was based
on single-nuclei transcriptomics assessed from healthy and injury kidney
allografts. Geneset enrichment was chosen specifically because it is a
rank-based test and does not require cutoffs based on p-values or fold
change. Nonetheless, here we decided apriori that we would pre-screen
genes based on their p-values from differential expression analyses
described in Section 3. The intent for this pre-screening was to limit
interpretation from GSEA to genes which had the greatest measured
difference between placebo and felzartamab arms over time. Note that
there are many ways for carrying out pathway and functional enrichment
analyses. Also note that geneset libraries are frequently updated and
you may find slightly different results based on the package version
being used. For example, here we encountered such a difference in the
DOSE library, where versions later than 3.30.1 identified different
pathways than version 3.30.1. Thus, we have included a tarball of the
package version (v3.30.1) used for the original analyses in this paper
for convenient installation and replication of the presented
findings.
References:
• Yu, G., Wang, L.G., Han, Y. & He, Q.Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284-287 (2012).
• Wu, T., et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2, 100141 (2021).
• Hinze, C., et al. Single-cell transcriptomics reveals common epithelial response patterns in human acute kidney injury. Genome Med 14, 103 (2022).
4.2 DATA WRANGLING
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse)
library(flextable)
library(officer)
library(msigdbr)
# Bioconductor libraries
# DOSE v 3.30.1 required for reproducibility
# install.packages("C:/R/CD38-effect-of-treatment/data/DOSE_3.30.1.tar.gz", repos = NULL, type = "source")
library(DOSE)
library(ReactomePA)
library(Biobase)
library(clusterProfiler)
library(org.Hs.eg.db)
library(pRoloc)
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# load DE data
load("C:/R/CD38-effect-of-treatment/results/deg.RData")
# Seed for reproducibility
set.seed(42)4.3 GENESET ENRICHMENT ANALYSES
Geneset
enrichment analyses using the ClusterProfiler package have simlar, but
not identical implementation. This section will carry out GSEA on all
described libraries, after which the results will be collated and
summarized. For all GSEA, genes are ranked by
t-value.
GO
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_go <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::gseGO(
gene = genes_gsea, ont = "BP", OrgDb = org.Hs.eg.db,
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# SIMPLIFY GO TERMS ####
set.seed(42)
gsea_go_simplified <- gsea_go %>%
mutate(gsea_simplified = map(gsea, clusterProfiler::simplify))
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_go_simplified_formatted <- gsea_go_simplified %>%
mutate(
gsea_tables = map(
gsea_simplified,
function(gsea_simplified) {
gsea_simplified %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_go_tables <- gsea_go_simplified_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_go <- gsea_go_tables %>%
mutate(db = "GO", .before = 1)
names(gsea_go$gsea_flextables) <- gsea_go$design
DOSE
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = clusterProfiler::bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_dose <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
DOSE::gseDO(
gene = genes_gsea,
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_dose_formatted <- gsea_dose %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"))
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_dose_tables <- gsea_dose_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_dose <- gsea_dose_tables %>%
mutate(db = "DOSE", .before = 1)
names(gsea_dose$gsea_flextables) <- gsea_dose$design
KEGG
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_kegg <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::gseKEGG(
gene = genes_gsea, organism = "hsa",
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_kegg_formatted <- gsea_kegg %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_kegg_tables <- gsea_kegg_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_kegg <- gsea_kegg_tables %>%
mutate(db = "kegg", .before = 1)
names(gsea_kegg$gsea_flextables) <- gsea_kegg$design
MSigDB
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# DEFINE MSigDB PATHWAYS ####
gene_sets <- msigdbr(category = "H")
map <- gene_sets[, c("gs_name", "entrez_gene")]
map$entrez_gene <- as.character(map$entrez_gene)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_msigdb <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::GSEA(
gene = genes_gsea, TERM2GENE = map,
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_msigdb_formatted <- gsea_msigdb %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_msigdb_tables <- gsea_msigdb_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
# slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_msigdb <- gsea_msigdb_tables %>%
mutate(db = "msigdb", .before = 1)
names(gsea_msigdb$gsea_flextables) <- gsea_msigdb$design
REACTOME
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_reactome <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
ReactomePA::gsePathway(
gene = genes_gsea, organism = "human",
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_reactome_formatted <- gsea_reactome %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_reactome_tables <- gsea_reactome_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_reactome <- gsea_reactome_tables %>%
mutate(db = "Reactome", .before = 1)
names(gsea_reactome$gsea_flextables) <- gsea_reactome$design
WikiPathways
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_wiki <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::gseWP(
gene = genes_gsea, organism = "Homo sapiens",
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_wiki_formatted <- gsea_wiki %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_wiki_tables <- gsea_wiki_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
# slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_wiki <- gsea_wiki_tables %>%
mutate(db = "wiki", .before = 1)
names(gsea_wiki$gsea_flextables) <- gsea_wiki$design
AKI
geneset
# load gene lists
load("C:/R/CD38-effect-of-treatment/data/hinze2022_injury_markers.RData")
# DEFINE INJURY PATHWAYS ####
map_aki <- injury_markers %>% dplyr::select(gs_name, entrez_gene)
map_aki_joined <- injury_markers %>%
dplyr::select(gs_name, entrez_gene) %>%
mutate(gs_name = "AKI")
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_aki <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::GSEA(
gene = genes_gsea, TERM2GENE = map_aki_joined,
minGSSize = 5, maxGSSize = Inf,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_aki_formatted <- gsea_aki %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = "Cellular response to acute kidney injury"
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_aki_tables <- gsea_aki_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_aki <- gsea_aki_tables %>%
mutate(db = "Hinze", .before = 1)
names(gsea_aki$gsea_flextables) <- gsea_aki$design4.4 SUMMARIZE THE RESULTS
GSEA
results from each library can now be combined and summarized. Note that
pathway annotations need to be interpreted in light of their
developement in combination with how GSEA determines their enrichment,
and the relevance of enriched pathways must be evaluated in the context
of transplantation. For our interpretation, we avail to the pathway
annotation as well as core enrichment genes contributing to enrichment
via GSEA. We have simplified the interpretation of enriched pathways
based on the overlap in biological function of their core enrichement
genes.
# JOIN THE GSEA RESULTS ####
gsea <- purrr::reduce(
list(
gsea_go,
gsea_msigdb,
gsea_wiki,
gsea_dose,
gsea_kegg,
gsea_reactome,
gsea_aki
), bind_rows
) %>%
mutate(keep = map_dbl(gsea_tables, nrow)) %>%
dplyr::filter(keep > 0) %>%
dplyr::select(db, design, genes_gsea, gsea_tables) %>%
mutate(gsea_tables = map(gsea_tables, mutate, Description = Description %>% as.character()))
# DEFINE INTERPREATION BASED ON PATHWAYS ####
pathway_interpretation <- gsea %>%
dplyr::select(design, gsea_tables) %>%
unnest(everything()) %>%
mutate(
interpretation = case_when(
ID %in% c(
"GO:0006952",
"GO:0009605",
"GO:0009607",
"GO:0043207",
"GO:0044419",
"GO:0051707",
"GO:0006955",
"GO:0002376",
"DOID:2914",
"DOID:77",
"DOID:1579",
"R-HSA-168256"
) & design == "Baseline_vs_Week24" ~ "immune response",
ID %in% c(
"R-HSA-168256",
"R-HSA-168249",
"R-HSA-1643685",
"R-HSA-5663205"
) & design == "Baseline_vs_Week52" ~ "immune-related response to injury",
ID %in% c("AKI", "R-HSA-109582") ~ "response to injury",
ID %in% c(
"R-HSA-9006934",
"R-HSA-597592",
"R-HSA-5653656",
"R-HSA-392499"
) ~ "response to injury",
ID %in% c(
"R-HSA-2262752",
"R-HSA-8953897"
) ~ "response to injury",
)
) %>%
dplyr::select(design, ID, Description, interpretation)
# WRANGLE THE GSEA TABLES ####
gsea_tables <- gsea %>%
mutate(
gsea_tables = pmap(
list(design, gsea_tables),
function(design, gsea_tables) {
gsea_tables %>%
mutate(
design = design,
Description = Description %>% as.character()
) %>%
left_join(pathway_interpretation, by = c("ID", "Description", "design")) %>%
dplyr::select(ID, Description, interpretation, setSize, NES, p.adjust, core_enrichment) %>%
distinct(ID, setSize, interpretation, .keep_all = TRUE)
}
)
) %>%
unnest(gsea_tables) %>%
arrange(p.adjust) %>%
dplyr::filter(p.adjust < 0.001) %>%
nest(.by = design, gsea_table = c(-design, -genes_gsea)) %>%
left_join(gsea %>% dplyr::select(design, genes_gsea) %>% distinct(design, .keep_all = TRUE),
by = "design"
)
# MAKE FLEXTABLES ####
gsea_flextables <- gsea_tables %>%
mutate(
flextable = pmap(
list(design, gsea_table),
function(design, gsea_table) {
cellWidths <- c(2, 3, 7, 4, 1.5, 1.5, 1.5, 13) # for individual tables up or down
if (design == "Baseline_vs_Week24") {
title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week24 (by FDR)", sep = "")
} else if (design == "Week24_vs_Week52") {
title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between week24 and week52 (by FDR)", sep = "")
} else if (design == "Baseline_vs_Week52") {
title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week52 (by FDR)", sep = "")
}
gsea_table %>%
mutate(
library = db,
pathway = Description,
"core enrichment genes" = core_enrichment %>% str_replace_all("/", ", "),
"n genes" = setSize,
NES = NES %>% round(2),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(library, ID, pathway, interpretation, "n genes", NES, FDR, "core enrichment genes") %>%
flextable::flextable() %>%
flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
flextable::merge_v(part = "header") %>%
flextable::merge_h(part = "header") %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::padding(padding = 0) %>%
flextable::padding(j = "core enrichment genes", padding.left = 5) %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core enrichment genes", align = "left", part = "body") %>%
flextable::font(fontname = "Arial", part = "all") %>%
flextable::fontsize(size = 10, part = "body") %>%
flextable::fontsize(size = 12, part = "header") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::width(width = cellWidths, unit = "cm")
}
)
)
Summary
of all enriched pathways at week 24 (compared to
baseline)
# PRINT FLETABLES TO POWERPOINT ####
gsea_flextables %>%
dplyr::filter(design == "Baseline_vs_Week24") %>%
pull(flextable) %>%
pluck(1) Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week24 (by FDR) | |||||||
|---|---|---|---|---|---|---|---|
library | ID | pathway | interpretation | n genes | NES | FDR | core enrichment genes |
GO | GO:0006952 | defense response | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0009605 | response to external stimulus | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0009607 | response to biotic stimulus | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0043207 | response to external biotic stimulus | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0044419 | biological process involved in interspecies interaction between organisms | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0051707 | response to other organism | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0006955 | immune response | immune response | 23 | -2.93 | 5e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, TRAV12-2, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0002376 | immune system process | immune response | 24 | -2.75 | 0.0002 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, TRAV12-2, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
Reactome | R-HSA-168256 | Immune System | immune response | 14 | -2.51 | 0.0005 | PIK3AP1, PTPRC, TRIM4, KLRF1, LYZ, FCGR1A, KLRD1, GNLY, IL18RAP, CXCL10, FCGR3A, GBP1, SH2D1B |
DOSE | DOID:2914 | immune system disease | immune response | 10 | -2.56 | 0.0007 | PRF1, FCGR1A, CXCL9, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB |
DOSE | DOID:77 | gastrointestinal system disease | immune response | 15 | -2.59 | 0.0007 | PTPRC, EPB41L3, LYZ, PRF1, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB |
DOSE | DOID:1579 | respiratory system disease | immune response | 10 | -2.56 | 0.0007 | PRF1, FCGR1A, CXCL9, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB |
Summary
of all enriched pathways at week 52 (compared to
baseline)
gsea_flextables %>%
dplyr::filter(design == "Baseline_vs_Week52") %>%
pull(flextable) %>%
pluck(1) Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week52 (by FDR) | |||||||
|---|---|---|---|---|---|---|---|
library | ID | pathway | interpretation | n genes | NES | FDR | core enrichment genes |
Reactome | R-HSA-168249 | Innate Immune System | immune-related response to injury | 20 | -3.18 | 6e-06 | LRRFIP1, C1QB, JUN, ATP6V1A, NCF2, CFI, EP300, FCGR1A, CAPZA1, CAPZA2, ADAM10, HSPA1A, STAT6, CALM2, SDCBP, PAFAH1B2, ACTB, RHOA |
Reactome | R-HSA-168256 | Immune System | immune-related response to injury | 35 | -3.08 | 2e-05 | LRRFIP1, C1QB, CANX, JUN, ATP6V1A, RNF213, NCF2, PIK3AP1, VIM, CFI, EP300, FCGR1A, CAPZA1, CAPZA2, ADAM10, HSPA1A, IFNGR1, STAT6, CALM2, SDCBP, COL1A2, PAFAH1B2, ACTB, COL1A1, RHOA |
Reactome | R-HSA-1643685 | Disease | immune-related response to injury | 32 | -2.92 | 2e-05 | LRRFIP1, PPFIBP1, CANX, JUN, RNF213, CD163, PIK3AP1, TPM3, SIN3A, EP300, FCGR1A, RAN, ADAM10, HSPA1A, IFNGR1, CALM2, CEBPB, SBSPON, COL1A2, ITGB3, ACTB, COL1A1, HNRNPK |
Hinze | AKI | Cellular response to acute kidney injury | response to injury | 39 | -2.75 | 4e-05 | SPTBN1, NRP1, KLF6, HLA-DRB1, DOCK11, DUSP1, CHSY1, LDHA, CDC42SE2, TPM4, PTBP3, ARHGAP29, SYNE2, SERP1, CORO1C, LRRFIP1, PPFIBP1, CANX, JUN, SCOC, RNF213, SPP1, PIK3AP1, VIM, SVIL, TPM3, RAN, ANXA5, ADAM10, HSPA1A, CEBPB, SPARC, ACTB, HNRNPK, SLFN5, RHOA |
Reactome | R-HSA-9006934 | Signaling by Receptor Tyrosine Kinases | response to injury | 19 | -2.59 | 0.0002 | ATP6V1A, SPP1, NCF2, SIN3A, EP300, ADAM10, STAT6, CALM2, COL1A2, SPARC, ITGB3, ACTB, COL1A1, RHOA |
Reactome | R-HSA-597592 | Post-translational protein modification | response to injury | 17 | -2.68 | 0.0002 | CANX, SPP1, PRMT3, PPP6R3, SIN3A, EP300, CAPZA1, CAPZA2, ADAM10, CALM2, SBSPON, ACTB, HNRNPK, RHOA |
Reactome | R-HSA-109582 | Hemostasis | response to injury | 17 | -2.54 | 0.0005 | SIN3A, CAPZA1, CAPZA2, ANXA5, CALM2, COL1A2, SPARC, ITGB3, ACTB, COL1A1, RHOA |
Reactome | R-HSA-392499 | Metabolism of proteins | response to injury | 19 | -2.42 | 0.0005 | PRMT3, PPP6R3, SIN3A, EP300, CAPZA1, CAPZA2, ADAM10, CALM2, SBSPON, ACTB, HNRNPK, RHOA |
Reactome | R-HSA-5653656 | Vesicle-mediated transport | response to injury | 15 | -2.47 | 0.0008 | SCOC, AP4B1, CD163, PPP6R3, CAPZA1, CAPZA2, CALM2, COL1A2, SPARC, PAFAH1B2, ACTB, COL1A1 |