查看: 505|回复: 8

[技术前沿] 5.7+生信文章复现(五):单因素cox+lasso筛选预后相关DEGs

[复制链接]

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
806
注册时间
2020.6.16
在线时间
439 小时

发表于 2022.9.19 10:11:32 | 显示全部楼层 |阅读模式
无代码的TCGA纯生信文章复现学习:
《5+纯生信文章全图复现!无代码新手友好》

基于R语言的TCGA纯生信文章(预后风险模型构建)复现学习:
《5.7+生信文章复现(一):TCGA表达矩阵+临床信息清洗》
《5.7+生信文章复现(二):ICGC表达矩阵清洗》
《5.7+生信文章复现(三):ICGC临床信息清洗》
《5.7+生信文章复现(四):DESeq2差异基因筛选》​

复现文章:

Frontiers in oncology, IF= 5.738,2022

在前几期中,我们学习了两个公共数据库(TCGA+ICGC)表达矩阵和临床信息清洗,以及DESeq2筛选DEHGs(差异缺氧相关基因),今天我们学习如何从DEHGs中进一步筛选与预后相关的DEHGs。

预后DEHGs筛选流程:
1. 将TCGA总队列1:1拆分为训练集和测试集;
2. 单因素cox比例风险回归,筛选候选的预后DEHGs;
3. Lasso回归,十折交叉检验筛选最佳λ值,并获取最佳模型使用基因作为最终预后DEHGs,用于下一步预后模型构建。

  1. #相关R包和数据载入:
  2. library(stringr)
  3. library(caret)
  4. library(survminer)
  5. library(survival)
  6. #载入TCGA表达矩阵和临床信息表格(第1期整理):
  7. load('02_exp_cl_fpkm.Rdata') #第1期我们整理的是counts表达矩阵,这里我们使用fpkm,数据下载和整理方式完全相同,因此不再赘述
  8. #载入筛选的66个DEHGs(第4期整理):
  9. load('03_DEHG.Rdata')
  10. #矩阵自检,如果是character则需转换为numeric:
  11. exp2 <- apply(exp,2,as.numeric)
  12. row.names(exp2) <- row.names(exp)
  13. #构建差异缺氧相关基因DEHG表达矩阵:
  14. exp <- exp2[DEHG,]
  15. exp <- log2(exp+1) #取log(FPKM+1)进行后续分析
复制代码

一、临床信息清洗补充

在正式今天的内容前,我们在第一期整理的临床数据上稍微做点小补充。

  1. #OS单位从天转换为年:
  2. cl2 <- cl
  3. cl2$OS <- round(cl$OS/365,2) #单位年,保留2位小数
  4. #T satge去掉细分:
  5. cl2$`T stage` <- substr(cl$`T stage`,1,2)
  6. #将总分期中的小写转换为大写:
  7. cl2$`Clinical stage` <- gsub(pattern = "i", replacement = "I", cl2$`Clinical stage`)
  8. cl2$`Clinical stage` <- gsub(pattern = "v", replacement = "V", cl2$`Clinical stage`)
  9. cl <- cl2
  10. exp <- exp
  11. View(exp)
复制代码


  1. View(cl)
复制代码


  1. save(exp,cl,file = c('02_exp_cl_fpkm2.Rdata')) #保存一下更新后的所需数据
复制代码

补充完毕!

二、TCGA训练集和测试集队列拆分

  1. #1.训练集和测试集队列拆分(1:1):
  2. set.seed(123)
  3. inA <- createDataPartition(cl2$status, p = 0.5,times = 1,list = F)
  4. #训练集:
  5. train_exp <- exp[,inA]
  6. train_cl <- cl2[inA,]
  7. #测试集:
  8. test_exp <- exp[,-inA]
  9. test_cl <- cl2[-inA,]
  10. #切割后两个队列患者生死比例保持基本一致:
  11. prop.table(table(train_cl$status))
  12. prop.table(table(test_cl$status))
复制代码


三、单因素cox筛选预后相关候选DEHGs

用一个循环对矩阵进行批量单因素cox:

  1. cox <- apply(
  2. train_exp,1,function(x){
  3. train_cl$gene <- as.numeric(x)
  4. cox_genes <- coxph(Surv(OS, status) ~ gene, data = train_cl)
  5. coef <- coef(cox_genes) #回归系数
  6. SE <- sqrt(diag(vcov(cox_genes))) #标准误
  7. HR <- exp(coef) #风险比
  8. cox_need <- cbind(HR = HR,
  9. HR.95L = exp(coef - qnorm(.975, 0, 1) * SE),
  10. HR.95H = exp(coef + qnorm(.975, 0, 1) * SE),
  11. pvalue = 1 - pchisq((coef/SE)^2, 1))
  12. return(cox_need['gene',])
  13. }
  14. )
  15. unicox <- t(cox)
  16. head(unicox)
复制代码


  1. #提取预后差异基因列表(p<0.05):
  2. diff_unicox <- unicox[unicox[,4]<0.05,]
  3. dim(diff_unicox) #从66个差异基因中筛选出了25个候选预后DEHGs
  4. table(diff_unicox[,1]<1) #18个风险因子,7个保护因子
  5. head(diff_unicox)
复制代码


  1. #保存所需数据并清空工作环境:
  2. save(train_exp,train_cl,test_exp,test_cl,diff_unicox,file = '04_diff_unicox.Rdata')
  3. rm(list = ls())
复制代码

四、Lasso回归:十折交叉检验筛选最佳λ并获取最终DEHGs

  1. #继续载入lasso回归所需R包:
  2. library(glmnet)
  3. #载入训练集队列及候选预后DEHGs:
  4. load('04_diff_unicox.Rdata')
  5. #简化训练集队列命名:
  6. exp <- train_exp
  7. cl <- train_cl
  8. #构建候选预后DEHGs表达矩阵:
  9. exp <- exp[rownames(diff_unicox),]
  10. #使用lasso回归进一步收缩单因素cox筛选的25个候选DEHGs:
  11. x <- t(exp)
  12. y <- data.matrix(Surv(time = cl$OS,event = cl$status))
  13. x[1:6,1:6]
  14. head(y)
复制代码


  1. #构建模型:
  2. fit <- glmnet(x, y, family = 'cox', type.measure = "deviance", nfolds = 10)
  3. plot(fit,xvar = 'lambda',label = T) #候选DEHGs的lasso系数
复制代码


  1. #十折交叉检验筛选最佳lambda:
  2. set.seed(007)
  3. lasso_fit <- cv.glmnet(x, y, family = 'cox', type.measure = 'deviance', nfolds = 10)
  4. plot(lasso_fit)
复制代码


  1. #提取最佳λ值(这里选择1se对应lambda):
  2. lambda.1se <- lasso_fit$lambda.1se
  3. lambda.1se
复制代码


  1. #使用1se的lambda重新建模:
  2. model_lasso_1se <- glmnet(x, y, family = 'cox', type.measure = 'deviance', nfolds = 10,lambda = lambda.1se)
  3. #拎出建模使用基因:
  4. gene_1se <- rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]#as.numeric后"."会转化为0
  5. gene_1se #筛选出7个
复制代码


使用Lasso收缩后获得的7个基因作为最终预后相关DEHGs,用于下一步预后模型构建。

  1. #保存所需数据:
  2. save(gene_1se,file = c('05_lasso_gene_1se.Rdata'))
复制代码

更多TCGA相关内容可戳#TCGA数据库,今天的分享就到这里!

参考文献
Ning X, Li N, Qi Y, et al. Identification of a Hypoxia-Related Gene Model for Predicting the Prognosis and Formulating the Treatment Strategies in Kidney Renal Clear Cell Carcinoma[J]. Frontiers in oncology, 2022, 11: 806264.



*未经许可,不得以任何方式复制或抄袭本篇文章之部分或全部内容。版权所有,侵权必究。

本文作者:基迪奥-喵酱

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!

最近看过此主题的会员

回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2021.6.22
在线时间
61 小时

发表于 2022.9.19 14:14:42 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2019.4.24
在线时间
17 小时

发表于 2022.9.20 09:10:02 | 显示全部楼层
新的一天
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2016.5.10
在线时间
24 小时

发表于 2022.9.20 14:02:37 | 显示全部楼层
加油十八天
新的一天加油!
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2019.4.24
在线时间
17 小时

发表于 2022.9.21 09:39:28 | 显示全部楼层
来了学习
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2016.5.10
在线时间
24 小时

发表于 2022.9.21 10:21:58 | 显示全部楼层
加油十九天
新的一天加油!
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2019.4.24
在线时间
17 小时

发表于 2022.9.22 09:55:43 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2016.5.10
在线时间
24 小时

发表于 2022.9.22 21:21:19 | 显示全部楼层
加油第二十天
新的一天加油!
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2019.4.24
在线时间
17 小时

发表于 2022.9.23 09:27:32 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表