加載表達(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)的有效七冲。
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)
plot(model_lasso, xvar="lambda", label=TRUE)
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)
再用得到的最佳的位置去建模
#再用得到的最佳的位置去建模,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é)果
判斷預(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 )
后面這段是補(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"
)