install.packages("pheatmap")install.packages("limma ")library(pheatmap)library("limma")

setwd("F:\ 96glycolysis\09.glyGeneExp\m6a\8")          inputFile="sampleExp.txt"                                       #输入文件pvalFilter=0.05                                                  #p值临界值logFCfilter=0                                                    #logFC临界值conNum=41                                                         #normal组样品数目treatNum=476                                                      #tumor组样品数目           #读取输入文件outTab=data.frame()grade=c(rep(1,conNum),rep(2,treatNum))rt=read.table(inputFile,sep="t",header=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)data=data[rowMeans(data)>0,]           #差异分析for(i in row.names(data)){        geneName=unlist(strsplit(i,"\|",))[1]        geneName=gsub("\/", "_", geneName)        rt=rbind(expression=data[i,],grade=grade)        rt=as.matrix(t(rt))        wilcoxTest<-wilcox.test(expression ~ grade, data=rt)        conGeneMeans=mean(data[i,1:conNum])        treatGeneMeans=mean(data[i,(conNum+1):ncol(data)])        logFC=log2(treatGeneMeans)-log2(conGeneMeans)        pvalue=wilcoxTest$p.value        conMed=median(data[i,1:conNum])        treatMed=median(data[i,(conNum+1):ncol(data)])        diffMed=treatMed-conMed        if( ((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){               outTab=rbind(outTab,cbind(gene=i,conMean=conGeneMeans,treatMean=treatGeneMeans,logFC=logFC,pValue=pvalue))        }}           #输出所有基因的差异情况write.table(outTab,file="all.txt",sep="t",row.names=F,quote=F)           #输出差异表格outDiff=outTab[(abs(as.numeric(as.vector(outTab$logFC)))>logFCfilter & as.numeric(as.vector(outTab$pValue))     <pvalFilter),]< span></pvalFilter),]<>write.table(outDiff,file="diff.xls",sep="t",row.names=F,quote=F)write.table(outDiff,file="diff.txt",sep="t",row.names=F,quote=F)           #输出差异基因的表达文件heatmap=rbind(ID=colnames(data[as.vector(outDiff[,1]),]),data[as.vector(outDiff[,1]),])write.table(heatmap,file="diffGeneExp.txt",sep="t",col.names=F,quote=F)           四、文件的导入  依然先将工作路径进行修改设置setwd("F:\ 96glycolysis\09.glyGeneExp\m6a\8")          rt=read.table("risk.txt",sep="t",header=T,row.names=1,check.names=F)  #读取输入文件rt=rt[order(rt$riskScore),]                                            #按照riskScore对样品排序#使用read.table()函数读取名为"risk.txt"的数据文件，并存储在变量rt中。

#风险曲线riskClass=rt[,"risk"]#从数据框rt中提取名为"risk"的列，将其赋值给riskClass变量。这个列应该包含风险等级的信息。lowLength=length(riskClass[riskClass=="low"])#计算风险等级为"low"的数量，将结果赋值给lowLength变量。highLength=length(riskClass[riskClass=="high"])#计算风险等级为"high"的数量，将结果赋值给highLength变量。line=rt[,"riskScore"]    #从数据框rt中提取名为"riskScore"的列，将其赋值给line变量。这个列应该包含风险评分的信息。line[line>10]=10     #将line中大于10的值设为10，限制风险评分的范围。(小果设定的值是10，小伙伴们可根据自己的数据集进行修改~)pdf(file="riskScore.pdf",width = 10,height = 5)#设置宽度、高度。（小伙伴们可根据自己需求修改参数）plot(line,     type="p",     pch=20,     xlab="Patients (increasing risk socre)",     ylab="Risk score",     col=c(rep("green",lowLength),     rep("red",highLength)))#绘制散点图，横坐标为患者，风险评分递增，纵坐标为风险评分，根据风险等级为"low"和"high"的数量来设定颜色（小伙伴们可根据自己的喜好修改~）#plot()是R语言中非常常用且功能强大的绘图内置函数，用于绘制各种类型的图形，如散点图、折线图、直方图等abline(v=lowLength,lty=5)#在图中添加一条垂直虚线legend("topleft",c("Highrisk","lowRisk"),bty="n",pch=19,col=c("red","green"),cex=1.2)#在左上角添加图例，分别表示"Highrisk"和"lowRisk"，颜色为红色和绿色。dev.off()#关闭PDF文件绘图设备，保存图表

#生存状态图color=as.vector(rt$fustat)color[color==1]="red"color[color==0]="blue"#创建一个向量color，其中的值与rt$fustat中的值相同。接着根据rt$fustat的值，将color向量中为1的元素设为"red"，为0的元素设为"blue"。pch=as.vector(rt$fustat)pch[pch==1]=15pch[pch==0]=19#创建一个向量pch，其中的值与rt$fustat中的值相同。接着根据rt$fustat的值，将pch向量中为1的元素设为15，为0的元素设为19。这里，小果说一下，pch是用于指定数据点的形状哦~pdf(file="survStat.pdf",width = 10,height = 5)#同样的，设置宽和高。（可根据自己需求修改参数哦~）plot(rt$futime,     pch=pch,     xlab="Patients (increasing risk socre)",     ylab="Survival time (years)",         col=color)#依然用到plot()函数绘制散点图，横坐标为rt$futime，纵坐标为rt$fustat，数据点的形状由pch向量指定，颜色由color向量指定。横坐标和纵坐标的标签分别为"Patients (increasing risk score)"和"Survival time (years)"。legend("topleft", c("Dead", "Alive"),bty="n",pch=c(15,19),col=c("red","blue"),cex=1.2)#在图中添加图例，位置在左上角，标签为"Dead"和"Alive"，对应的数据点形状为15和19，颜色分别为红色和蓝色。abline(v=lowLength,lty=2)#同样在图中添加一条垂直虚线，位置为lowLengthdev.off()#关闭并保存

#风险热图rt1=rt[c(3:(ncol(rt)-2))]#从数据框rt中选取第3列到倒数第3列的数据，并将其存储在新的数据框rt1中。（小伙伴们如果选用的数据集和小果不同的话记得要修改哦~）rt1=t(rt1)          #对数据框rt1进行转置，将行列交换，以便后续生成热图时行表示样本，列表示特征。rt1=log2(rt1+1)#对数据框rt1中的所有元素加1后取对数（以2为底），对数据进行对数变换，以便更好地展示数据的差异。annotation=data.frame(type=rt[,ncol(rt)])#创建一个名为annotation的数据框，其中包含一个名为type的列，列中的值来自rt数据框的最后一列。rownames(annotation)=rownames(rt)#将annotation数据框的行名设置为rt数据框的行名，以便后续在热图中添加注释信息。pdf(file="heatmap.pdf",width = 10,height = 5)#和上面的一样，设置宽和高。pheatmap(rt1,         annotation=annotation,          cluster_cols = FALSE,         cluster_rows = FALSE,         fontsize_row=11,         show_colnames = F,         fontsize=7,         fontsize_col=3,         color = colorRampPalette(c("green", "black", "red"))(50) )#用pheatmap()函数生成热图，其中：rt1——要展示的数据；annotation——包含了样本类别信息；cluster_cols和cluster_rows——是否对列和行进行聚类，这里小果选择了FALSE，即生成热图时，行和列的顺序将保持原始数据中的顺序，不会经过聚类算法重新排列；fontsize_row和fontsize_col——行和列标签的字体大小；show_colnames——是否显示列名，color——指定热图的颜色渐变，这里小果使用了从绿色到红色的渐变，共50个颜色。dev.off()#关闭并保存