######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("ggplot2")
#install.packages("ggpubr")


#引用包
library(limma)
library(reshape2)
library(ggplot2)
library(ggpubr)

expFile="symbol.txt"           #表达数据文件
ClusterFile="cluster.txt"      #分型结果文件
geneFile="gene.txt"            #免疫检查点的基因列表文件
setwd("C:\\biowolf\\Necroptosis\\38.clusterCheckpoint")     #设置工作目录

#读取基因表达文件,并对数据进行处理
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
	
#读取基因列表文件，提取免疫检查点相关基因的表达量
gene=read.table(geneFile, header=F, sep="\t", check.names=F)
sameGene=intersect(row.names(data), as.vector(gene[,1]))
data=t(data[sameGene,])
data=log2(data+1)

#删除正常样品
group=sapply(strsplit(row.names(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[group==0,]
row.names(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", row.names(data))
data=avereps(data)
	
#合并数据
Cluster=read.table(ClusterFile, sep="\t", header=T, check.names=F, row.names=1)
sameSample=intersect(row.names(data), row.names(Cluster))
rt1=cbind(data[sameSample,], Cluster[sameSample,])
rt1=rt1[,c(sameGene, "Cluster")]

#提取显著差异的基因
sigGene=c()
for(i in colnames(rt1)[1:(ncol(rt1)-1)]){
	if(sd(rt1[,i])<0.001){next}
	if(length(levels(factor(rt1[,"Cluster"])))>2){
		test=kruskal.test(rt1[,i] ~ rt1[,"Cluster"])
	}else{
		test=wilcox.test(rt1[,i] ~ rt1[,"Cluster"])
	}
	pvalue=test$p.value
	if(pvalue<0.05){
		sigGene=c(sigGene, i)
	}
}
sigGene=c(sigGene, "Cluster")
rt1=rt1[,sigGene]

#把数据转换成ggplot2输入文件
rt1=melt(rt1,id.vars=c("Cluster"))
colnames(rt1)=c("Cluster","Gene","Expression")

#绘制箱线图
bioCol=c("#0066FF","#FF9900","#00DB00","#FF0000","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length(unique(rt1$Cluster))]
boxplot=ggboxplot(rt1, x="Gene", y="Expression", fill="Cluster",
				  xlab="",
				  ylab="Gene expression",
				  legend.title="Cluster",
				  width=0.8,
				  palette = bioCol)+
				  rotate_x_text(50)+
	stat_compare_means(aes(group=Cluster),
		symnum.args=list(cutpoints=c(0, 0.001, 0.01, 0.05, 1), symbols=c("***", "**", "*", "")), label="p.signif")
	
#输出图片
pdf(file="checkpoint.diff.pdf", width=7.6, height=5)
print(boxplot)
dev.off()


######Video source: https://ke.biowolf.cn
######生信自学网: https://www.biowolf.cn/
######微信公众号：biowolf_cn
######合作邮箱：biowolf@foxmail.com
######答疑微信: 18520221056

