各位解螺旋的小夥伴們大家周日好，我是先鋒宇，歡迎大家來到每周日的先鋒宇專欄。上一期的推文我們打開了基因組學的大門，不知道大家對於好看又好吃的棒棒糖圖有沒有留下深刻的印象呢？棒棒糖雖好，但是也不要貪吃哦。今天小編繼續給大家帶來另外一份大餐，我們來進行生存曲線再一次的批量繪製，這一次我們把目光從轉錄組學移到基因組學，我們在基因組學進行基因突變生存曲線的批量繪製，小夥伴們有沒有很心動呢？那懷著愉悅的心情和小編一起來學習一下吧。
生存曲線的繪製對於生信分析是非常重要的，往往在一篇文章中我們看到肯定不止繪製一條生存曲線，而是有很多條生存曲線。先鋒宇老師很喜歡說一句話：機械重複的工作拿給機器來完成。所以我們應該擺脫機械重複的工作，而應該用更多的時間來享受生活。那麼今天我們就一起來看看基因組學的生存曲線批量繪製能不能提升我們的科研效率呢？
經過上一期推文的講解，我們已經從Xena下載了胃癌的基因組學數據以及過濾掉了不引起胺基酸改變的突變，得到了如下的數據框(這一步不清楚的同學可以去複習一下先鋒宇老師上一期的推文：小白重頭學習基因組學分析，從繪製棒棒糖圖開始)
那接下來為了繪製每個基因突變與否的批量生存曲線，我們首先需要去構建突變矩陣，指明基因在哪些樣本是突變的，哪些樣本是不突變的。所以我們首先需要提取上面這個數據框的前兩列(這裡不懂突變矩陣沒關係，看完了你就都明白了)
mutation_matrix <- verscan2_code %>% dplyr::select(1,2)可以看到這兩列數據是一個明顯的長數據，裡面只有樣本對應有突變的基因，所以我們先添加一列mutation，在後面標記上這些都是樣本裡面的有突變的基因。
mutation_matrix$group <- c("mutation")可以看到我們添加了一列信息來指明這些樣本中的基因是突變的，這在我們接下來構建突變矩陣的時候非常有用。接下來我們需要對突變矩陣去重，因為一個樣本的同一個基因發生不同的突變我們在突變矩陣裡面就只需要知道它是mutation，不需要管具體的突變方式。所以我們這裡使用unique函數進行去重複。
mutation_matrix <- unique(mutation_matrix)可以看到去重之後我們去掉了20000多行，那麼接下來我們就需要把這個長數據變成寬數據，也就是我們習慣的矩陣的形式，這樣才便於我們進行作圖。長寬數據轉換在R裡面用到的也比較多，還不會的小夥伴們也可以參考我之前的推文，我詳細介紹了長寬數據轉換的三種方法，這裡我們使用比較常用的pivot_wider函數進行轉換，函數作用和名字一樣，把數據變wider。
mutation_matrix1 <- pivot_wider(mutation_matrix, names_from = "Sample_ID", values_from = "group")可以看到變為寬數據之後，列名為樣本名，然後我剛剛加的一列標籤(mutation)在這裡就起作用啦，當數據變為寬數據之後，我們就能看到哪些基因在哪裡樣本裡面是突變的，而沒有突變的在這裡就是NA,那麼接下來我們再把NA全部變為wild(野生型)就可以啦。
mutation_matrix1 <- mutation_matrix1 %>% column_to_rownames("gene")mutation_matrix1[is.na(mutation_matrix1)] <- c("wild")可以看到突變矩陣就構建好了，那麼接下來我們為了要繪製生存曲線，我們接下來就需要把TCGA樣本的生存數據引入進來，前面先鋒宇專欄的推文已經講解了生存數據的下載和讀取，這裡就不加贅述啦，我們接下來讀取TCGA胃癌的生存數據。
stad_survival <- read_tsv("TCGA-STAD.survival.tsv.gz") %>% dplyr::select(1,2,4)colnames(stad_survival) <- c("sample", "fustat", "futime")stad_survival$futime <- stad_survival$futime/365可以看到我們把TCGA胃癌病人的生存信息讀入了，那麼接下來我們就把突變矩陣和生存信息進行整合，根據TCGA樣本講兩個數據框整理為一個數據框，我們這裡仍然不使用merge函數，我們使用更為簡單的inner_join函數來進行操作
mutation_matrix2 <- t(mutation_matrix1) %>% as.data.frame() %>% rownames_to_column("sample")mutation_survival <- inner_join(stad_survival, mutation_matrix2, by = "sample")可以看到我們把胃癌病人的生存信息與基因突變信息整合起來，這樣接下來我們就可以進行基因突變生存曲線的批量繪製啦。
首先為了降低難度，我們仍然先進行一條生存曲線的繪製，所以我們提取前面4列，只保留矩陣的第一個基因和生存信息。
gene_survival <- mutation_survival %>% dplyr::select(1,2,3,4)接下來我們進行生存曲線的繪製，我們首先打開survminer和survival包，方便我們調用survdiff函數來計算pvalue。
library(survminer)library(survival)diff=survdiff(Surv(futime, fustat) ~ gene_survival[[4]],data = gene_survival)pValue=1-pchisq(diff$chisq,df=1)接下來我們使用survfit來構建一個生存函數，並使用ggsurvplot函數來進行生存曲線繪製。
fit <- survfit(Surv(futime, fustat) ~ gene_survival[[4]],data = gene_survival)ggsurvplot(fit, data=gene_survival, conf.int=TRUE, pval=pValue, pval.size=6, legend.labs=c("mutation", "wild"), legend.title="CLSTN1", xlab="Time (years)", ylab="Overall survival", break.time.by = 1, risk.table.title="", palette=c("#d7191c", "#1a9641"), risk.table=T, risk.table.height=.25)可以看到一張可以用於發表的基因組的生存曲線就繪製好啦。
一張生存曲線繪製好了，但是這還遠遠不敢，成年人的世界裡我們不做選擇，而是全部都要，所以我們就要對生存曲線進行批量出圖，其實總的思想跟前面講的轉錄組也是一樣的，主要的區別就是這裡不是基因的表達量，而是基因的突變信息，同樣是按照列對每個基因分別進行生存曲線的繪製。
for (gene in colnames(mutation_survival[4:17413])) { diff=survdiff(Surv(futime, fustat) ~ mutation_survival[[gene]],data = mutation_survival) pValue=1-pchisq(diff$chisq,df=1) fit <- survfit(Surv(futime, fustat) ~ mutation_survival[[gene]],data = mutation_survival) surPlot = ggsurvplot(fit, data=gene_survival, conf.int=TRUE, pval=pValue, pval.size=6, legend.labs=c("mutation", "wild"), legend.title=gene, xlab="Time (years)", ylab="Overall survival", break.time.by = 1, risk.table.title="", palette=c("#d7191c", "#1a9641"), risk.table=T, risk.table.height=.25) pdf(file=paste0("survival/",gene,".pdf"), onefile=FALSE, width=6.5, height=5.5) print(surPlot) dev.off()}這裡我們使用for循環對於列進行循環，最終我們可以得到10000多個基因突變信息的生存曲線圖，大家只需要慢慢挑選自己想要的基因的生存曲線圖進行展示即可。
機械重複的工作我們應該儘可能交給計算機去完成，如果我們每天就是去做機械重複的工作，比如這裡我們去重複幾百次甚至上萬次一樣的操作去繪製想要的生存曲線，那這樣科研的效率就會大打折扣，有些時候學會適當的偷懶，把這些時間用來享受生活，這才是我們應該經歷的人生，而不是永遠沉浸在科研的無底洞裡面，看不到盡頭。
好啦，今天的生存曲線批量繪製就到這裡結束啦，希望這道大餐能給冬日冷冷的天氣裡面增加一絲暖意。我是先鋒宇，我們下周日繼續相約挑圈聯靠的先鋒宇專欄，我們不見不散。
往期傳送門：
讓生信工作者失業的神器——DrBioRight，真捨不得告訴你！
高效數據清洗！這個R包太強大了！你一定要試試！（文末附贈小彩蛋）
一站式分析R包來了，承包了生信各種分析！太全能了！
搞定這一步，說明你學R有天賦！TCGA數據從下載到差異分析（附代碼）
別走，baby，COX森林需要你
生存曲線居然能夠批量繪製了
ROC居然能夠批量篩選
TCGA轉錄組和臨床信息聯合分析，結果發現.
小白重頭學習基因組學分析，從繪製棒棒糖圖開始
歡迎大家關註解螺旋生信頻道-挑圈聯靠公號~
2021年國自然申請指南已公布
2021申請國自然還很迷茫？
國自然申請指南找不到重點？
酸談為你解讀最新國自然申請指南！
本周四酸談老談直播主題：
---2021年國自然申請指南解讀---
下方視頻號點擊預約直播
酸菜老談助你衝刺國自然！周四酸談直播基本信息
主題：《2021年國自然申請指南解讀》
日期：1月21日 18:00—20:00