######Video source: https://ke.biowolf.cn
######??????ѧ??: https://www.biowolf.cn/
######΢?Ź??ںţ?biowolf_cn
######???????䣺biowolf@foxmail.com
######????΢??: 18520221056

#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma")

#install.packages("pheatmap")


#???ð?
library(limma)
library(pheatmap)

clusterFile="cluster.txt"      #???͵Ľ????ļ?
immFile="infiltration_estimation_for_tcga.csv"     #????ϸ???????????ļ?
setwd("C:\\Users\\LZD\\Desktop\\THCA\\PANoptosis\\37.clusterImmCor")      #???ù???Ŀ¼

#??ȡ???ͽ????ļ?
cluster=read.table(clusterFile, header=T, sep="\t", check.names=F, row.names=1)

#??ȡ????ϸ???????ļ?
immune=read.csv(immFile, header=T, sep=",", check.names=F, row.names=1)
immune=as.matrix(immune)
rownames(immune)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*)", "\\1\\-\\2\\-\\3", rownames(immune))
immune=avereps(immune)

#???˷??ͺ?????ϸ???ϲ?
sameSample=intersect(row.names(cluster), row.names(immune))
cluster=cluster[sameSample, "Cluster", drop=F]
immune=immune[sameSample, , drop=F]
data=cbind(cluster, immune)

#???͵?????ϸ??????????
outTab=data.frame()
sigCell=c("Cluster")
for(i in colnames(data)[2:ncol(data)]){
	if(sd(data[,i])<0.001){next}
	if(length(levels(factor(data[,"Cluster"])))>2){
		test=kruskal.test(data[,i] ~ data[,"Cluster"])
	}else{
		test=wilcox.test(data[,i] ~ data[,"Cluster"])
	}
	pvalue=test$p.value
	if(pvalue<0.05){
		outTab=rbind(outTab,cbind(immune=i, pvalue))
		sigCell=c(sigCell, i)
	}
}
write.table(file="immuneCor.txt", outTab, sep="\t", quote=F, row.names=F)

#??ͼ????
data=data[,sigCell]
#???͵?ע??
data=data[order(data[,"Cluster"]),]
annCol=data[,1,drop=F]
annCol[,"Cluster"]=factor(annCol[,"Cluster"], unique(annCol[,"Cluster"]))
data=t(data[,(2:ncol(data))])
#???????͵?ע??
annRow=sapply(strsplit(rownames(data),"_"), '[', 2)
annRow=as.data.frame(annRow)
row.names(annRow)=row.names(data)
colnames(annRow)=c("Methods")
annRow[,"Methods"]=factor(annRow[,"Methods"], unique(annRow[,"Methods"]))
gapCol=as.vector(cumsum(table(annCol[,"Cluster"])))
gapRow=as.vector(cumsum(table(annRow[,"Methods"])))

#??????ͼע?͵???ɫ
bioCol=c("#0066FF","#FF9900","#00DB00","#FF0000","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length(unique(annCol[,"Cluster"]))]
Cluster=bioCol
names(Cluster)=levels(factor(annCol[,"Cluster"]))
ann_colors=list(Cluster=Cluster)

#??ͼ???ӻ?
pdf("immHeatmap.pdf", width=9, height=6)
pheatmap(data,
         annotation=annCol,
         annotation_row=annRow,
         annotation_colors = ann_colors,
         color = colorRampPalette(c(rep("#0099CC",5), "white", rep("#CC0033",5)))(100),
         cluster_cols =F,
         cluster_rows =F,
         gaps_row=gapRow,
         gaps_col=gapCol,
         scale="row",
         show_colnames=F,
         show_rownames=T,
         fontsize=6,
         fontsize_row=5,
         fontsize_col=6)
dev.off()


######Video source: https://ke.biowolf.cn
######??????ѧ??: https://www.biowolf.cn/
######΢?Ź??ںţ?biowolf_cn
######???????䣺biowolf@foxmail.com
######????΢??: 18520221056

