-
Notifications
You must be signed in to change notification settings - Fork 26
/
IPS.R
143 lines (130 loc) · 8.35 KB
/
IPS.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
####################################################
##
## This R-script can be used to calculate Immunophenoscore (IPS) and generate Immunophenogram from "EXPR.txt" and "IPS_genes.txt"
## (C) ICBI, Medical University of Innsbruck, Biocenter, Division of Bioinformatics
## Version 1.0 08.07.2016
## Needs packages ggplot2,grid,gridExtra
##
####################################################
library(ggplot2)
library(grid)
library(gridExtra)
## calculate Immunophenoscore
ipsmap<- function (x) {
if (x<=0) {
ips<-0
} else {
if (x>=3) {
ips<-10
} else {
ips<-round(x*10/3, digits=0)
}
}
return(ips)
}
## Assign colors
my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 1000)
mapcolors<-function (x) {
za<-NULL
if (x>=3) {
za=1000
} else {
if (x<=-3) {
za=1
} else {
za=round(166.5*x+500.5,digits=0)
}
}
return(my_palette[za])
}
my_palette2 <- colorRampPalette(c("black", "white"))(n = 1000)
mapbw<-function (x) {
za2<-NULL
if (x>=2) {
za2=1000
} else {
if (x<=-2) {
za2=1
} else {
za2=round(249.75*x+500.5,digits=0)
}
}
return(my_palette2[za2])
}
## Read expression data from tab-delimited text file, with official human gene symbols (HGNC) in the first columns
## and expression values (i.e. log2(TPM+1)) for each sample in the other columns
gene_expression<-read.table("EXPR.txt",row.names=1,header=TRUE, sep="\t", dec = ".",check.names=FALSE)
sample_names<-names(gene_expression)
## Read IPS genes and corresponding weights from tab-delimited text file "IPS_genes.txt"
# For different
IPSG<-read.table("IPS_genes.txt",header=TRUE, sep="\t", dec = ".",check.names=FALSE)
unique_ips_genes<-as.vector(unique(IPSG$NAME))
IPS<-NULL
MHC<-NULL
CP<-NULL
EC<-NULL
SC<-NULL
AZ<-NULL
# Gene names in expression file
GVEC<-row.names(gene_expression)
# Genes names in IPS genes file
VEC<-as.vector(IPSG$GENE)
# Match IPS genes with genes in expression file
ind<-which(is.na(match(VEC,GVEC)))
# List genes missing or differently named
MISSING_GENES<-VEC[ind]
dat<-IPSG[ind,]
if (length(MISSING_GENES)>0) {
cat("differently named or missing genes: ",MISSING_GENES,"\n")
}
for (x in 1:length(ind)) {
print(IPSG[ind,])
}
for (i in 1:length(sample_names)) {
GE<-gene_expression[[i]]
mGE<-mean(GE)
sGE<-sd(GE)
Z1<-(gene_expression[as.vector(IPSG$GENE),i]-mGE)/sGE
W1<-IPSG$WEIGHT
WEIGHT<-NULL
MIG<-NULL
k<-1
for (gen in unique_ips_genes) {
MIG[k]<- mean(Z1[which (as.vector(IPSG$NAME)==gen)],na.rm=TRUE)
WEIGHT[k]<- mean(W1[which (as.vector(IPSG$NAME)==gen)])
k<-k+1
}
WG<-MIG*WEIGHT
MHC[i]<-mean(WG[1:10])
CP[i]<-mean(WG[11:20])
EC[i]<-mean(WG[21:24])
SC[i]<-mean(WG[25:26])
AZ[i]<-sum(MHC[i],CP[i],EC[i],SC[i])
IPS[i]<-ipsmap(AZ[i])
## Plot Immunophenogram
data_a <- data.frame (start = c(0,2.5,5,7.5,10,15,seq(20,39),0,10,20,30), end = c(2.5,5,7.5,10,15,seq(20,40),10,20,30,40), y1=c(rep(2.6,26),rep(0.4,4)),y2=c(rep(5.6,26),rep(2.2,4)),z=c(MIG[c(21:26,11:20,1:10)],EC[i],SC[i],CP[i],MHC[i]),vcol=c(unlist(lapply(MIG[c(21:26,11:20,1:10)],mapcolors)), unlist(lapply(c(EC[i],SC[i],CP[i],MHC[i]),mapbw))), label = c(unique_ips_genes[c(21:26,11:20,1:10)],"EC","SC","CP","MHC"))
data_a$label <- factor(data_a$label, levels=unique(data_a$label))
plot_a1<-ggplot() + geom_rect(data=data_a, mapping=aes(xmin=start, xmax=end, ymin=y1, ymax=y2, fill=label), size=0.5,color="black", alpha=1) + coord_polar() + scale_y_continuous(limits = c(0, 6)) + scale_fill_manual(values =as.vector(data_a$vcol),guide=FALSE) + theme_bw() + theme(panel.margin = unit(0, 'mm'), panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.border = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "white"), axis.text=element_blank(), axis.ticks= element_blank()) + geom_text(aes(x=5, y=1.3, label="EC"), size=4) + geom_text(aes(x=15, y=1.3, label="SC"), size=4) + geom_text(aes(x=25, y=1.3, label="CP"), size=4) + geom_text(aes(x=35, y=1.3, label="MHC"), size=4)
plot_a2<-plot_a1+geom_text(aes(x=1.25, y=4.1, label="+ Act CD4"), angle=78.75, size=4)+geom_text(aes(x=3.75, y=4.1, label="+ Act CD8"),angle=56.25, size=4)+geom_text(aes(x=6.25, y=4.1, label="+ Tem CD4"), angle=33.75,size=4)+geom_text(aes(x=8.75, y=4.1, label="+ Tem CD8"), angle=11.25,size=4)+geom_text(aes(x=17.5, y=4.1, label="- MDSC"), angle=-67.5,size=4)+geom_text(aes(x=12.5, y=4.1, label="- Treg"), angle=-22.5,size=4)
plot_a3<-plot_a2+geom_text(aes(x=20.5, y=4.1, label="PD-1 -"), angle=85.5, size=4)+geom_text(aes(x=21.5, y=4.1, label="CTLA4 -"), angle=76.5, size=4)+geom_text(aes(x=22.5, y=4.1, label="LAG3 -"), angle=67.5, size=4)+geom_text(aes(x=23.5, y=4.1, label="TIGIT -"), angle=58.5, size=4)+geom_text(aes(x=24.5, y=4.1, label="TIM3 -"), angle=49.5, size=4)+geom_text(aes(x=25.5, y=4.1, label="PD-L1 -"), angle=40.5, size=4)+geom_text(aes(x=26.5, y=4.1, label="PD-L2 -"), angle=31.5, size=4)+geom_text(aes(x=27.5, y=4.1, label="CD27 +"), angle=22.5, size=4)+geom_text(aes(x=28.5, y=4.1, label="ICOS +"), angle=13.5, size=4)+geom_text(aes(x=29.5, y=4.1, label="IDO1 -"), angle=4.5, size=4)
plot_a4<-plot_a3+geom_text(aes(x=30.5, y=4.1, label="B2M +"), angle=-4.5, size=4)+geom_text(aes(x=31.5, y=4.1, label="TAP1 +"), angle=-13.5, size=4)+geom_text(aes(x=32.5, y=4.1, label="TAP2 +"), angle=-22.5, size=4)+geom_text(aes(x=33.5, y=4.1, label="HLA-A +"), angle=-31.5, size=4)+geom_text(aes(x=34.5, y=4.1, label="HLA-B +"), angle=-40.5, size=4)+geom_text(aes(x=35.5, y=4.1, label="HLA-C +"), angle=-49.5, size=4)+geom_text(aes(x=36.5, y=4.1, label="HLA-DPA1 +"), angle=-58.5, size=4)+geom_text(aes(x=37.5, y=4.1, label="HLA-DPB1 +"), angle=-67.5, size=4)+geom_text(aes(x=38.5, y=4.1, label="HLA-E +"), angle=-76.5, size=4)+geom_text(aes(x=39.5, y=4.1, label="HLA-F +"), angle=-85.5, size=4)
plot_a5<-plot_a4+geom_text(aes(x=0, y=6, label=paste("Immunophenoscore: ",IPS[i],sep="")), angle=0,size=6,vjust=-0.5)+ theme(axis.title=element_blank())
plot_a <-plot_a5 + theme(plot.margin=unit(c(0,0,0,0),"mm")) + geom_text(vjust=1.15,hjust=0,aes(x=25.5, y=6,label="\n\n\n\n MHC: Antigen Processing EC: Effector Cells\n CP: Checkpoints | Immunomodulators SC: Suppressor Cells\n\n", hjust = 0), size=4)
## Legend sample-wise (averaged) z-scores
data_b <- data.frame (start = rep(0,23), end = rep(0.7,23), y1=seq(0,22,by=1), y2=seq(1,23,by=1),z=seq(-3,3,by=6/22),vcol=c(unlist(lapply(seq(-3,3,by=6/22),mapcolors))), label = LETTERS[1:23])
data_b_ticks <- data.frame(x = rep(1.2, 7), value = seq(-3,3, by=1), y = seq(0,6, by=1)*(22/6) +0.5)
legendtheme <- theme(plot.margin = unit(c(2,0,2,0),"inch"), panel.margin = unit(0,"null"), panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.border = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "white"), axis.text=element_blank(), axis.ticks= element_blank(), axis.title.x=element_blank())
plot_b<-ggplot(hjust=0) + geom_rect(data=data_b, mapping=aes(xmin=start, xmax=end, ymin=y1, ymax=y2, fill=label), size=0.5,color="black", alpha=1) + scale_x_continuous(limits = c(0, 1.5),expand = c(0,0)) + scale_fill_manual(values =as.vector(data_b$vcol),guide=FALSE) + geom_text(data=data_b_ticks, aes(x=x, y=y, label=value),hjust="inward", size=4) +theme_bw() + legendtheme + ylab("Sample-wise (averaged) z-score")
## Legend weighted z-scores
data_c <- data.frame (start = rep(0,23), end = rep(0.7,23), y1=seq(0,22,by=1), y2=seq(1,23,by=1),z=seq(-2,2,by=4/22),vcol=c(unlist(lapply(seq(-2,2,by=4/22),mapbw))), label = LETTERS[1:23])
data_c_ticks <- data.frame(x = rep(1.2, 5), value = seq(-2,2, by=1), y = seq(0,4, by=1)*(22/4) +0.5)
plot_c<-ggplot() + geom_rect(data=data_c, mapping=aes(xmin=start, xmax=end, ymin=y1, ymax=y2, fill=label), size=0.5,color="black", alpha=1) + scale_x_continuous(limits = c(0, 1.5),expand = c(0,0)) + scale_fill_manual(values =as.vector(data_c$vcol),guide=FALSE) + geom_text(data=data_c_ticks, aes(x=x, y=y, label=value),hjust="inward", size=4) +theme_bw() + legendtheme + ylab("Weighted z-score")
## Save plot to file (1 pdf file for each sample)
file_name<-paste("IPS_",sample_names[i],".pdf",sep="")
pdf(file_name, width=10, height=8)
grid.arrange(plot_a,plot_b,plot_c, ncol=3, widths=c(0.8,0.1,0.1))
dev.off()
}
DF<-data.frame(SAMPLE=sample_names,MHC=MHC,EC=EC,SC=SC,CP=CP,AZ=AZ,IPS=IPS)
write.table(DF,file="IPS.txt",row.names=FALSE, quote=FALSE,sep="\t")
barplot(DF[[7]],DF[[1]],col=c(rep("yellow",7),rep("green",4)))