--- a
+++ b/Signatures/IPS.R
@@ -0,0 +1,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)))
\ No newline at end of file