TCGA數(shù)據(jù)挖掘九:lasso回歸

加載表達(dá)矩陣和生存數(shù)據(jù)

rm(list=ls())
options(stringsAsFactors = F)

Rdata_dir='Rdata/'
Figure_dir='figures/'
# 加載上一步從RTCGA.miRNASeq包里面提取miRNA表達(dá)矩陣和對(duì)應(yīng)的樣本臨床信息。
load( file = 
        file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata')
)
dim(expr)
dim(meta)
# 可以看到是 537個(gè)病人入宦,但是有593個(gè)樣本蛙埂,每個(gè)樣本有 552個(gè)miRNA信息。
# 當(dāng)然驰弄,這個(gè)數(shù)據(jù)集可以下載原始測(cè)序數(shù)據(jù)進(jìn)行重新比對(duì)麻汰,可以拿到更多的miRNA信息

# 這里需要解析TCGA數(shù)據(jù)庫(kù)的ID規(guī)律,來(lái)判斷樣本歸類(lèi)問(wèn)題戚篙。
group_list=ifelse(as.numeric(substr(colnames(expr),14,15)) < 10,'tumor','normal')

table(group_list)
exprSet=na.omit(expr)

table(group_list)
exprSet=na.omit(expr)
dim(exprSet)
load(  file = 
         file.path(Rdata_dir,'TCGA-KIRC-miRNA-survival_input.Rdata')
)
dim(exprSet) ## remove the nomral
head(phe)
exprSet[1:4,1:4]
head(colnames(exprSet))
head(phe$ID)
## 必須保證生存資料和表達(dá)矩陣五鲫,兩者一致
all(substring(colnames(exprSet),1,12)==phe$ID)

進(jìn)行l(wèi)asso回歸,找到最佳建模位置

library(lars) 
library(glmnet) 
x=t(log2(exprSet+1))#歸一化
y=phe$event
#用基因的表達(dá)情況預(yù)測(cè)生死
model_lasso <- glmnet(x, y, family="binomial", nlambda=50, alpha=1)#拉手回歸模型
print(model_lasso)
# 列%Dev代表了由模型解釋的殘差的比例岔擂,對(duì)于線性模型來(lái)說(shuō)就是模型擬合的R^2(R-squred)位喂。
# 它在0和1之間浪耘,越接近1說(shuō)明模型的表現(xiàn)越好,
# 如果是0塑崖,說(shuō)明模型的預(yù)測(cè)結(jié)果還不如直接把因變量的均值作為預(yù)測(cè)值來(lái)的有效七冲。
QQ截圖20190715194723.jpg
head(coef(model_lasso, s=c(model_lasso$lambda[29],0.009)))
#找到最有意義的一個(gè)點(diǎn)納入模型,就是說(shuō)沒(méi)有把所有基因都放入模型里面规婆,只是找到了放入基因數(shù)使得模型最好的那個(gè)點(diǎn)
#這個(gè)29是用眼睛看來(lái)的澜躺,我們可以用下面第三個(gè)代碼劃線來(lái)看
plot.glmnet(model_lasso, xvar = "norm", label = TRUE)
image.png
plot(model_lasso, xvar="lambda", label=TRUE)
image.png
cv_fit <- cv.glmnet(x=x, y=y, alpha = 1, nlambda = 1000)
plot.cv.glmnet(cv_fit)
# 兩條虛線分別指示了兩個(gè)特殊的λ值:
c(cv_fit$lambda.min,cv_fit$lambda.1se) 
image.png

再用得到的最佳的位置去建模

#再用得到的最佳的位置去建模,lASSO可以防止過(guò)度擬合
model_lasso <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.1se)
lasso.prob <- predict(cv_fit, newx=x , s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
re=cbind(y ,lasso.prob)
#得到預(yù)測(cè)結(jié)果
dat=as.data.frame(re[,1:2])
colnames(dat)=c('event','prob')
dat$event=as.factor(dat$event)#畫(huà)圖時(shí)需要factor
library(ggpubr) 
p <- ggboxplot(dat, x = "event", y = "prob",
               color = "event", palette = "jco",
               add = "jitter")
#  Add p-value
p + stat_compare_means()
#得出預(yù)測(cè)結(jié)果
QQ截圖20190715201226.jpg
image.png

判斷預(yù)測(cè)結(jié)果的準(zhǔn)確性

#判斷預(yù)測(cè)結(jié)果的準(zhǔn)確性
library(ROCR)
library(glmnet)
library(caret)
# calculate probabilities for TPR/FPR for predictions
pred <- prediction(re[,2], re[,1])
perf <- performance(pred,"tpr","fpr")
performance(pred,"auc") # shows calculated AUC for model
plot(perf,colorize=FALSE, col="black") # plot ROC curve
lines(c(0,1),c(0,1),col = "gray", lty = 4 )
image.png

后面這段是補(bǔ)充抒蚜,可以不用看了掘鄙。。嗡髓。

fit <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.1se)#構(gòu)建模型
head(fit$beta)
#一倍SE內(nèi)的更簡(jiǎn)潔的模型,是22個(gè)miRNA
#fit <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
#head(fit$beta)# 這里是40個(gè)miRNA
choose_gene=rownames(fit$beta)[as.numeric(fit$beta)!=0]
#選出用于建模的那些基因
length(choose_gene)
myexpr=x[,choose_gene]
mysurv=phe[,c("days","event")]
mysurv$days[mysurv$days< 1] = 1 
# 詳細(xì)代碼參見(jiàn)這個(gè)網(wǎng)站https://github.com/jeffwong/glmnet/blob/master/R/coxnet.R#
fit <- glmnet( myexpr, Surv(mysurv$days,mysurv$event), 
              family = "cox") 
#用包自帶的函數(shù)畫(huà)圖
plot(fit, xvar="lambda", label = TRUE)
plot(fit, label = TRUE)
## 如果需要打印基因名操漠,需要修改函數(shù),這里不展開(kāi)器贩。

library(pheatmap) 
choose_matrix=expr[choose_gene,]
choose_matrix[1:4,1:4]
n=t(scale(t(log2(choose_matrix+1))))  #scale()函數(shù)去中心化和標(biāo)準(zhǔn)化
#對(duì)每個(gè)探針的表達(dá)量進(jìn)行去中心化和標(biāo)準(zhǔn)化
n[n>2]=2 #矩陣n中歸一化后颅夺,大于2的項(xiàng),賦值使之等于2(相當(dāng)于設(shè)置了一個(gè)上限)
n[n< -2]= -2 #小于-2的項(xiàng)蛹稍,賦值使之等于-2(相當(dāng)于設(shè)置了一個(gè)下限)
n[1:4,1:4]

## http://www.bio-info-trainee.com/1980.html
annotation_col = data.frame( group_list=group_list  )
rownames(annotation_col)=colnames(expr)

pheatmap(n,show_colnames = F,annotation_col = annotation_col,
         filename = 'lasso_genes.heatmap.png')

library(ggfortify)
df=as.data.frame(t(choose_matrix))
df$group=group_list
png('lasso_genes.pca.png',res=120)
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')+theme_bw()
dev.off()

## 也可以嘗試其它主成分分析的R包吧黄,視頻就不繼續(xù)沒(méi)完沒(méi)了的講解了。


library("FactoMineR")
library("factoextra")  
## 這里的PCA分析唆姐,被該R包包裝成一個(gè)簡(jiǎn)單的函數(shù)拗慨,復(fù)雜的原理后面講解。
dat.pca <- PCA(t(choose_matrix), graph = FALSE) #'-'表示“非”
fviz_pca_ind(dat.pca,repel =T,
             geom.ind = "point", # show points only (nbut not "text")只顯示點(diǎn)不顯示文本
             col.ind =  group_list, # color by groups 顏色組
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses 集中成橢圓
             legend.title = "Groups"
)

最后

感謝jimmy的生信技能樹(shù)團(tuán)隊(duì)奉芦!

感謝導(dǎo)師岑洪老師赵抢!

感謝健明、孫小潔声功,慧美等生信技能樹(shù)團(tuán)隊(duì)的老師一路以來(lái)的指導(dǎo)和鼓勵(lì)烦却!

文中代碼來(lái)自生信技能樹(shù)jimmy老師!

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個(gè)濱河市先巴,隨后出現(xiàn)的幾起案子其爵,更是在濱河造成了極大的恐慌,老刑警劉巖伸蚯,帶你破解...
    沈念sama閱讀 206,126評(píng)論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件摩渺,死亡現(xiàn)場(chǎng)離奇詭異,居然都是意外死亡剂邮,警方通過(guò)查閱死者的電腦和手機(jī)摇幻,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,254評(píng)論 2 382
  • 文/潘曉璐 我一進(jìn)店門(mén),熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái),“玉大人绰姻,你說(shuō)我怎么就攤上這事枉侧。” “怎么了狂芋?”我有些...
    開(kāi)封第一講書(shū)人閱讀 152,445評(píng)論 0 341
  • 文/不壞的土叔 我叫張陵棵逊,是天一觀的道長(zhǎng)。 經(jīng)常有香客問(wèn)我银酗,道長(zhǎng)辆影,這世上最難降的妖魔是什么? 我笑而不...
    開(kāi)封第一講書(shū)人閱讀 55,185評(píng)論 1 278
  • 正文 為了忘掉前任黍特,我火速辦了婚禮蛙讥,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘灭衷。我一直安慰自己次慢,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 64,178評(píng)論 5 371
  • 文/花漫 我一把揭開(kāi)白布翔曲。 她就那樣靜靜地躺著迫像,像睡著了一般。 火紅的嫁衣襯著肌膚如雪瞳遍。 梳的紋絲不亂的頭發(fā)上闻妓,一...
    開(kāi)封第一講書(shū)人閱讀 48,970評(píng)論 1 284
  • 那天,我揣著相機(jī)與錄音掠械,去河邊找鬼由缆。 笑死,一個(gè)胖子當(dāng)著我的面吹牛猾蒂,可吹牛的內(nèi)容都是我干的均唉。 我是一名探鬼主播,決...
    沈念sama閱讀 38,276評(píng)論 3 399
  • 文/蒼蘭香墨 我猛地睜開(kāi)眼肚菠,長(zhǎng)吁一口氣:“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼舔箭!你這毒婦竟也來(lái)了?” 一聲冷哼從身側(cè)響起蚊逢,我...
    開(kāi)封第一講書(shū)人閱讀 36,927評(píng)論 0 259
  • 序言:老撾萬(wàn)榮一對(duì)情侶失蹤层扶,失蹤者是張志新(化名)和其女友劉穎,沒(méi)想到半個(gè)月后时捌,有當(dāng)?shù)厝嗽跇?shù)林里發(fā)現(xiàn)了一具尸體怒医,經(jīng)...
    沈念sama閱讀 43,400評(píng)論 1 300
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡炉抒,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,883評(píng)論 2 323
  • 正文 我和宋清朗相戀三年奢讨,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 37,997評(píng)論 1 333
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡拿诸,死狀恐怖扒袖,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情亩码,我是刑警寧澤季率,帶...
    沈念sama閱讀 33,646評(píng)論 4 322
  • 正文 年R本政府宣布,位于F島的核電站描沟,受9級(jí)特大地震影響飒泻,放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜吏廉,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,213評(píng)論 3 307
  • 文/蒙蒙 一泞遗、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧席覆,春花似錦史辙、人聲如沸。這莊子的主人今日做“春日...
    開(kāi)封第一講書(shū)人閱讀 30,204評(píng)論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽(yáng)。三九已至生巡,卻和暖如春耙蔑,著一層夾襖步出監(jiān)牢的瞬間,已是汗流浹背孤荣。 一陣腳步聲響...
    開(kāi)封第一講書(shū)人閱讀 31,423評(píng)論 1 260
  • 我被黑心中介騙來(lái)泰國(guó)打工纵潦, 沒(méi)想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留,地道東北人垃环。 一個(gè)月前我還...
    沈念sama閱讀 45,423評(píng)論 2 352
  • 正文 我出身青樓邀层,卻偏偏與公主長(zhǎng)得像,于是被迫代替她去往敵國(guó)和親遂庄。 傳聞我的和親對(duì)象是個(gè)殘疾皇子寥院,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,722評(píng)論 2 345