|
在剛剛進(jìn)入生信領(lǐng)域的時(shí)候,，我想做的事情就是三個(gè),， 第一知道任何我想研究的基因在組織中的表達(dá)情況， 第二,，我選的基因?qū)δ[瘤的生存有無(wú)影響,， 第三這個(gè)基因可能的作用是什么？ 這是來(lái)自臨床醫(yī)生的視角,，研究疾病,，最終希望能夠服務(wù)臨床，臨床離不開(kāi)診斷和治療,，假設(shè)一個(gè)基因的表達(dá)對(duì)腫瘤的預(yù)后有影響,，他很可能就是我的盤(pán)中餐,。 有一大堆網(wǎng)頁(yè)工具可以實(shí)現(xiàn)生存分析，但是你看看jimmy已經(jīng)寫(xiě)的帖子 都可以批量做生存分析了,，還要網(wǎng)頁(yè)工具干嘛,？ 一拳把人打翻在地，扶都扶不起來(lái),。 但是沒(méi)有辦法他是群主,。即使會(huì)隨時(shí)被朝陽(yáng)群眾扭送到派出所，他依然是群主,，他的排版有問(wèn)題,，但是架不住他的內(nèi)容有深度，群主喜聞樂(lè)見(jiàn),。 那么問(wèn)題來(lái)了,，上面的問(wèn)題也可以反過(guò)來(lái)問(wèn)，既然別人已經(jīng)準(zhǔn)備好了網(wǎng)頁(yè)工具給你用,，你還寫(xiě)代碼做什么,？,！ 四個(gè)字,，批量，自由
大氣一點(diǎn)就是高端定制,。 就這么簡(jiǎn)單,。 那開(kāi)始我們的表演： 今天我們要對(duì)TCGA里面的任意基因做生存分析，最關(guān)鍵的我們要批量做生存分析,，然后選取生存差異最顯著的基因,。 首先安裝需要的包： # Load the bioconductor installer.
source('https:///biocLite.R')
options(BioC_mirror='https://mirrors.ustc.edu.cn/bioc/')
# Install the main RTCGA package
biocLite('RTCGA')
# Install the clinical and mRNA gene expression data packages
biocLite('RTCGA.clinical')
biocLite('RTCGA.mRNA')
然后是加載包 library(RTCGA)
#了解數(shù)據(jù)
infoTCGA <> infoTCGA() #這個(gè)命令會(huì)返回一個(gè)數(shù)據(jù)框，可以知道有哪些數(shù)據(jù)可被下載
#獲得臨床數(shù)據(jù)：
# Create the clinical data
library(RTCGA.clinical)
clin <> survivalTCGA(BRCA.clinical) #到這里臨床部分的信息已經(jīng)獲得啦
得到數(shù)據(jù)后我們先看一下他是什么結(jié)構(gòu) [1] 'data.frame'
再看一下前面幾行數(shù)據(jù) times bcrpatientbarcode patient.vital_status 1 3767 TCGA-3C-AAAU 0 2 3801 TCGA-3C-AALI 0 3 1228 TCGA-3C-AALJ 0 4 1217 TCGA-3C-AALK 0 5 158 TCGA-4H-AAAK 0 6 1477 TCGA-5L-AAT0 0
簡(jiǎn)單說(shuō)就是三列,，TCAGid,，生存時(shí)間，和發(fā)生的事件 獲得gene表達(dá)數(shù)據(jù)： library(RTCGA.mRNA) #加載數(shù)據(jù)包
class(BRCA.mRNA) #查看數(shù)據(jù)類(lèi)型發(fā)現(xiàn)是個(gè)數(shù)據(jù)框
dim(BRCA.mRNA) #看一下數(shù)據(jù)維度發(fā)現(xiàn)有590個(gè)樣本,，17815個(gè)基因
BRCA.mRNA[1:5, 1:5] #看一下部分?jǐn)?shù)據(jù)樣子
bcr_patient_barcode ELMO2 CREB3L1 RPS11 PNMA1
1 TCGA-A1-A0SD-01A-11R-A115-07 0.5070833 1.43450 0.765000 0.52600 2 TCGA-A1-A0SE-01A-11R-A084-07 0.1814167 0.89075 0.716000 0.13175 3 TCGA-A1-A0SH-01A-11R-A084-07 0.4615000 2.25925 0.417125 0.32500 4 TCGA-A1-A0SJ-01A-11R-A084-07 0.8770000 0.43775 0.115000 0.75775 5 TCGA-A1-A0SK-01A-12R-A084-07 1.4123333 -0.63725 0.492875 0.94325
好了,，行是樣本，列為gene 下面我們挑選幾個(gè)基因的表達(dá)數(shù)據(jù)出來(lái),，融合到生存的數(shù)據(jù)上去,，因?yàn)楸磉_(dá)量數(shù)據(jù)中TCGA的id號(hào)要長(zhǎng)一點(diǎn) 所以在融合前，需要先裁剪一下,，為了避免產(chǎn)生過(guò)多的中間變量,，我們使用管道符號(hào)%>%,他的作用是 把前一個(gè)計(jì)算得到結(jié)果，作為第二個(gè)函數(shù)的參數(shù),，示例如下： library(dplyr)
exprSet <> BRCA.mRNA %>%
# then make it a tibble (nice printing while debugging)
as_tibble() %>%
# then get just a few genes,這里是測(cè)試用
select(bcr_patient_barcode, PAX8, GATA3, ESR1) %>%
# then trim the barcode (see head(clin), and substr)
mutate(bcr_patient_barcode = substr(bcr_patient_barcode, 1, 12)) %>%
# then join back to clinical data
inner_join(clin, by='bcr_patient_barcode')
看一下數(shù)據(jù),，發(fā)現(xiàn)就是表達(dá)量加上time,，event兩列，所以記住這規(guī)律,，最終批量的時(shí)候需要減掉這兩列
 mark
開(kāi)始做生存分析 library(survival)
library(survminer)
對(duì)需要做生存分析的樣本分組,，把連續(xù)變量變成分類(lèi)變量，這里選擇測(cè)試的基因是GATA3 group <> ifelse(esprSet$GATA3>median(esprSet$GATA3),'high','low')
構(gòu)建生存對(duì)象,，并且進(jìn)行數(shù)據(jù)處理,，這里是兩步，我只是融合在了一起 sfit <> survfit(Surv(times, patient.vital_status)~group, data=esprSet)
sfit
summary(sfit)
直接繪圖 ggsurvplot(sfit, conf.int=F, pval=TRUE)
 mark
單個(gè)基因繪制成功,，我看一下能不能小規(guī)模循環(huán)操作 首先得到得到生存對(duì)象my.surv my.surv <> Surv(esprSet$times, esprSet$patient.vital_status)
使用apply循環(huán),對(duì)數(shù)據(jù)的2到4列進(jìn)行操作,，實(shí)際上就是PAX8, GATA3, ESR1這三個(gè)基因， 這里面使用了survdiff,，用來(lái)比較差異大小,，獲得p值 log_rank_p <> apply(esprSet[,2:4], 2, function(values1){
group=ifelse(values1>median(values1),'high','low')
kmfit2 <> survfit(my.surv~group,data=esprSet)
#plot(kmfit2)
data.survdiff=survdiff(my.surv~group)
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
})
運(yùn)行完了之后返回的找出小于0.05的P.val，在這個(gè)例子里面因?yàn)闆](méi)有基因的p.val小于0.05,，所以篩選不出來(lái) log_rank_p <> log_rank_p[log_rank_p0.05]
篩選后排序,，并獲得基因名 genediff <->rank_p))
好了生存數(shù)據(jù)已經(jīng)獲取， 表達(dá)數(shù)據(jù)小規(guī)模試驗(yàn)可行,， 批量操作也能實(shí)現(xiàn),， 下面就進(jìn)入實(shí)戰(zhàn)環(huán)節(jié)，前戲?qū)嵲谔L(zhǎng),，但是你要相信你只要不睡著,，這些都是值得的 獲得完整的表達(dá)量數(shù)據(jù) library(dplyr)
esprSet <> BRCA.mRNA %>%
# then make it a tibble (nice printing while debugging)
as_tibble() %>%
# then trim the barcode (see head(clin), and ?substr)
mutate(bcr_patient_barcode = substr(bcr_patient_barcode, 1, 12)) %>%
# then join back to clinical data
inner_join(clin, by='bcr_patient_barcode')
構(gòu)建生存對(duì)象my.surv library(survival)
my.surv <> Surv(esprSet$times, esprSet$patient.vital_status)
在進(jìn)行下一步之前，我居然突發(fā)奇想,，我想看一看,，這個(gè)表達(dá)數(shù)據(jù)里面哪些基因的NA值最多，有沒(méi)有NA值多過(guò)樣本數(shù)量一般的基因呢,？ 先構(gòu)建了一個(gè)函數(shù),，他對(duì)數(shù)據(jù)的列起作用，統(tǒng)計(jì)NA值的個(gè)數(shù),，最終返回成一個(gè)數(shù)據(jù)框 rem <> function(x){
r <>as.numeric(apply(x,2,function(i) sum(is.na(i))))
return(data.frame(geneName=names(x)[which(r > 0)],na_num=r[which(r > 0)]))
}
然后對(duì)表達(dá)量數(shù)據(jù)進(jìn)行統(tǒng)計(jì) 最終發(fā)現(xiàn)NA最多的基因是LCE1B,，有17個(gè)，所以數(shù)據(jù)不需要特殊處理啦 na_count <> dplyr::arrange(na_count,desc(na_num))
那就開(kāi)始批量運(yùn)算了,，一開(kāi)始就用apply,，發(fā)現(xiàn)大概需要運(yùn)行50分鐘以上，所以嘗試使用并行化處理 R語(yǔ)言里面的并行化有個(gè)專(zhuān)門(mén)的項(xiàng)目就是給apply的,，使用起來(lái)也是很方便 #嘗試使用并行運(yùn)算
library(parallel)
#detectCores()檢查當(dāng)前電腦可用核數(shù)
cl.cores <> detectCores()
#makeCluster(cl.cores)使用剛才檢測(cè)的核并行運(yùn)算,，我的服務(wù)器是28核56線程，我就用50吧
cl <> makeCluster(50)
#這是坑,，parApply里面用到的函數(shù)以及變量都需要申明,，不聲明就必須用模塊
clusterExport(cl,c('esprSet','my.surv'))
#length(names(esprSet))-2,，為什么減去2，因?yàn)橹靶∫?guī)模測(cè)試時(shí),，我們知道最后兩個(gè)是time和event,，不是表達(dá)量
#數(shù)據(jù)從25開(kāi)始，原因是從2開(kāi)始會(huì)報(bào)錯(cuò),，暫時(shí)無(wú)法解決,還有要注意是parApply,，A要大寫(xiě)的
log_rank_p <> parApply(cl,esprSet[,25:length(names(esprSet))-2],2,function(values){
group=ifelse(values>median(na.omit(values)),'high','low')
kmfit2 <> survival::survfit(my.surv~group,data=esprSet)
#plot(kmfit2)
data.survdiff=survival::survdiff(my.surv~group)
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
})
這個(gè)運(yùn)行的時(shí)間可能就是2分鐘 終止并行化 找出小于0.05的P.val log_rank_p <> log_rank_p[log_rank_p0.05]
篩選后排序，并獲得基因名 gene_diff <> as.data.frame(sort(log_rank_p))
最終得到的gene是2153個(gè),，這個(gè)數(shù)據(jù)是我之前留下的,，本次寫(xiě)貼時(shí)要帶孩子，在家沒(méi)法運(yùn)行運(yùn)算,。 數(shù)據(jù)保存可以這樣： save(gene_diff,file = 'gene_df.Rda')
如果想用的時(shí)候就這樣： load(file = 'gene_df.Rda')
沒(méi)有必要先寫(xiě)成txt格式,，然后要用的時(shí)候再讀取進(jìn)來(lái)，直接保存成R語(yǔ)言的格式即可,， 你能想象每次用完word保存成圖片,，下次再使用時(shí)用OCR識(shí)別圖片變成文字再編輯的狀態(tài)么？ 既然到了這一步,，我們隨便選取一個(gè)基因來(lái)作圖,，閉著眼睛都知道，都是有差異的 library(survminer)
group <> ifelse(esprSet$LRRC8D>median(esprSet$LRRC8D),'high','low')
sfit <> survfit(Surv(times, patient.vital_status)~group, data=esprSet)
ggsurvplot(sfit, conf.int=FALSE, pval=TRUE)
 mark好吧效果很不錯(cuò)嘛 這時(shí)候把癌和癌旁的數(shù)據(jù)作差異分析,，得到的基因與今天獲得的基因取交集,，就可以獲得又差異表達(dá),，又對(duì)生存有影響的基因了,。 今天我們是用的別人已經(jīng)下載好的數(shù)據(jù)，明天我們來(lái)嘗試自己下載并且清理數(shù)據(jù),，而且只用一種語(yǔ)言,，就是R語(yǔ)言 要知道，今天往后的ceRNA網(wǎng)絡(luò)構(gòu)建,，單基因GSEA都是基于這些表達(dá)量數(shù)據(jù)的,。 今天很美好，明天更有用,，但是大部分人都,。。,。,。馬云說(shuō)的。
|
