Master Regulators Driving response in ICB-based immunotherapy
Signatures
source("code/MR.R")
load("Data/LungMR_GRN.RData")
Tx = read.table(file = "Signature/Gene_signature-CTL.txt", row.names = NULL, header = FALSE)
CTL = as.character(Tx[, 1])
Tx = read.table(file = "Signature/Gene_signature-dysfunction.txt", row.names = 1, header = TRUE)
dysfunction = list(pos = rownames(filter(Tx, Correlation == "Positive")),
neg = rownames(filter(Tx, Correlation == "Negtive")) )Master regulators
T1 = mrMulti(GEX, GRN, gs = list(CTL = CTL), nsample = 50)
write.table(T1, file = "out/LungGRN_CTL.txt", row.names = FALSE, col.names = TRUE, sep = "\t")
T2 = mrMulti(GEX, GRN, gs = dysfunction, nsample = 50)
write.table(T2, file = "out/LungGRN_dysfunction.txt", row.names = FALSE, col.names = TRUE, sep = "\t")T1 = read.table(file = "out/LungGRN_CTL.txt", row.names = NULL, header = TRUE, sep = "\t")
T1$signature = "CTL"
T2 = read.table(file = "out/LungGRN_dysfunction.txt", row.names = NULL, header = TRUE, sep = "\t")
T2$signature = "dysfunction"
mergedT = rbind(T1, T2)
uTable = filter(mergedT, p.value < 0.01, NES > 0)
dTable = filter(mergedT, p.value < 0.01, NES < 0)Activated master regulators
Comparison
x = table(uTable$Regulon, uTable$signature)
y = as.data.frame.matrix(x)
z = melt(updateTable(table(y[,1], y[,2]), space = 0:7, N00 = TRUE)) %>% filter(Var1 + Var2 > 0)
p = ggplot(z, aes(Var1, Var2)) + theme_bw() +
labs(x = "# of data sets\n(T-cell abundance)",
y = "# of data sets\n(T-cell dysfunction)", title = "") +
geom_tile(aes(fill = value)) +
geom_text(aes(label = round(value, 1))) +
scale_fill_gradient(low = "white", high = "red")
print(p)
write.table(y, file = "out/Immunotherapy_MR_Activated.txt",
row.names = TRUE, col.names = TRUE, sep = "\t")
filter(y, CTL >=6, dysfunction >=6) CTL dysfunction
ARHGAP25 7 7
BTK 6 6
IKZF1 6 6
IL16 7 7
ITGAL 6 6
STK10 6 6Pathway enrichment
source("/sibcb2/bioinformatics/Script/PathwayEnrichment.R")
load("/sibcb2/bioinformatics/KnowledgeBase/MSigDB/MSigDB.RData")
PG = as.character(read.table(file = "Data/VIPER_Regulators.txt", row.names = NULL, col.names=TRUE,sep = "\t")[,1])
pwList = MSigDB$MSigDB_c5_bp_v5
TopCTL = rownames(filter(as.data.frame.matrix(x), CTL >= 4))
TopDysfunction = rownames(filter(as.data.frame.matrix(x), dysfunction >= 4))
EnrichCTL = pHyperList(TopCTL, pwList, PG) %>% arrange(pvalue)
EnrichDysfunction = pHyperList(TopDysfunction, pwList, PG) %>% arrange(pvalue)
write.table(EnrichCTL, file = "out/Immunotherapy_CTL_pathway_Activated.txt",
row.names = TRUE, col.names = TRUE, sep = "\t")
EnrichCTL[1:10, 1:5] sg sr n pvalue qvalue
IMMUNE_SYSTEM_PROCESS 119 197 32 1.351041e-22 1.114609e-19
IMMUNE_RESPONSE 119 132 21 4.224908e-15 1.742775e-12
T_CELL_ACTIVATION 119 31 9 1.843440e-10 5.069459e-08
DEFENSE_RESPONSE 119 147 17 3.457886e-10 7.131889e-08
T_CELL_DIFFERENTIATION 119 8 5 1.379221e-09 2.275715e-07
LYMPHOCYTE_ACTIVATION 119 43 9 6.534982e-09 8.985601e-07
LEUKOCYTE_ACTIVATION 119 48 9 2.052431e-08 2.418936e-06
CELL_ACTIVATION 119 50 9 3.118522e-08 3.026801e-06
HEMOPOIESIS 119 50 9 3.118522e-08 3.026801e-06
HEMOPOIETIC_OR_LYMPHOID_ORGAN_DEVELOPMENT 119 52 9 4.646523e-08 3.833382e-06write.table(EnrichDysfunction, file = "out/Immunotherapy_Dysfunction_pathway_Activated.txt",
row.names = TRUE, col.names = TRUE, sep = "\t")
EnrichDysfunction[1:10,1:5] sg sr n pvalue qvalue
IMMUNE_SYSTEM_PROCESS 90 197 21 4.505862e-14 3.717336e-11
T_CELL_ACTIVATION 90 31 7 1.038583e-08 4.284156e-06
LEUKOCYTE_ACTIVATION 90 48 8 2.465437e-08 6.779951e-06
T_CELL_DIFFERENTIATION 90 8 4 3.511420e-08 7.242303e-06
CELL_ACTIVATION 90 50 8 3.594750e-08 5.392125e-06
HEMOPOIESIS 90 50 8 3.594750e-08 5.392125e-06
HEMOPOIETIC_OR_LYMPHOID_ORGAN_DEVELOPMENT 90 52 8 5.152011e-08 6.072012e-06
IMMUNE_SYSTEM_DEVELOPMENT 90 53 8 6.130925e-08 6.322516e-06
LYMPHOCYTE_ACTIVATION 90 43 7 1.650571e-07 1.513023e-05
LEUKOCYTE_DIFFERENTIATION 90 27 5 2.101731e-06 1.733928e-04shownPathway = names(pwList)
EnrichMAP = data.frame(pathway = shownPathway,
CTL = EnrichCTL[shownPathway,"qvalue"],
dysfunction = EnrichDysfunction[shownPathway, "qvalue"])
p = ggplot(EnrichMAP, aes(-log10(CTL), -log10(dysfunction))) + theme_bw() +
geom_point(color = "black", alpha = 0.5, shape = 1) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_vline(xintercept = -log10(0.05), linetype = "longdash", color = "black", size = 0.5) +
geom_hline(yintercept = -log10(0.05), linetype = "longdash", color = "black", size = 0.5)
print(p)
Repressed master regulators
Comparison
x = table(dTable$Regulon, dTable$signature)
y = as.data.frame.matrix(x)
z = melt(updateTable(table(y[,1], y[,2]), space = 0:7, N00 = TRUE)) %>% filter(Var1 + Var2 > 0)
p = ggplot(z, aes(Var1, Var2)) + theme_bw() +
labs(x = "# of data sets\n(T-cell abundance)",
y = "# of data sets\n(T-cell dysfunction)", title = "") +
geom_tile(aes(fill = value)) +
geom_text(aes(label = round(value, 1))) +
scale_fill_gradient(low = "white", high = "red")
print(p)
write.table(y, file = "out/Immunotherapy_MR_Repressed.txt",
row.names = TRUE, col.names = TRUE, sep = "\t")
filter(y, dysfunction >=6) CTL dysfunction
COPS5 0 6
DDX1 0 6
GGCT 2 6
MSH2 0 6
TAF2 0 6
TFB2M 2 7
ZNHIT3 0 6Pathway enrichment
pwList = MSigDB$MSigDB_h_all_v5
TopDysfunction = rownames(filter(as.data.frame.matrix(x), dysfunction >= 4))
EnrichDysfunction = pHyperList(TopDysfunction, pwList, PG) %>% arrange(pvalue)
write.table(EnrichDysfunction, file = "out/Immunotherapy_Dysfunction_pathway_Repressed.txt",
row.names = TRUE, col.names = TRUE, sep = "\t")
filter(EnrichDysfunction, qvalue < 0.05) sg sr n pvalue qvalue
HALLMARK_MYC_TARGETS_V1 73 59 11 2.645560e-12 1.322780e-10
HALLMARK_E2F_TARGETS 73 92 7 1.252881e-05 3.132203e-04
HALLMARK_MYC_TARGETS_V2 73 16 3 3.174254e-05 5.290424e-04
HALLMARK_OXIDATIVE_PHOSPHORYLATION 73 21 2 1.914255e-03 2.392819e-02
HALLMARK_FATTY_ACID_METABOLISM 73 23 2 2.505319e-03 2.505319e-02
HALLMARK_MTORC1_SIGNALING 73 57 3 4.746314e-03 3.955262e-02
overlapVector
HALLMARK_MYC_TARGETS_V1 CBX3 COPS5 DEK HDAC2 HSPE1 ILF2 ORC2 PRDX3 SRPK1 YWHAQ PTGES3
HALLMARK_E2F_TARGETS ASF1A DEK IPO7 MSH2 ORC2 RAD21 RBBP7
HALLMARK_MYC_TARGETS_V2 CBX3 HSPE1 TFB2M
HALLMARK_OXIDATIVE_PHOSPHORYLATION FH PRDX3
HALLMARK_FATTY_ACID_METABOLISM FH HSP90AA1
HALLMARK_MTORC1_SIGNALING COPS5 EEF1E1 HSPE1gs = gsub("HALLMARK_", "", rownames(EnrichDysfunction))[1:6]
gs = factor(gs, levels = gs)
mTable = data.frame(gs = gs, LogFDR = -log10(EnrichDysfunction$qvalue)[1:6])
p <- ggplot(mTable, aes(x = gs, y = LogFDR)) +
geom_bar(stat="identity", fill = "black") + theme_bw() +
theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
ylab("-log10(FDR)") + xlab("") + theme(axis.text.x=element_text(angle = 60, hjust = 1, size = 9)) +
geom_hline(yintercept = -log10(0.05), color = "red", linetype = "longdash")
print(p)