查看: 655|回复: 10

[R语言] 案例复现(一):表达数据的下载与差异分析

[复制链接]

管理员

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

主题
909
注册时间
2020.6.16
在线时间
485 小时

发表于 2022.10.26 09:48:54 | 显示全部楼层 |阅读模式
基于数据库进行数据挖掘,之前在《GEO数据库介绍与数据检索》《如何从GEO下载数据辅助文章发表?》 《表达谱数据的整理与差异分析》 等推文中已经为大家介绍过GEO数据库数据的查找、下载、分析等。

下面就以一篇今年发表的文章(如下图)为例,看下如何以具体的案例学习基于GEO数据库进行数据挖掘


从文献中找到GEO数据库的登录号(这里为GSE74089),在NCBI的GEO数据库进行搜索。点击查到的数据记录,可以查看该数据集相关的详细描述信息,如文章标题、发表时间、平台信息、样本信息(比如这里有8个样本)等。然后,我们只需点击蓝色的 “Analysis with GEO2R” 按钮即可开始对数据进行分析,如下图。


最终的分析结果如下,结果中给出了火山图、MA plot、Umap聚类图等常见的组学图表,最重要的是下面的差异基因信息表格,如下图。我们可以点击Download full table 下载完整的表格。当然,也可以点击Select columns 选择添加其他所需的注释信息,比如Gene symbol、GO ID、GO注释信息等。


使用GEO2R工具进行分析比较简单,但如何使用R脚本进行个性化分析呢?下面仍以文献中的数据集(GSE74089)为例,详细为大家演示如何使用R script 进行分析吧。

1. 数据下载

  1. #GEOquery包的安装方法如下:
  2. library(BiocManager)
  3. install("GEOquery")
  4. #加载所需R包;
  5. library(GEOquery)
  6. library(limma)
  7. library(umap)
  8. #使用getGEO函数下载GEO数据集和平台数据;
  9. gset <- getGEO("GSE74089", GSEMatrix =TRUE, AnnotGPL=TRUE)
  10. #查看基因集对象中的数据集个数;
  11. num <- length(gset)
  12. num
  13. #[1] 1
  14. #获取gset列表中与平台数据匹配的数据集;
  15. if (length(gset) > 1) {
  16. idx <- grep("GPL13497", attr(gset, "names"))
  17. }else {idx <- 1}
  18. #从list中提取匹配的数据集;
  19. gset <- gset[[idx]]
  20. #保存数据;
  21. save(gset,file = "GSE74089_data.Rdata")
  22. #载入数据;
  23. load("GSE74089_data.Rdata")
  24. #查看数据(列名)标签;
  25. fvarLabels(gset)
  26. #提取并调整变量标签便于用作矩阵列名,比如将":"、空格等转换成".";
  27. fvarLabels(gset) <- make.names(fvarLabels(gset))
  28. fvarLabels(gset)
复制代码


  1. #使用featureData函数获取基因注释信息;
  2. #主要包含varMetadata 和 data 两个文件,
  3. fD <- featureData(gset)
  4. #使用fData函数获取基因注释,即fd2中的data文件;
  5. fd <- fData(gset)
  6. #查看部分基因注释信息;
  7. fd[12:19,5:8]
复制代码


  1. #使用pData函数获取样本信息;
  2. pd <- pData(gset)
  3. #预览样本信息;
  4. pd[1:2]
复制代码


  1. #使用exprs函数提取基因表达量数据;
  2. ex <- exprs(gset)
  3. #预览数据;
  4. head(ex)
复制代码


  1. #计算表达量数据的分位数;
  2. qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
  3. qx
  4. #[1] 2.323333 3.562519 5.325553 6.968367 12.570155 17.036762
  5. #确定数据是否需要对数转换;
  6. LogC <- (qx[5] > 100) || (qx[6]-qx[1] > 50 && qx[2] > 0)
  7. LogC
  8. #[1] FALSE
  9. #对表达量数据进行log2转换(log2 transformation);
  10. #根据上面的结果,这里并没有进行log2转换;
  11. if (LogC) { ex[which(ex <= 0)] <- NaN
  12. exprs(gset) <- log2(ex) }
复制代码

2. 差异分析(基于limma)

  1. #根据样本信息指定分组组合(1为control组,0为patient组);
  2. gsms <- "10010011"
  3. #分割成单个字符;
  4. sml <- strsplit(gsms, split="")[[1]]
  5. #[1] "1" "0" "0" "1" "0" "0" "1" "1"
  6. #指定样本分组;
  7. gs <- factor(sml)
  8. gs
  9. #[1] 1 0 0 1 0 0 1 1
  10. #Levels: 0 1
  11. #生成分组标签;
  12. groups <- make.names(c("patient","control"))
  13. groups
  14. #设置分组水平;
  15. levels(gs) <- groups
  16. gs
  17. #[1] control patient patient control patient patient control control
  18. #Levels: patient control
  19. #在基因集对象中添加分组信息;
  20. gset$group <- gs
  21. #生成模型设计矩阵;
  22. design <- model.matrix(~group + 0, gset)
  23. #替换列名;
  24. colnames(design) <- levels(gs)
  25. #线性模型拟合 (fit linear model);
  26. fit <- lmFit(gset, design)
  27. #查看设计矩阵;
  28. design
复制代码


  1. #生成比较标签;
  2. cts <- paste(groups[1], groups[2], sep="-")
  3. #构建比较矩阵;
  4. cont.matrix <- makeContrasts(contrasts=cts, levels=design)
  5. cont.matrix
  6. #计算estimated coefficients和标准误;
  7. fit2 <- contrasts.fit(fit, cont.matrix)
  8. #差异分析(Empirical Bayes Statistics,经验贝叶斯);
  9. fit2 <- eBayes(fit2, 0.01)
  10. #提取TOP500基因表格;
  11. tT <- topTable(fit2, adjust="fdr", sort.by="B", number=500)
  12. #选择相关的列;
  13. tT <- subset(tT, select=c("ID","adj.P.Val","P.Value","logFC","GENE_SYMBOL","ENSEMBL_ID"))
  14. head(tT)
  15. #查看表格行列数;
  16. dim(tT)
复制代码


  1. #提取所有基因表格;
  2. tT2 <- topTable(fit2, adjust="fdr", sort.by="p", number=Inf)
  3. #保存表格;
  4. write.csv(tT2, file=paste0(cts,"_diff.csv"), row.names=F)
复制代码

3. 数据的可视化

  1. #绘制直方图;
  2. hist(tT2$adj.P.Val, col = "gold",
  3. border = "white", xlab = "P-adj",
  4. ylab = "Number of genes",
  5. main = "P-adj value distribution")
复制代码


  1. #鉴定差异表达基因;
  2. dT <- decideTests(fit2, adjust.method="fdr", lfc=1.5,p.value=0.05)
  3. #1对应"up", -1对应"down" ,0对应 "not expressed";
  4. tail(dT,8)
  5. #参看比较组名称;
  6. colnames(fit2)
  7. #选择对照组;
  8. ct <- 1
  9. #绘制火山图;
  10. volcanoplot(fit2, coef=ct, main=colnames(fit2)[ct], pch=20,
  11. highlight=length(which(dT[,ct]!=0)), names=rep('+', nrow(fit2)))
复制代码


  1. #Mean-difference plot (MD-plot);
  2. plotMD(fit2, column=ct, status=dT[,ct], legend=F, pch=20, cex=1)
  3. abline(h=0,col="blue")
复制代码


好啦,本次的分享就到这里啦!

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


本文作者:基迪奥-莫北

本帖子中包含更多资源

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

x
新的一天加油!

最近看过此主题的会员

回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2021.4.18
在线时间
23 小时

发表于 2022.10.27 09:04:23 | 显示全部楼层
看看
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.10.27 09:29:52 | 显示全部楼层
学习
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.10.27 16:07:26 | 显示全部楼层
新的一轮加油
新的一天加油!
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2021.4.18
在线时间
23 小时

发表于 2022.10.28 08:37:04 | 显示全部楼层
看看
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.10.28 08:58:50 | 显示全部楼层
打卡
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.10.28 22:08:28 | 显示全部楼层
加油第二天啊
新的一天加油!
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.10.29 08:44:45 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2021.4.18
在线时间
23 小时

发表于 2022.10.29 10:51:54 | 显示全部楼层
看看
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.10.30 09:35:39 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

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

发表于 2022.10.30 10:17:30 | 显示全部楼层
10月倒数第二天
新的一天加油!
回复 支持 反对

使用道具 举报

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

本版积分规则

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