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

#install.packages('survival')#install.packages("survminer")

library(limma)library(survival)library(survminer)
expFile="cor.MirnaExp.txt"            #共表达miRNA的表达文件cliFile="TCGA-STAD.survival.tsv"      #生存数据文件#指定了miRNA表达数据文件和生存数据文件并将它们分别赋值给变量"expFile"和"cliFile"#读取表达文件，并对输入文件整理data=read.table(expFile, header=T, sep="t", check.names=F, row.names=1)#使用read.table函数读取名为"cor.MirnaExp.txt"的miRNA表达文件。
#删掉正常样品group=sapply(strsplit(colnames(data),"\-"), "[", 4)group=sapply(strsplit(group,""), "[", 1)group=gsub("2", "1", group)data=t(data[,group==0,drop=F])rownames(data)=gsub("(.*?)\-(.*?)\-(.*?)\-(.*?)\-.*", "\1\-\2\-\3-\4", rownames(data))#对读取的miRNA表达数据进行预处理。首先，从样本名中提取数字部分以便区分正常样本#和肿瘤样本。#然后，根据"group"变量中的值，将正常样本对应的列从数据中删除。最后，根据正则表达#式，将行名的样本编号进行调整。#读取生存数据cli=read.table(cliFile, header=T, sep="t", check.names=F, row.names=1)cli=cli[,c(3,1)]cli=na.omit(cli)colnames(cli)=c("futime","fustat")cli$futime=cli$futime/365#读取名为"TCGA-STAD.survival.tsv"的生存数据文件。文件以制表符分隔，第一行是列名，#所以header=T表示第一行是列名。sep="t"表示文件中的字段是用制表符分隔的。check.names=F表示不检查列名的合法性。row.names=1表示使用第一列作为行名。#生存数据文件中的列意义分别是：第1列是样本ID，第3列是生存时间，#第2列是生存状态（0表示生存，1表示死亡）。接下来的代码将读取的生存数据进行处理：#提取第1列和第3列（生存时间和生存状态），删除可能包含缺失值的行，#并将列名改为"futime"和"fustat"。最后，将生存时间转换为年为单位（除以365）。#数据合并sameSample=intersect(row.names(data), row.names(cli))data=data[sameSample,,drop=F]cli=cli[sameSample,,drop=F]rt=cbind(cli, data)#读取名为"TCGA-STAD.survival.tsv"的生存数据文件。文件以制表符分隔，第一行是列名，所#以header=T表示第一行是列名。sep="t"表示文件中的字段是用制表符分隔的。check.names=F表示不检查列名的合法性。row.names=1表示使用第一列作为行名。
#生存数据文件中的列意义分别是：第1列是样本ID，第3列是生存时间，第2列是生存状#态（0表示生存，1表示死亡）。#接下来的代码将读取的生存数据进行处理：提取第1列和第3列（生存时间和生存状态），#删除可能包含缺失值的行，#并将列名改为"futime"和"fustat"。最后，将生存时间转换为年为单位（除以365）。#对miRNA进行循环，找出预后相关的miRNAoutTab=data.frame()#创建一个空的数据框"outTab"来保存每个miRNA的分析结果。for(gene in colnames(rt)[3:ncol(rt)]){  #提取miRNA信息  if(sd(rt[,gene])<0.1){next}  data=rt[,c("futime", "fustat", gene)]  colnames(data)=c("futime", "fustat", "gene")
  #获取最优cutoff,并对样品进行分组  res.cut=surv_cutpoint(data, time = "futime", event = "fustat", variables =c("gene"))  res.cat=surv_categorize(res.cut)  fit=survfit(Surv(futime, fustat) ~gene, data = res.cat)
  #比较高低表达组生存差异  diff=survdiff(Surv(futime, fustat) ~gene, data =res.cat)  pValue=1-pchisq(diff$chisq, df=1)  outVector=cbind(gene, res.cut$cutpoint[1], pValue)  outTab=rbind(outTab,outVector)  if(pValue<0.001){    pValue="p<0.001"  }else{    pValue=paste0("p=",sprintf("%.03f",pValue))  }
  #绘制生存曲线  surPlot=ggsurvplot(fit,             data=res.cat,             pval=pValue,             pval.size=6,             legend.title=gene,             legend.labs=c("high","low"),             xlab="Time(years)",             ylab="Overall survival",             palette=c("red", "blue"),             break.time.by=1,             conf.int=T,             risk.table=F,             risk.table.title="",             risk.table.height=.25)
  #输出图形  pdf(file=paste0("sur.",gene,".pdf"),onefile = FALSE,    width = 5.5,       #图片的宽度    height =5)         #图片的高度  print(surPlot)  dev.off()}#循环对每个miRNA进行生存分析，并绘制生存曲线。主要步骤如下： # 使用"for"循环遍历从第3列开始到最后一列的所有miRNA的列名。在这里，每个miRNA#的表达数据和生存数据都存储在变量"rt"中#  首先，检查当前miRNA表达数据的标准差是否小于0.1，如果是，则跳过该miRNA的#分析，继续下一个miRNA。#提取当前miRNA的相关数据（"futime"为生存时间，"fustat"为生存状态，"gene"为miRNA#的表达数据）。# 使用"surv_cutpoint"函数计算最优截断点，将患者根据miRNA的表达值分为高和低表达组。#使用"surv_categorize"函数将数据进行分组并转换为适合生存分析的格式。"survfit"函数拟#合生存曲线。"survdiff"函数比较高低表达组之间的生存差异，得到p值。#    根据p值的大小，将p值的显示格式改为"p<0.001"或"p=xxx"（小数点后三位）。使用#"ggsurvplot"函数绘制生存曲线，其中包括p值的显示、图例等。将生存曲线保存为PDF文#件，文件名类似"sur.<miRNA>.pdf"。 #   循环继续处理下一个miRNA。循环结束后，将"outTab"数据框按照p值的大小进行排#序，并将结果保存为名为"survival.txt"的制表符分隔文本文件。outTab=outTab[order(as.numeric(as.vector(outTab[,"pValue"]))),]write.table(outTab, file="survival.txt", sep="t", row.names=F, quote=F)