Immunologic response to liver allograft injury

V.I Log-transformation

df1 <- read.csv("402df.csv")

par(mfrow = c(3, 3))
hist(df1$CYTO1, main="CYTO1")
hist(df1$CYTO2, main="CYTO2")
hist(df1$CYTO3, main="CYTO3")
hist(df1$CYTO10, main="CYTO10")
hist(df1$CYTO11, main="CYTO11")
hist(df1$CYTO12, main="CYTO12")
hist(df1$CYTO18, main="CYTO18")
hist(df1$CYTO19, main="CYTO19")
hist(df1$CYTO20, main="CYTO20")

As we can see from the histograms of some of the histograms, there is a strong right skewness for all of the cytokine levels. Therefore, there exist a motivation to log-transform the data.

cols <- c("CYTO1","CYTO2", "CYTO3","CYTO4", "CYTO5","CYTO6", "CYTO7","CYTO8",
          "CYTO9","CYTO10", "CYTO11","CYTO12", "CYTO13","CYTO14", "CYTO15","CYTO16",
          "CYTO17","CYTO18", "CYTO19","CYTO20", "CYTO21","CYTO22", "CYTO23")
df1[cols] <- log(df1[cols])

par(mfrow = c(3, 3))
hist(df1$CYTO1, main="CYTO1")
hist(df1$CYTO2, main="CYTO2")
hist(df1$CYTO3, main="CYTO3")
hist(df1$CYTO10, main="CYTO10")
hist(df1$CYTO11, main="CYTO11")
hist(df1$CYTO12, main="CYTO12")
hist(df1$CYTO18, main="CYTO18")
hist(df1$CYTO19, main="CYTO19")
hist(df1$CYTO20, main="CYTO20")

Now the selected histograms have improved significantly. We are going to proceed out analysis with the transformed dataset.

V.II Wilcoxon test

#Prepare data
df <- read.csv("402df.csv")
df_visual <- df %>% select(CYTO1:CYTO23) 
df <- df %>% mutate(IRI=as.factor(IRI))
df[,paste0("CYTO",c(1:23))] <- log(df[,paste0("CYTO",c(1:23))])
df_dif <- df
df_dif[which(df_dif$Stype=="PO"), paste0("CYTO",c(1:23))] <-
  df_dif[which(df_dif$Stype=="LF"), paste0("CYTO",c(1:23))] - 
  df_dif[which(df_dif$Stype=="PO"), paste0("CYTO",c(1:23))]
df_dif[which(df_dif$Stype=="LF"), paste0("CYTO",c(1:23))] <-
  df_dif[which(df_dif$Stype=="W1"), paste0("CYTO",c(1:23))] - 
  df_dif[which(df_dif$Stype=="LF"), paste0("CYTO",c(1:23))]
df_dif[which(df_dif$Stype=="W1"), paste0("CYTO",c(1:23))] <-
  df_dif[which(df_dif$Stype=="W1"), paste0("CYTO",c(1:23))] - 
  df_dif[which(df_dif$Stype=="PO"), paste0("CYTO",c(1:23))]
df_dif$Stype[which(df_dif$Stype=="W1")] <- "W1-PO"
df_dif$Stype[which(df_dif$Stype=="LF")] <- "W1-LF"
df_dif$Stype[which(df_dif$Stype=="PO")] <- "LF-PO"

Firstly, we evaluated each cytokine for its correlation with IRI outcome individually based on the Wilcoxon nonparametric rank-sum test.

# Wilcoxon nonparametric rank-sum test
df_test <- rbind(df, df_dif)
stat.test <- NULL
for(k in paste0("CYTO",c(1:23))){
  for(i in unique(df_test$Stype)){
  temp <- df_test %>%
    filter(Stype==i) %>%
    rename("CYTO" = k) %>%
    wilcox_test(CYTO ~ IRI) %>%
    add_significance() %>%
    mutate(Stype=i,var=k)
  stat.test <- rbind(stat.test, temp)
  }
}

stat.test[,-1] %>% 
  select(var, Stype, group1:p.signif) %>%
  filter(p.signif!="ns") %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")

As we can tell from the summary table, CYTO16 value in W1 stage and the change of CYTO16 value between W1 and LF stage are significant different between IRI+ and IRI- groups. In addition, CYTO22 value in LF stage are significant different between IRI+ and IRI- groups.

VI.I Main effects

# Data preparation
df_1 <- df1 %>% 
  group_by(PID) %>% 
  filter(Stype == 'PO') %>%
  select(c(4, 8:30))

df_2 <- df1 %>% 
  group_by(PID) %>% 
  filter(Stype == 'LF') %>%
  select(c(4, 8:30))

df_3 <- df1 %>% 
  group_by(PID) %>% 
  filter(Stype == 'W1') %>%
  select(c(4, 8:30))

# Get their differences for each patient
df_1num <- df_1[, 4:25]
df_2num <- df_2[, 4:25]
df_3num <- df_3[, 4:25]
diff1 <- df_2num - df_1num
diff2 <- df_3num - df_2num
diff3 <- df_3num - df_1num

df_diff1 <- df_1
df_diff1[, 4:25] <- diff1

df_diff2 <- df_1
df_diff2[, 4:25] <- diff2

df_diff3 <- df_1
df_diff3[, 4:25] <- diff3

Then we can fit the model with main effects only:

# Before
mod1_before <- glm(IRI ~ .-PID, family = binomial, df_1)
  # No significant cytokine

# During
mod1_during <- glm(IRI ~ .-PID, family = binomial, df_2)
  # CTYO01, 02, 014 were significant

# After
mod1_after <- glm(IRI ~ .-PID, family = binomial, df_3)
  # CTYO22

# Difference 1
mod1_diff1 <- glm(IRI ~ .-PID, family = binomial, df_diff1)
  # CTYO04, 09, 014 were significant

# Difference 2
mod1_diff2 <- glm(IRI ~ .-PID, family = binomial, df_diff2)
    # CTYO09

# Difference 3
mod1_diff3 <- glm(IRI ~ .-PID, family = binomial, df_diff3)
  # No significant cytokine

After fitting the logistic models with only the main effect, CTYO01, 02, 09, 014, 022 have shown significance.

VI.II With interactions

Then we include all the two-way interactions between cytokines to test whether cytokines influence IRI together. In order to arrive at a more parsimonious model, we simplified the model by a elastic net selection algorithm. Elastic net linear regression applied the penalties that from both the lasso and ridge techniques to regularize regression models. It combines both the lasso and ridge regression methods and improve the shortcomings of these two methods. The technique is appropriate where the dimensional data is greater than the number of samples used. Note that we build the model based on each Stype (PO, LF, W1, LF-PO, W1-LF, W1-PO) independently.

# log the data
df_test2 <- model.matrix(~.^2-1,df_test[,paste0("CYTO",c(1:23))])
glm_results <- NULL
for(i in unique(df_test$Stype)){
  # logit regression
  set.seed(1996)
  
  y <- df_test %>% filter(Stype==i) %>% select(IRI) %>% as.matrix
  x <- df_test2[which(df_test$Stype==i),] %>% as.matrix
   
  cv.lasso <- cv.glmnet(x, y, family = "binomial")
  #plot(cv.lasso)
  model <- glmnet(x, y, family = "binomial",
                  lambda = cv.lasso$lambda.1se)
  ind <- abs(coef(model)[,ncol(coef(model))])>0.000000000
  vars.selected <- rownames(coef(model))[ind]
  vars.selected <- vars.selected[-which(vars.selected == "(Intercept)")]
  
  df_final <- data.frame(as.numeric(y), x[,vars.selected]) 
  names(df_final) <- c("y", vars.selected)
  glm(y~., data=df_final, 
      family=binomial) %>% summary %>%
    .$coefficients %>% as.data.frame -> temp
  temp$var <- row.names(temp)
  glm_results <- rbind(glm_results, data.frame(Stype=i, temp))
}
colnames(glm_results)[c(3,5)] <- c("sd","p-value")
glm_results <- glm_results[which(grepl("(Intercept)", 
                                       row.names(glm_results))==F),]
glm_results[,2:5] <- glm_results[,2:5] %>% round(3)
row.names(glm_results)<-NULL

Following is the summary table for models,

glm_results %>%
  select("var","Stype","Estimate","sd","z.value","p-value") %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")

As we can tell from the model summary, CYTO23 is proved to be significant different between IRI+ and IRI- groups at PO stage (p = .035). CYTO22 appears to be different between IRI+ and IRI- groups at LF stage (p = .065). We also found that there exists significant interaction effect between the CYTO5 (change between LF and PO stage) and CYTO8 (change between LF and PO stage) between IRI+ and IRI- groups (p = .022). Interaction effect between the CYTO12 (change between W1 and LF stage) and CYTO14 (change between W1 and LF stage) between IRI+ and IRI- groups (p = .038) is also significant.

VI.III NEC/INF

For the next section, we analyzied the correlation between re-transplant cytokines and after-transplant by spearman correlation

#data processing
df <- read.csv("402df.csv")
df[,"INF"] <- as.factor(df[,"INF"])
df[,"NEC"] <- as.factor(df[,"NEC"])
pre_index <- seq(1,300,3)
after_index <- seq(3,300,3)
df_test <- NULL
df_test$INF_after <- df[after_index,"INF"]
df_test$NEC_after <- df[after_index,"NEC"]
df_test <- cbind(df_test, df[pre_index,paste0("CYTO",c(1:23))])
df_test[,paste0("CYTO",c(1:23))] <- log(df_test[,paste0("CYTO",c(1:23))])

We apply the spearman correlation (Koelman) with INF and NEC as factors, we see no strong correlations between INF and NEC after transplant with Cytokines before transplant nor are they significant.

correlation::correlation(df_test[1:2], df_test[3:25], method = "spearman",
                         include_factors = TRUE)

Visualize the correlation to see the weak correlations between cytokines and INF/NEC score.
Note here NEC and INC are treated as numeric variables.

df_test[,"INF_after"] <- as.numeric(df_test[,"INF_after"])
df_test[,"NEC_after"] <- as.numeric(df_test[,"NEC_after"])
corrplot(cor(df_test,method = "spearman"),
  method = "color")