查看: 514|回复: 10

[技术前沿] 5.7+生信文章复现(四):DESeq2差异基因筛选

[复制链接]

管理员

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

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

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

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


复现文章:

Frontiers in oncology, IF= 5.738,2022

前三期我们对文章使用到的两个公共数据库(TCGA+ICGC)肾透明细胞癌表达矩阵和临床数据进行了整理,下面就可以开始进入数据挖掘啦!本期为第一步:差异表达缺氧基因(DEHGs)的筛选。

DEHGs筛选步骤:
1. 分子特征数据库(MSigDB)获取缺氧相关基因列表;
2. DEseq2差异分析筛选DEGs(TCGA+ICGC两个队列);
3. 取二者交集即为差异缺氧相关基因(DEHGs)

一、缺氧相关基因列表获取

进入MSigDBhttp://www.gsea-msigdb.org/gsea/msigdb/index.jsp)首页,在左侧Human Collections处点击Search进入检索页。


通过关键词搜索,这里我们检索HALLMARK_HYPOXIA即可获取目标基因集。或网页直接进入:https://www.gsea-msigdb.org/gsea ... et/HALLMARK_HYPOXIA


Download gene set处下载所需格式整理即可。


或者也可以在Show members处展开基因列表,直接复制到excel。



二、DESeq2差异分析

PS:
文章中写着差异分析是使用的DESeq包,这里实操我们使用当下更为主流的DESeq2进行。

  1. #载入所需R包:
  2. library(stringr)
  3. library(dplyr)
  4. library(DESeq2)
  5. #1.TCGA表达矩阵DESeq2差异分析:
  6. #载入TCGA表达矩阵:
  7. load('01_exp.Rdata') #数据来自第一期整理
  8. View(exp) #TCGA表达矩阵
复制代码


  1. #低表达基因预过滤过滤:
  2. exp2 <- exp[apply(exp, 1, function(x) sum(x > 1) > 607*0.5), ]
  3. dim(exp2)
复制代码


  1. #矩阵自检,如果是character则需转换为numeric:
  2. exp3 <- apply(exp2,2,as.numeric)
  3. row.names(exp3) <- row.names(exp2)
  4. #癌/癌旁分组:
  5. group_list <- ifelse(as.numeric(str_sub(colnames(exp3),14,15)) < 10,'tumor','normal')
  6. #指定顺序(转化为因子):
  7. group_list <- factor(group_list,levels = c("normal","tumor")) #需指定对照组在前
  8. table(group_list)
复制代码


  1. #保存所需数据并清空工作环境:
  2. save(exp3,group_list,file = c('03_exp_group_TCGA.Rdata'))
  3. rm(list = ls())
  4. #DESeq2流程:
  5. #创建colData(样本名和分组对应):
  6. load('03_exp_group_TCGA.Rdata')
  7. colData <- data.frame(row.names = colnames(exp3),
  8. condition = group_list)
  9. head(colData)
复制代码


  1. #构建dds对象:
  2. dds <- DESeqDataSetFromMatrix(countData = exp3,
  3. colData = colData,
  4. design = ~ condition)
  5. #跑DESeq2差异分析流程:
  6. dds <- DESeq(dds)
  7. #从dds提取差异分析结果表:
  8. ##顺序指定:tumor在前,normal在后(即tumor vs mormal,肿瘤样本相较于正常样本基因的上下调);
  9. res <- results(dds, contrast = c("condition",rev(levels(group_list)))) #和因子顺序相反
  10. res
复制代码


  1. #按P值排序并转化为dataframe:
  2. RES<- as.data.frame(res[order(res$pvalue),] )
  3. head(RES)
复制代码


  1. #添加上下调基因分组标签:
  2. ##筛选阈值:p<0.05,|log2FC|>1
  3. DESeq2 <- RES
  4. DESeq2$group <- case_when(
  5. DESeq2$log2FoldChange > 1 & DESeq2$pvalue < 0.05 ~ "up",
  6. DESeq2$log2FoldChange < -1 & DESeq2$pvalue < 0.05 ~ "down",
  7. abs(DESeq2$log2FoldChange) <= 1 ~ "none",
  8. DESeq2$pvalue >= 0.05 ~ "none"
  9. )
  10. #查看差异基因及数量:
  11. up <- rownames(DESeq2)[DESeq2$group=="up"]#差异上调基因
  12. down <- rownames(DESeq2)[DESeq2$group=="down"]#差异下调基因
  13. diff <- c(up,down)#所有差异基因
  14. length(up)
  15. length(down)
复制代码


  1. #和HYPOXIA基因集(缺氧相关基因)取交集:
  2. hypoxia <- read.csv('HALLMARK_HYPOXIA.v7.5.1.csv',header = F)
  3. head(hypoxia)
  4. length(hypoxia$V1) #200个缺氧相关基因
复制代码


  1. diff_hypoxia <- intersect(diff,hypoxia$V1)
  2. length(diff_hypoxia) #TCGA队列筛选出80个差异缺氧基因(DEHG)
复制代码


  1. #保存差异分析结果并清空工作环境:
  2. DESeq2_TCGA <- DESeq2
  3. diff_hypoxia_TCGA <- diff_hypoxia
  4. save(DESeq2_TCGA,diff_hypoxia_TCGA,file = c('03_DESeq2_TCGA.Rdata'))
  5. rm(list = ls())
复制代码

下面我们继续对ICGC队列完成DESeq2差异分析的流程,流程代码相同。

  1. #2.ICGC表达矩阵DESeq2差异分析:
  2. #载入ICGC表达矩阵和分组文件(数据来自第二期整理):
  3. load('ICGC_exp_group2.Rdata')
  4. View(exp_sym) #ICGC表达矩阵
复制代码


  1. head(group) #样本分组文件
复制代码


  1. #低表达基因预过滤过滤:
  2. exp <- exp_sym
  3. exp2 <- exp[apply(exp, 1, function(x) sum(x > 1) > 136*0.5), ]
  4. #指定分组顺序(转化为因子):
  5. group_list <- group$group
  6. group_list <- factor(group_list,levels = c("normal","tumor")) #指定对照组在前
  7. table(group_list)
复制代码


  1. #保存所需数据并清空工作环境:
  2. save(exp2,group,group_list,file = c('03_exp_group_ICGC.Rdata'))
  3. rm(list = ls())
  4. #DESeq2流程:
  5. #创建样本名和分组对应的colData:
  6. load('03_exp_group_ICGC.Rdata')
  7. colData <- data.frame(row.names = colnames(exp2),
  8. condition = group_list)
  9. #构建dds对象:
  10. dds <- DESeqDataSetFromMatrix(countData = exp2,
  11. colData = colData,
  12. design = ~ condition)
  13. #DESeq2差异分析流程:
  14. dds <- DESeq(dds)
  15. #从dds提取差异分析结果表:
  16. res <- results(dds, contrast = c("condition",rev(levels(group_list))))
  17. #按P值排序并转化为dataframe:
  18. RES<- as.data.frame(res[order(res$pvalue),] )
  19. #添加上下调基因分组标签:
  20. ##筛选阈值:p<0.05,|log2FC|>1
  21. DESeq2 <- RES
  22. DESeq2$group <- case_when(
  23. DESeq2$log2FoldChange > 1 & DESeq2$pvalue < 0.05 ~ "up",
  24. DESeq2$log2FoldChange < -1 & DESeq2$pvalue < 0.05 ~ "down",
  25. abs(DESeq2$log2FoldChange) <= 1 ~ "none",
  26. DESeq2$pvalue >= 0.05 ~ "none"
  27. )
  28. head(DESeq2)
复制代码


  1. #查看差异基因及数量:
  2. up <- rownames(DESeq2)[DESeq2$group=="up"]
  3. down <- rownames(DESeq2)[DESeq2$group=="down"]
  4. diff <- c(up,down)
  5. length(up)
  6. length(down)
复制代码


  1. #和HYPOXIA基因集取交集:
  2. hypoxia <- read.csv('HALLMARK_HYPOXIA.v7.5.1.csv',header = F)
  3. diff_hypoxia <- intersect(diff,hypoxia$V1)
  4. length(diff_hypoxia) #ICGC队列筛选出91个差异缺氧基因(DEHG)
复制代码


  1. #保存差异分析结果并清空工作环境:
  2. DESeq2_ICGC <- DESeq2
  3. diff_hypoxia_ICGC <- diff_hypoxia
  4. save(DESeq2_ICGC,diff_hypoxia_ICGC,file = c('03_DESeq2_ICGC.Rdata'))
  5. rm(list = ls())
复制代码

三、差异缺氧相关基因(DEHGs)获取

  1. #TCGA和ICGC差异分析结果取交集:
  2. load('03_DESeq2_TCGA.Rdata')
  3. load('03_DESeq2_ICGC.Rdata')
  4. length(intersect(diff_hypoxia_TCGA,diff_hypoxia_ICGC)) #66个共有差异缺氧相关基因
复制代码


  1. DEHG <- intersect(diff_hypoxia_TCGA,diff_hypoxia_ICGC) #差异缺氧相关基因筛选完成

  2. #保存差异缺氧相关基因列表:
  3. save(DEHG,file = c('03_DEHG.Rdata'))
  4. rm(list = ls())
复制代码

这样我们就完成了差异缺氧相关基因的筛选工作,原文直接用筛出的基因进入下一步预后筛选。我们额外做一小点补充,因为可能有小伙伴想要学习用火山图或韦恩图可视化上述差异分析筛选结果。简单过一遍代码:

  1. ##补充内容:
  2. #载入所需数据:
  3. load('03_DESeq2_TCGA.Rdata')
  4. load('03_DESeq2_ICGC.Rdata')
  5. hypoxia <- read.csv('HALLMARK_HYPOXIA.v7.5.1.csv',header = F)

  6. #载入绘图R包:
  7. library(ggplot2)
  8. library(VennDiagram)

  9. #TCGA队列:
  10. #将group转换为因子指定绘图顺序;
  11. df <- DESeq2_TCGA
  12. df$'-log10(pvalue)' <- -log10(df$pvalue)
  13. df$group <- factor(df$group, levels = c("up","down","none"))

  14. #绘图;
  15. colnames(df)
  16. mycol <- c("#FF9999","#99CC00","gray80") #自定义颜色
  17. p <- ggplot(df,aes(x = log2FoldChange,y = -log10(pvalue), color = group)) + #建立映射
  18. geom_point(size=1.6) + #绘制散点
  19. scale_colour_manual(name="",values=alpha(mycol,0.9)) +
  20. scale_y_continuous(expand=expansion(add = c(2, 0)),
  21. limits = c(0, 300), breaks = c(0,50,100,150,200,250,300))+
  22. scale_x_continuous(limits = c(-10, 10), breaks = c(-10,-5,0,5,10)) + #限制坐标轴
  23. geom_hline(yintercept = c(-log10(0.05)),size = 0.7,color = "grey",lty = "dashed")+
  24. geom_vline(xintercept = c(-1,1),size = 0.7,color = "grey",lty = "dashed") + #阈值辅助线
  25. theme_classic()
  26. p
复制代码


  1. #ICGC队列(代码相同):
  2. #将group转换为因子指定绘图顺序;
  3. df <- DESeq2_ICGC
  4. df$'-log10(pvalue)' <- -log10(df$pvalue)
  5. df$group <- factor(df$group, levels = c("up","down","none"))

  6. #绘图;
  7. colnames(df)
  8. mycol2 <- c("#ffcd05","#6acdff","gray80") #自定义颜色
  9. p1 <- ggplot(df,aes(x = log2FoldChange,y = -log10(pvalue), color = group)) + #建立映射
  10. geom_point(size=1.6) + #绘制散点
  11. scale_colour_manual(name="",values=alpha(mycol,0.9)) +
  12. scale_y_continuous(expand=expansion(add = c(2, 0)),
  13. limits = c(0, 300), breaks = c(0,50,100,150,200,250,300))+
  14. scale_x_continuous(limits = c(-10, 10), breaks = c(-10,-5,0,5,10)) + #限制坐标轴
  15. geom_hline(yintercept = c(-log10(0.05)),size = 0.7,color = "grey",lty = "dashed")+
  16. geom_vline(xintercept = c(-1,1),size = 0.7,color = "grey",lty = "dashed") + #阈值辅助线
  17. theme_classic()
  18. p1
复制代码


  1. #韦恩图展示两个队列差异基因和缺氧基因共有交集:
  2. #生成差异基因list:
  3. genelist <- list(TCGA = c(rownames(DESeq2_TCGA)[DESeq2_TCGA$group =="up"],rownames(DESeq2_TCGA)[DESeq2_TCGA$group =="down"]),
  4. ICGC = c(rownames(DESeq2_ICGC)[DESeq2_ICGC$group =="up"],rownames(DESeq2_ICGC)[DESeq2_ICGC$group =="down"]),
  5. hypoxia = hypoxia$V1)
  6. str(genelist)
复制代码


  1. filename <- c("diff_venn.tiff")#输出图片的文件名/文件格式设定
  2. mycol3 <- c("#fb0007","#139177","#ed9e08")
  3. venn.diagram(genelist,
  4. filename = filename,
  5. col = mycol3,
  6. fill = mycol3,
  7. main = "Differential Expression Hypoxia Genes",
  8. main.cex = 1.5,
  9. alpha = 0.3)
复制代码


筛选出66个差异缺氧相关基因,over。

更多绘图教程可搜索关键词查看往期。好啦,今天的分享就到这里!

参考文献
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: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
118 小时

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

使用道具 举报

迅猛龙

Rank: 8Rank: 8

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

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

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
1
注册时间
2016.8.25
在线时间
99 小时

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

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.9.20 09:12:53 | 显示全部楼层
学习学习
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.2.12
在线时间
119 小时

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

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.3.26
在线时间
13 小时

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

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.9.20 13:55:16 | 显示全部楼层
加油第十八天
新的一天加油!
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

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

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
118 小时

发表于 2022.9.21 11:29:08 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

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

使用道具 举报

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

本版积分规则

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