生信分析第三步：生存曲線批量繪制

各位解螺旋的小伙伴大家好，我是先鋒宇，歡迎大家來到每周日的先鋒宇專欄，經過前兩期推文的學習，很多小伙伴都私信我說從先鋒宇助教的專欄很接地氣，自己能夠開始慢慢處理數據，并且希望先鋒宇助教能夠繼續把這條線走通。聽到解螺旋小伙伴積極正向的反饋，小編心理也是非常開心，那么今天咱們繼續往下走，我們在前兩期推文中完成數據的下載以及差異分析和單因素COX回歸，那么今天小編就帶大家進行生存曲線的繪制，先鋒宇助教還是本著科研效率第一的原則，我們當然不是一次就繪制一條生存曲線，這樣與我的風格不符，這次我就來教大家直接批量出圖，收獲滿滿的批量生存曲線圖。

寫在前面

相信很多小伙伴在看文獻的時候總是能夠看到作者拿一張figure放置多個生存曲線圖，不知道大家想過沒有，如果作者一張小圖一張小圖的畫，那可能圖還沒有畫完就直接開始拍桌子了。肯定為了提升科研的效率，我們還是希望把這些重復的工作都交給計算機不厭其煩地去做，然后留更多的時間給我們自己享受生活。好啦，言歸正傳，我們開始今天批量生存曲線繪制的講解。

代碼演示

經過前面兩期專欄的處理(沒有跟上的同學趕快去前兩期專欄看看，打牢基礎才能走得更遠)，我們現在已經得到了單因素COX回歸的結果。

接下來我們篩選p值小于0.05的基因進行保留，這里我們使用dplyr包中的filter函數進行過濾：

uniTab % dplyr::filter(pvalue < 0.05)

接著我們把單因素COX回歸有意義的基因再提取出來，因為剛剛得到的是數據庫，我們把第一列取出來即可：

unicox_gene

得到了單因素COX回歸有顯著統計學意義的基因之后，我們就要開始進行生存曲線繪制即K-M分析，關于K-M分析和COX回歸到底有啥區別，大家可以參考風師兄在生信下篇段位3有詳細的講解。如果用我自己的實用的理解那就是，COX分析篩選變量太多的話我們就再加上K-M分析再篩選一次，相當于雙重過濾標準；但是如果你COX回歸篩選之后就只有幾個基因了，那就沒有必要再用K-M分析去篩選了，因為篩選了完了你可能就沒有基因了。

接下來我們從單因素COX回歸的數據框中把pvalue小于0.05的基因提取出來，這里我們使用dplyr包中的select函數，注意這里要記得使用all_off函數把向量放在函數里面，這樣才能提取對應的列：

unicoxSig % dplyr::select(1,2,3,all_of(unicox_gene))unicoxSig$futime

接下來我們進行生存曲線的繪制，首先繪制生存曲線，我們首先需要解壓兩個強大的包，survival包和survminer包。

library(survival)library(survminer)

首先為了降低難度，我們先來進行一條生存生存曲線的繪制，我們先提取一個基因的表達量：

single_gene % dplyr::select(1,2,3,4)

然后我們構建一個分組文件，根據基因表達量的中位值進行高低兩組的劃分：

group median(single_gene[[4]]), "high", "low")

然后我們計算高低表達兩組之間的生存的p值大小：

diff=survdiff(Surv(futime, fustat) ~ group2, data = gene_surv)pValue=1-pchisq(diff$chisq,df=1)if(pValue<0.001){ pValue="p<0.001"}else{ pValue=paste0("p=",sprintf("%.03f",pValue))}

接下來擬合一個生存函數，這里我們使用survfit函數進行擬合：

fit

最后我們使用ggsurvplot函數來繪制生存曲線，代碼參考來自于生信體系課下篇，需要進一步學習的同學可以參看我們的生信體系課，里面有更多豐富的知識等待大家。

ggsurvplot(fit, data=single_gene, conf.int=TRUE, pval=pValue, pval.size=5, legend.labs=c("High", "Low"), xlab="Time (years)", ylab="Overall survival", break.time.by = 1, risk.table., palette=c("#d7191c", "#2b83ba"), risk.table=T, risk.table.height=.25)

一張可用于文章發表的生存曲線圖就繪制好了：

雖然一張繪制好了，但是我們本期的問題還沒有解決，我們不僅要一張，我們要很多張。套用一句經典的話就是，只要小孩子才做選擇，成年人是全部都要，而且越多越好，圖多了我們才有選擇的余地。

接下來我們進行批量繪制，批量繪制的原理無非就是循環，而循環就是一列一列循環，然后每一列繪制一個生存曲線。

for(gene in colnames(unicoxSig)[4:ncol(unicoxSig)]){group median(unicoxSig[[gene]]),???????????????"high",?"low")diff=survdiff(Surv(futime, fustat) ~ group, data = unicoxSig)pValue=1-pchisq(diff$chisq,df=1)if(pValue<0.001){ pValue="p<0.001"}else{ pValue=paste0("p=",sprintf("%.03f",pValue))}fit?surPlot=ggsurvplot(fit, data=unicoxSig, conf.int=TRUE, pval=pValue, pval.size=5, legend.labs=c("High", "Low"), legend.title=gene, xlab="Time (years)", ylab="Overall survival", break.time.by = 1, risk.table., palette=c("#d7191c", "#2b83ba"), risk.table=T, risk.table.height=.25)pdf(file=paste0("surv/",gene,".pdf"), onefile=FALSE, width=6.5, height=5.5)print(surPlot)dev.off()}

寫在最后

先鋒宇助教每次在跑循環的時候總感覺就是在收獲財富，因為每跑一張圖，就有可能放到論文里面，然后構成一個完成的figure，希望大家也能和我有同樣的感受！希望繼續繼續關注挑圈聯靠公眾號，繼續關注我的專欄，希望大家都能在這里學有所獲，收到滿滿的干貨，好啦，這期的內容就到這里啦，我們下周日再見~

往期傳送門：

讓生信工作者失業的神器——DrBioRight，真舍不得告訴你！

高效數據清洗！這個R包太強大了！你一定要試試！(文末附贈小彩蛋)

一站式分析R包來了，承包了生信各種分析！太全能了！

搞定這一步，說明你學R有天賦！TCGA數據從下載到差異分析(附代碼)

別走，baby，COX森林需要你

歡迎大家關注解螺旋生信頻道-挑圈聯靠公號~

—END—撰文丨先鋒宇排版丨四金兄值班 | 弘 ?毅主編丨小雪球

總結

以上是生活随笔為你收集整理的r 保留之前曲线_生存曲线居然能够批量绘制了的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。