|
有一些五六年前的學(xué)生們都成長為了各個生物信息學(xué)相關(guān)公司的小領(lǐng)導(dǎo),，而且他們都有了自己的公眾號,，知乎號,，也算是一番人物,。最近他們跟我反饋面試找不到或者說很難直接考核篩選到認(rèn)真干活的生信工程師,，挺有意思的,。讓我想起來了早在生信技能樹論壇創(chuàng)立之初我為了引流,，而規(guī)劃的200個生信工程師面試題。值得繼續(xù)分享： 生存分析簡介在生物醫(yī)學(xué)研究中,，生存分析是非常重要和常見的分析方法,。生存分析經(jīng)常用在癌癥等疾病的研究中，例如在對某種抗癌藥物做臨床實驗時,，會首先篩選一部分癌癥患者隨機(jī)分為兩組,，一組服用該實驗藥物,，一組用對照藥物,，服藥之后開始統(tǒng)計每個患者從服藥一直到死亡的生存時間，通過考察兩組之間的病人在生存時間上是否有統(tǒng)計學(xué)差異來判斷試驗藥物是否有效,。在這種情況下,，死亡是整個實驗中重點觀測的事件,，即event。對于每個病人,，需要記錄它們發(fā)生該事件的具體時間,。因此，生存分析可以抽象的概述為,，研究在不同情況下,，特定時間發(fā)生的時間的關(guān)系是否存在差異。這些具體事件可以是死亡,，也可以是腫瘤轉(zhuǎn)移,、復(fù)發(fā)、病人出院,、重新入院等任何可以明確識別的事件,，而不同條件即為不同的分組依據(jù)，可以是年齡,、性別,、地獄、某個基因表達(dá)量的高低,、某個突發(fā)的攜帶是否等等,。兩個重要的概念： - 生存概率：即Survival probability，指的是研究對象從試驗開始知道某個時間點仍然存活的概率,，可見它是一個對時間t的函數(shù),，我們將其定義為S(t)
- 風(fēng)險概率：Hazard probability，指的是研究對象從試驗開始到某個特定時間t之間的存活,，但在t時間點發(fā)生觀測事件如死亡的概率,，它也是對時間t的函數(shù)，定義為H(t),。
我覺得這個關(guān)于生存分析的教程寫得真心不錯,，像我這樣第一次接觸生存分析的人也能理解：http://thisis./blog/index.php/2020/04/06/survival-analysis/ 首先使用模擬數(shù)據(jù)來了解下生存分析和風(fēng)險分析生存分析 目的，挑選出表達(dá)高低與生存差異存在顯著相關(guān)的基因 #首先生存模擬數(shù)據(jù),，50個病人,，200個基因
dat = matrix(rnorm(10000,mean = 8,sd=4),nrow = 200)
rownames(dat) = paste0('gene_',1:nrow(dat)) #行為基因
colnames(dat) = paste0('sample_', 1:ncol(dat)) #列
os_years = abs(floor(rnorm(ncol(dat),mean = 50,sd=20))) #年齡
os_status = sample(rep(c(0,1),25)) #生存狀態(tài),0代表死亡，1代表存活
library(survival)
#創(chuàng)建一個生存對象,，
mysurv = Surv(os_years,os_status )
## [1] 40 11 24+ 5 52+ 31+ 40 46 47+ 57 41+ 77+ 48+ 5+ 47+ 25 5 29+ 73+
#根據(jù)公式（Kaplan-Meier）,、之前擬合的Cox模型或者先前擬合的加速失效時間模型創(chuàng)建生存曲線
fit.KM = survfit(mysurv~1)
#根據(jù)基因的表達(dá)情況進(jìn)行分組，高表達(dá)組和低表達(dá)組,，然后確定每個基因的高表達(dá)組和低表達(dá)組的病人之間是否存在生存的顯著差異,，計算P值
log_rank_p <- apply(dat, 1, function(values1){
group=ifelse(values1> median(values1), 'high', 'low')
kmfit2 <- survfit(mysurv~group)
#plot(kmfit2)
data.survdiff <- survdiff(mysurv~group)
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) -1)
})
#挑選出表達(dá)高低能夠顯著區(qū)分生存的基因，這些基因可以用來繪制生存曲線圖
names(log_rank_p[log_rank_p<0.05])
挑選差異最顯著的兩個基因進(jìn)行生存分析繪圖 #將p值設(shè)置為0.01的時候,，只能篩選出兩個基因 ,，我們使用這兩個基因進(jìn)行分析
diff_genes<-names(log_rank_p[log_rank_p<0.01])
new_dat <- as.data.frame(t(dat[diff_genes,]))
new_dat$DSS.time <- os_years*30
new_dat$DSS <- os_status
new_dat$patient <- rownames(new_dat)
library(survminer)
#確定連續(xù)表達(dá)量的最優(yōu)切割點
new.cut <- surv_cutpoint(
new_dat,
time = 'DSS.time',
event = 'DSS',
variables = c("gene_37","gene_91")
)
summary(new.cut)
plot(new.cut,'gene_37',palette = 'npg')
#以切割位點為閾值,，將信息分為high，low
new.cut <- surv_categorize(new.cut)#使用兩個基因顎高低表達(dá)信息來判斷是否與病人的生存有關(guān)
library(survival)
fit <- survfit(Surv(DSS.time,DSS) ~ gene_37+gene_91, data = new.cut)
ggsurvplot(fit,pval = TRUE, conf.int = FALSE, title = 'METABRIC',
xlim = c(0,10000), break.time.by = 2000, risk.table.y.text.col=T,
risk.table.y.text = FALSE)
 
風(fēng)險分析 目的：挑選出能夠與患病風(fēng)險有相關(guān)性的基因 gender= ifelse(rnorm(ncol(dat))>1, 'male','female')
age = abs(floor(rnorm(ncol(dat), mean = 50, sd=20)))
cox_results <- apply(dat, 1, function(values1){
group = ifelse(values1 > median(values1), 'high','low')
survival_dat <- data.frame(group=group, gender=gender, age= age, stringsAsFactors = F)
#擬合比例風(fēng)險回歸模型,，使用Andersen和Gill技術(shù)公示將時間因變量,，時間因式，每個主題的多個時間以及其他或者的內(nèi)容都納入計算
m = coxph(mysurv~age + gender + group, data = survival_dat)
beta <- coef(m)
se = sqrt(diag(vcov(m)))
HR <- exp(beta)
HRse = HR * se
tmp = round(cbind(coef = beta, se=se, z = beta/se,
p = 1-pchisq((beta/se)^2,1),
HR = HR, HRse= HRse,
HEz = (HR-1)/HRse, HRp=1-pchisq(((HR-1)/HRse)^2, 1),
HRCILL = exp(beta - qnorm(.975, 0,1)*se),
HRCIUL = exp(beta + qnorm(.975, 0, 1)*se)),3)
return(tmp['grouplow',])
})
#挑選出能夠能夠顯著預(yù)測是否有致病風(fēng)險的基因
cox_results[,cox_results[4,]<0.05]
使用CGDS數(shù)據(jù)庫的乳腺癌表達(dá)數(shù)據(jù)進(jìn)行生存分析library(cgdsr)
#install.packages("cgdsr")
#DT是DataTables的R接口,，使得R的數(shù)據(jù)對象可以在HTML界面中以表格的形式展示
library(DT)
#cgdsr是一個用來訪問MSKCC Cancer Genomics Data Server（CGDS）中的數(shù)據(jù)的包
#從CGDS的終端URL創(chuàng)建一個連接對象
mycgds = cgdsr::CGDS("http://www./")
#為CGDA函數(shù)的調(diào)用設(shè)置詳細(xì)的日志級別
setVerbose(mycgds, TRUE)
#從CGDS中獲取可用的癌癥研究
all_TCGA_studies = getCancerStudies(mycgds)
#從CGDS中查詢到的所有可用的數(shù)據(jù)信息以網(wǎng)頁的形式展示
DT::datatable(all_TCGA_studies)
#從某個特定的癌癥研究中獲得可用的遺傳數(shù)據(jù)
all_gen_pro = getGeneticProfiles(mycgds, 'brca_metabric')
#獲得特定癌癥研究中可用的病例列表
DT::datatable(all_gen_pro)
all_tables <- getCaseLists(mycgds, 'brca_metabric')
#通過查看all_gen_pro可以看到很多可用的數(shù)據(jù),，我們挑選mRNA的表達(dá)數(shù)據(jù)
my_dataset <- 'brca_metabric_mrna'
my_table <- 'brca_metabric_mrna'
#檢索基因以及遺傳特征的基因組數(shù)據(jù),第二個參數(shù)為目標(biāo)基因，是必須要求的,，這里以
#"BCL11A","MBNL1"兩個基因為例
exp<-getProfileData(mycgds,c("BCL11A","MBNL1"), my_dataset,my_table)
exp$patient = rownames(exp)
#獲取臨床信息
cil_dat = getClinicalData(mycgds, my_table)
head(cil_dat)
DT::datatable(cil_dat,
extensions = "FixedColumns",
options = list(
scrollX=TRUE,
fixedColumns = TRUE
))
#查看臨床信息,，然后取出重要的兩列
colnames(cil_dat)
cil_dat_1 <- cil_dat[,c("OS_MONTHS","VITAL_STATUS")]
cil_dat_1$patient <- rownames(cil_dat)
#將表達(dá)數(shù)據(jù)和臨床信息根據(jù)病人進(jìn)行合并
tmp <- merge(exp,cil_dat_1,by='patient')
tmp = na.omit(tmp)
#將OS_MONTHS替換為DSS.time, VITAL_STATUS替換為DSS
colnames(tmp)[4] = 'DSS.time'
colnames(tmp)[5] = 'DSS'
tmp$DSS = ifelse(tmp$DSS=='Died of Disease',1,0)
tmp$DSS.time = tmp$DSS.time*30
library(survminer)
#確定連續(xù)表達(dá)量的最優(yōu)切割點
tmp.cut <- surv_cutpoint(
tmp,
time = 'DSS.time',
event = 'DSS',
variables = c("BCL11A","MBNL1")
)
summary(tmp.cut)
plot(tmp.cut,'BCL11A',palette = 'npg')
#以切割位點為閾值，將信息分為high,，low
tmp.cut <- surv_categorize(tmp.cut)#使用兩個基因顎高低表達(dá)信息來判斷是否與病人的生存有關(guān)
library(survival)
fit <- survfit(Surv(DSS.time,DSS) ~ BCL11A+MBNL1, data = tmp.cut)
ggsurvplot(fit,pval = TRUE, conf.int = FALSE, title = 'METABRIC',
xlim = c(0,10000), break.time.by = 2000, risk.table.y.text.col=T,
risk.table.y.text = FALSE)
DT可以將很長的數(shù)據(jù)表格以網(wǎng)頁的形式進(jìn)行展示,， 生存分析結(jié)果如下： 
文末友情推薦
|
