手把手教你miRNA生存相關分析繪圖
微小RNA的表達數據與患者的生存數據進行分析，以尋找與患者生存相關的miRNA。生存分析和繪制生存曲線，找出與患者生存相關的miRNA。這對于了解miRNA在癌癥等疾病預后中的作用以及尋找潛在的生物標志物具有重要意義。生存分析是一種常用的統(tǒng)計分析方法，可用于評估不同基因或生物分子在患者預后中的潛在影響。通過該腳本，研究人員可以快速而全面地對大規(guī)模miRNA表達數據進行生存相關分析，從而輔助癌癥等疾病的研究和診斷。
代碼中對選定的相關miRNA的表達量進行生存分析。下面是要研究的多個mirna的樣本表達量，以及下載的STAD癌癥的生存時間和狀態(tài)。
#if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") #BiocManager::install("limma") 檢查是否已安裝"BiocManager"包，如果沒有安裝，則通過install.packages函數安裝它。 "BiocManager"包用于管理和安裝生物信息學相關的其他包。 #install.packages('survival') #install.packages("survminer") 安裝"survminer"包，它是一個用于生存數據可視化和統(tǒng)計分析的包。 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表達數據進行預處理。首先，從樣本名中提取數字部分以便區(qū)分正常樣本#和腫瘤樣本。 #然后，根據"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列是生存狀態(tài)（0表示生存，1表示死亡）。接下來的代碼將讀取的生存數據進行處理： #提取第1列和第3列（生存時間和生存狀態(tài)），刪除可能包含缺失值的行， #并將列名改為"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列是生存狀#態(tài)（0表示生存，1表示死亡）。 #接下來的代碼將讀取的生存數據進行處理：提取第1列和第3列（生存時間和生存狀態(tài)），#刪除可能包含缺失值的行， #并將列名改為"futime"和"fustat"。最后，將生存時間轉換為年為單位（除以365）。
#對miRNA進行循環(huán)，找出預后相關的miRNA
outTab=data.frame() #創(chuàng)建一個空的數據框"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") #獲取最優(yōu)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()
}
#循環(huán)對每個miRNA進行生存分析，并繪制生存曲線。主要步驟如下：
?# 使用"for"循環(huán)遍歷從第3列開始到最后一列的所有miRNA的列名。在這里，每個miRNA
#的表達數據和生存數據都存儲在變量"rt"中
#?首先，檢查當前miRNA表達數據的標準差是否小于0.1，如果是，則跳過該miRNA的#分析，繼續(xù)下一個miRNA。
#提取當前miRNA的相關數據（"futime"為生存時間，"fustat"為生存狀態(tài)，"gene"為miRNA
#的表達數據）。
# 使用"surv_cutpoint"函數計算最優(yōu)截斷點，將患者根據miRNA的表達值分為高和低表達組。
#使用"surv_categorize"函數將數據進行分組并轉換為適合生存分析的格式。 "survfit"函數擬#合生存曲線。"survdiff"函數比較高低表達組之間的生存差異，得到p值。
#???根據p值的大小，將p值的顯示格式改為"p<0.001"或"p=xxx"（小數點后三位）。使用#"ggsurvplot"函數繪制生存曲線，其中包括p值的顯示、圖例等。將生存曲線保存為PDF文#件，文件名類似"sur.
