手把手教你miRNA生存相关分析绘图
微小RNA的表达数据与患者的生存数据进行分析，以寻找与患者生存相关的miRNA。生存分析和绘制生存曲线，找出与患者生存相关的miRNA。这对于了解miRNA在癌症等疾病预后中的作用以及寻找潜在的生物标志物具有重要意义。生存分析是一种常用的统计分析方法，可用于评估不同基因或生物分子在患者预后中的潜在影响。通过该脚本，研究人员可以快速而全面地对大规模miRNA表达数据进行生存相关分析，从而辅助癌症等疾病的研究和诊断。
代码中对选定的相关miRNA的表达量进行生存分析。下面是要研究的多个mirna的样本表达量，以及下载的STAD癌症的生存时间和状态。
#if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") #BiocManager::install("limma") 检查是否已安装"BiocManager"包，如果没有安装，则通过install.packages函数安装它。 "BiocManager"包用于管理和安装生物信息学相关的其他包。 #install.packages('survival') #install.packages("survminer") 安装"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进行循环，找出预后相关的miRNA
outTab=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.
