数据挖掘学习笔记——GEO数据库：芯片数据分析

已于 2023-10-22 11:18:29 修改
阅读量2.9k 收藏 78
点赞数 13
分类专栏： 生物信息学 文章标签： 数据挖掘 学习 数据库
于 2023-02-13 15:22:52 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_45794091/article/details/128888235
版权

数据挖掘

数据挖掘学习笔记——GEO数据库：芯片数据分析

文章目录

虽然目前市场上高通量测序成为了主流，不过芯片测序也发展了十几年，在公共数据库中留下了大量的值得被数据挖掘的数据集，所以学习一下GEO数据挖掘也是有意义的。

一、芯片基础知识

1.1、背景

  • 意义：基因芯片分析就是为了通过生物信息学方法从这些芯片数据中发现可能对生物效应起作用的关键基因，从中寻找特定模式并对每个基因给予注释，从而挖掘出隐含的生物学过程并抽提出生物学的或功能层面上的意义。

  • 根据芯片的使用目的，一张芯片可能包含数十、数百甚至数十万的不同序列。

    • 探针：被排列成矩阵的DNA片段
    • 靶标：样本RNA则被成为靶标

  • 芯片实验过程：样本mRNA首先被反转录成cDNA（在过程中同时被荧光标记），后与芯片上的核酸探针混合，互补杂交的cDNA就结合到芯片上，而未被杂交的样本被洗脱掉。芯片被一个荧光扫描仪扫描后，芯片上某个位置探针结合上了样本中互补的核酸，就在该位置显出了一个荧光点，此位置提示基因的身份，而荧光强度则提示了原始样本中该mRNA水平的高低。

  • 芯片技术中的两种基本方法：单染色技术和双染色技术。

    • 单染色技术：将一个样本经一种荧光标记后单独杂交的一张芯片上（使用/数据集最多）。产生的芯片数据为单通道信号数据
      • 优点：方便在多张芯片之间进行比较。
      • 缺点：数据变异大，需要通过重复实验来减少误差。
    • 双染色技术：将两个样本用不同荧光标记后一起杂交到同一张芯片上。用于检测两种不同条件下基因表达的差异情况，如样品组与对照组。两个样本被两种不同荧光标记（Cy5（红色染料）标记，Cy3（绿色染料）标记）后，与芯片上的探针竞争杂交。产生双通道信号数据。
      • 优点：此数据便于两样本间的直接比较，有助于减少数据变异性，提高组间差异表达分析的准确性，同时减少了芯片的使用量，节约了成本。
      • 缺点：使用这种技术需要确定好实验设计，因而无法与其他样本进行比较。

  • 主流平台：Affymetric公司、Agilent公司和Illumina公司。

  • 数据下载：NCBI的GEO和EBI的ArrayExpress是目前最大的公开资源数据库，用于存储和发布与MIAME相容的芯片数据。

二、GEO数据库概述

2.1、基础简介

  • 简介：GENE EXPRESSION OMNIBUS（GEO）数据库，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。起先只是为表达芯片数据准备的，后期纳入了各种NGS组学数据。这些数据包括基于单通道和双通道微阵列的实验，检测mRNA，基因组DNA和蛋白质丰度，以及非阵列技术，如基因表达系列分析（SAGE），质谱蛋白质组学数据和高通量测序数据。
  • 四种基本实体类型：
    • GEO Platform (GPL) 平台
    • GEO Sample (GSM) 样本
    • GEO Series (GSE) 系列
    • GEO Dataset (GDS) 数据集
  • 前三个（样本，平台和系列）由用户提供；数据集由GEO工作人员根据用户提交的数据进行编译和策划。
  • 四者关联：一篇文章可以有一个或者多个GSE数据集，一个GSE里面可以有一个或者多个GSM样本。多个研究的GSM样本可以根据研究目的整合为一个GDS，不过GDS本身用的很少。而每个数据集都有着自己对应的芯片平台，就是GPL。

2.2、检索页面展示

  • 页面：

  • 检索：

  • 情况一：只能确定大体研究方向（如肺癌lung cancer）

  • 在NCBI官网中All Databases下拉框中选择GEO DataSets后搜索lung cancer
    在这里插入图片描述

  • 情况二：阅读文献后获得具体的GSE号（如GSE42872）

    • 方法一：浏览器中键入https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84498（通过修改GSE号检索不同GSE）
    • 方法二：GEO首页内搜索GSE号
      在这里插入图片描述
      在这里插入图片描述

​soft和miniml都是表示该platform的基础信息，比如GPL编号，上传日期等。

三、GSE项目的三种下载方式

3.1、主页下载原始数据

在这里插入图片描述

由于针对原始数据使用的不同芯片，分析处理的方法也不同，所以不推荐这种方法。

3.2、主页下载表达矩阵

  • 下载：在这里插入图片描述
    在这里插入图片描述

  • 读取：

    • 函数：read.table()
    • 示例：
gset <- read.table("GSE42872_series_matrix.txt.gz",comment.char = "!",header = T,row.names = 1)

在这里插入图片描述

3.3、GEOquery包下载表达量

利用GEOquery包读取GEF号。和第二种方法的结果是一样的，都是将压缩包下载下来。
在这里插入图片描述

  • 函数：GEOquery包的getGEO函数
  • 语法：
getGEO(
  GEO = NULL, # 输入GDS/GPL/GSE号
  destdir = tempdir(), # 指定下载地址
  GSEMatrix = TRUE, # 默认为 TRUE，表示不会进行 SOFT 格式解析
  AnnotGPL = FALSE, # 表示是否使用注释 GPL 信息
  getGPL = TRUE, # 是否下载 GPL 信息（SOFTT文件），指定F可以节省时间（文件大）
  parseCharacteristics = TRUE # 是否解析 GSE 矩阵文件的特征信息
)
  • 示例：

1、根据GDS号来下载数据，下载soft文件

gds858 <- getGEO(‘GDS858’, destdir=“.”，getGPL = F)

2、根据GPL号下载的是芯片设计的信息！

gpl96 <- getGEO(‘GPL96’, destdir=“.”，getGPL = F)

3、根据GSE号下载数据，下载series_matrix.txt.gz

gse1009 <- getGEO(‘GSE1009’, destdir=“.”，getGPL = F)
  • 下载数据并获得表达矩阵
# 加载R包，suppressPackageStartupMessages函数取消了加载R包时的输出内容
suppressPackageStartupMessages(library(GEOquery))
# 下载数据，如果文件夹中有会直接读入（这边以GSE42872为例）
gset = getGEO('GSE42872', destdir=".",getGPL = F)
# 此时的gset 还是一个对象（list），需要将里面的东西取出来才能进行下一步分析。
# gset[[1]]就是list的第一层，也就是第一个GPL，由于该GSE只有一个GPL所以也只有一层list。
gset=gset[[1]]
# 获取表达矩阵
expr <- exprs(gset)

在这里插入图片描述

四、基因名与探针ID的转换技巧

—般分析芯片数据都需要把探针的ID切换成基因名

4.1、获取对照关系

4.1.1、利用Bioconductor中汇总的R包

我将常用的GPL平台对应的BIOconductor的R包整理为了表格，直接进行导入就可以了。

# 根据GPL查询R包
platformMap <- read.csv("platformMap.csv",header = T,row.names = 1)
index <- "GPL6244" # 改GPL号
package <- platformMap$bioc_package[grep(index,platformMap$gpl)] 
package_db <- paste0(package,".db")

# 下载并调用R包
BiocManager::install(package_db)
library(hugene10sttranscriptcluster.db) # 改library包

# 导出探针ID与基因名的对照矩阵
probe2symbol_df <- toTable(get(paste0(package,"SYMBOL")))

在这里插入图片描述

4.1.2、利用平台的数据

不是所有的平台数据都被包装为了R包，剩下的这些需要我们手动绘制ID转换表格。

4.1.2.1、直接利用GEOquery包下载
  • 1、下载数据
# 下载GPL数据（也就是soft文件）
library(GEOquery)
GPL6244 <-getGEO('GPL6244',destdir =".")

# 读取对象
GPL6244_anno <- Table(GPL6244)

基因名称在了gene_assignment这一列的中间
在这里插入图片描述

  • 2、提取数据，绘制为ID对照表格
    我们需要将探针ID（第一列）和gene_assignment列中的基因名提取出来
library(dplyr) # 函数filter
library(tidyr) # 函数separate
probe2symbol_df1 <- GPL6244_anno %>% # 用管道进行连续操作
  select(ID,gene_assignment) %>%  # select函数将需要的列提取出来
  filter(gene_assignment != "---") %>% # filter函数将gene_assignment列中= "---"的行全部过滤掉
  separate(gene_assignment,c("drop","symbol"),sep="//") %>% # separate函数将gene_assignment按照//作为分隔符，保留前两项并命名
  select(-drop)# 基因名在第二列，将第一列删除，并将结果传递给probe2symbol_df

在这里插入图片描述

4.1.2.2、从NCBI中下载GPL的soft包（适用网络不好时）

在这里插入图片描述
在这里插入图片描述

在这里插入图片描述

  • 流程：
GPL6244_anno <-data.table::fread("GSE42872_family.soft.gz",skip ="ID")
# skip ="ID"不能漏

library(dplyr) # 函数filter
library(tidyr) # 函数separate
probe2symbol_df2 <- GPL6244_anno %>% # 用管道进行连续操作
  select(ID,gene_assignment) %>%  # select函数将需要的列提取出来
  filter(gene_assignment != "---") %>% # filter函数将gene_assignment列中= "---"的行全部过滤掉
  separate(gene_assignment,c("drop","symbol"),sep="//") %>% # separate函数将gene_assignment按照//作为分隔符，保留前两项并命名
  select(-drop)# 基因名在第二列，将第一列删除，并将结果传递给probe2symbol_df

在这里插入图片描述

4.1.2.3、从NCBI中下载GPL的full table（适用网络不好时）

这里读取到的表格也可以在网页上查看
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

GPL6244_anno1 <-data.table::fread("GPL6244-17930.txt")

library(dplyr) # 函数filter
library(tidyr) # 函数separate
probe2symbol_df3 <- GPL6244_anno1 %>% # 用管道进行连续操作
  select(ID,gene_assignment) %>%  # select函数将需要的列提取出来
  filter(gene_assignment != "---") %>% # filter函数将gene_assignment列中= "---"的行全部过滤掉
  separate(gene_assignment,c("drop","symbol"),sep="//") %>% # separate函数将gene_assignment按照//作为分隔符，保留前两项并命名
  select(-drop)# 基因名在第二列，将第一列删除，并将结果传递给probe2symbol_df

这一个例子中，利用R包导入的基因量少于平台数据整理出来的数据量
在这里插入图片描述

4.2、ID转换

现在我们希望将表达矩阵中的探针ID改成基因名
只需要将两个数据框按照探针名称合并就可以了， 同时在这里，我们还需要把重复的探针去掉，方法就是计算每一个探针的平均值，留下最大的那个。
在这里插入图片描述

# 调整两个数据框的探针名字和类型，方便合并
colnames(exprSet)[1] <- colnames(probe2symbol_df)[1]
exprSet <- as.data.frame(exprSet)

library(dplyr)

exprSet1 <- exprSet %>% 
  inner_join(probe2symbol_df1,by="probe_id") %>% #合并探针的信息
  select(-probe_id) %>% #去掉多余信息
  select(symbol, everything()) #重新排列，
 
#求出平均数(rowMean),并创建新的列（mutate函数） 
  a <- exprSet1[,c(grep("GSM", colnames(exprSet1)))]
for (i in 1:ncol(a)) {
  a[,i] <- as.numeric(unlist(a[,i]))
}
a <- cbind(rowMeans(a),exprSet1) 

  exprSet <- a %>%
  arrange(desc("rowMeans(a)")) %>% #把表达量的平均值按从大到小排序
  distinct(symbol,.keep_all = T) %>% # symbol留下第一个
  select(-"rowMeans(a)") %>% #反向选择去除rowMean这一列
  tibble::column_to_rownames(colnames(.)[1]) # 把第一列变成行名并删除

此时得到的表达矩阵已经可以用于后续任何需要的分析了。

五、差异化分析——limma（筛选差异基因）

一般文章中会将差异基因的数量控制在200-300的范畴（当然越少越好）
这边介绍一下最适合分析芯片数据的差异化分析包——limma

5.1、简介

  • 概述：Linear Models for Microarray Data（limma）基于线性模型建模在这里插入图片描述
  • 官方帮助文档：链接（第8章 Linear Models Overview是重点）

5.2、前期准备（input文件）

  • 表达矩阵 (exprSet)：芯片数据可以通过exprs() ，也就是我们上述步骤得到的表达矩阵，常规的转录组可以通过read.csv()等导入。
  • 分组矩阵 (design) ：将实验样本（表达矩阵的列）按照实验设计分成几组（如样本组-对照组（case - control））（实验分组与实验样本数量无关）【通过model.matrix()得到】
  • 比较矩阵 (contrast)：指定函数如何进行组间比较【通过makeContrasts()得到】
  • 示例：
# 表达矩阵（exprSet）
# BiocManager::install('CLL')
suppressPackageStartupMessages(library(CLL))
data(sCLLex) # sCLLex是CLL包中的一个对象（例子）
sCLLex
samples=sampleNames(sCLLex)
exprSet=exprs(sCLLex)# 获取表达矩阵,列是各个样本，行是每个探针ID  
pdata=pData(sCLLex)# 获取分组信息
group_list=as.character(pdata[,2])# 描述了样本的状态，需要根据它来进行样本分类

# 分组矩阵(design)
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)


# 差异比较矩阵（contrast）
# 比较矩阵的组之间用-连接
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较
# 如果分组时有三个及三个以上，进行两两对比时，比较矩阵内就需要注明对象与顺序了
# contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)

5.3、主要流程

  • lmFit()：线性拟合模型构建【需要：表达矩阵和分组矩阵】 ，得到的结果再和差异比较矩阵一起导入contrasts.fit()函数
  • eBayes()：利用上一步contrasts.fit()的结果
  • topTable()：利用上一步eBayes()的结果，并最终导出差异分析结果
# limma差异化分析
# step1，lmFit()
fit <- lmFit(exprSet,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
# step2，eBayes()
fit2 <- eBayes(fit2)
# step3，topTable()
tempOutput = topTable(fit2, coef=1, n=Inf)
# 使用coef参数，这里设为1，也就是表示根据上面比较矩阵的第一个比较信息来提取结果，由于这个例子中只有两个分组，所以只会有一个比较信息
# n=Inf：展示的最大值设为无限大，可以不需要该参数
# 参数adjust="BH"表示对p值的校正方法，包括了："none", "BH", "BY" and "holm"，默认为BH，适用与多数情况，所以这边没有修改这个参数。
nrDEG = na.omit(tempOutput) # 删除NA值
write.csv(nrDEG,"C:/Users/16748/Desktop/limma_notrend.results.csv",quote = F) # 输出差异分析矩阵
  • P值校正：p值是针对单次检验，假设采用的p值为小于0.05，我们通常认为这个基因在两组样本中的表达是有显著差异的，但是仍旧有5%的概率表示这个基因并不是差异基因。但是，当两组样本中有20000个基因采用同样的检验方式进行统计检验时，就会遇到一个问题：单次犯错的概率为0.05， 如果进行20000次检验，那么就会有0.05*20000=1000 个基因在组间的差异被错误估计。简单来说，P值校正是为了降低大样本中假阳性的数量，提升数据可信度。

5.4、结果解读

在这里插入图片描述
内容包括AveExpr值（所有样本的平均表达量）、logFC值（差异倍数）、t值、P值、q值（即adj.P.Val值）和B值。一般logFC值、P值、q值和AveExpr值用来作为差异表达的判断标准，比如差异倍数在2倍以上、平均表达在10以上、P值小于0.01等。
t值是统计量，
• 第一列：“logFC”是两组表达值之间以2为底对数化的的变化倍数（Fold change, FC），由于基因表达矩阵本身已经取了对数，这里实际上只是两组基因表达值均值之差；
• 第二列：“AveExpr”是该探针组所在所有样品中的平均表达值；
• 第三列：“t”是贝叶斯调整后的两组表达值间 T 检验中的 t 值，根据t值可以判断p值进而判断模型是否显著；
• 第四列：“P.Value”是贝叶斯检验得到的 P 值（越小越好）；
• 第五列：“adj.P.Val”是调整后的 P 值（也称为q值，越小越好）.
• 第六列：“B”是经验贝叶斯得到的标准差的对数化值。

5.5、关于比较矩阵

上述提到差异分析需要三个矩阵，表达矩阵包含样本信息，而样本必须进行分组，才能分析，所以这两个矩阵肯定需要，但比较矩阵则不一定需要。
先看下示例

  • 不需要差异比较矩阵：
library(CLL)
data(sCLLex)
library(limma)
design <- model.matrix(~factor(sCLLex$Disease))
fit <- lmFit(sCLLex,design)
fit <- eBayes(fit)
topTable(fit,coef=2,adjust='BH') 

# 结果
           logFC AveExpr      t   P.Value adj.P.Val     B
39400_at  1.0285   5.621  5.836 8.341e-06   0.03344 3.234
36131_at -0.9888   9.954 -5.772 9.668e-06   0.03344 3.117
33791_at -1.8302   6.951 -5.736 1.049e-05   0.03344 3.052
1303_at   1.3836   4.463  5.732 1.060e-05   0.03344 3.044
36122_at -0.7801   7.260 -5.141 4.206e-05   0.10619 1.935
36939_at -2.5472   6.915 -5.038 5.362e-05   0.11283 1.737
41398_at  0.5187   7.602  4.879 7.824e-05   0.11520 1.428
32599_at  0.8544   5.746  4.859 8.207e-05   0.11520 1.389
36129_at  0.9161   8.209  4.859 8.212e-05   0.11520 1.389
37636_at -1.6868   5.697 -4.804 9.355e-05   0.11811 1.282
  • 需要差异比较矩阵：
library(CLL)
data(sCLLex)
library(limma)
design <- model.matrix(~0+factor(sCLLex$Disease))
colnames(design) <- c('progres.','stable')
fit <- lmFit(sCLLex,design)
cont.matrix <- makeContrasts('stable-progres.',levels = design)
fit2 <- contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)
topTable(fit2,adjust='BH')


# 结果
           logFC AveExpr      t   P.Value adj.P.Val     B
39400_at  1.0285   5.621  5.836 8.341e-06   0.03344 3.234
36131_at -0.9888   9.954 -5.772 9.668e-06   0.03344 3.117
33791_at -1.8302   6.951 -5.736 1.049e-05   0.03344 3.052
1303_at   1.3836   4.463  5.732 1.060e-05   0.03344 3.044
36122_at -0.7801   7.260 -5.141 4.206e-05   0.10619 1.935
36939_at -2.5472   6.915 -5.038 5.362e-05   0.11283 1.737
41398_at  0.5187   7.602  4.879 7.824e-05   0.11520 1.428
32599_at  0.8544   5.746  4.859 8.207e-05   0.11520 1.389
36129_at  0.9161   8.209  4.859 8.212e-05   0.11520 1.389
37636_at -1.6868   5.697 -4.804 9.355e-05   0.11811 1.282
  • 解析：差异比较矩阵的需要与否，主要看分组矩阵如何制作的
    • design=model.matrix(~factor(sCLLex$Disease))
    • design=model.matrix(~0+factor(sCLLex$Disease))
  • 第一种方法把需要比较的组做出到了一列，需要比较多次，就有多少列，第一列是截距不需要考虑，第二列开始往后用coef这个参数可以把差异分析结果一个个提取出来。
  • 第二种方法，仅仅是分组而已，组之间需要如何比较，需要自己再制作差异比较矩阵，通过makeContrasts函数来控制如何比较。
  • 个人看法：用第二种方法（需要比较矩阵）
    • 第一种方法在算法上优于第二种，在构建矩阵时其将常数项都列为1，在运算时就可以忽略常数项（关于线性方程的算法，可以参考下这篇文章 ）
    • 但是第二种方法可以自行定义组间比较的方式，其更加方便。第一种方法由于是写死的，固定为矩阵中的“1”表示的样本减去矩阵中的“0”表示的样本，其logFC与t值可能出现符号相反值的情况（由原先的正表达变为负表达了）
      在这里插入图片描述

5.6、关于分组矩阵

差异分析的一个使用难点就是构建分组矩阵。具体使用方法在帮助文档的第九章（p42）都有讲解，这边就不过多赘述了。

六、总流程代码

# 加载R包
suppressPackageStartupMessages(library(GEOquery))
# 下载数据，如果文件夹中有会直接读入
gset = getGEO('GSE42872', destdir=".",getGPL = F)
# 获取ExpressionSet对象，包括的表达矩阵和分组信息
gset=gset[[1]]

# 获取分组信息
pdata <- pData(gset)
SampleID <- row.names(pdata) 
Type <- c(rep("Control",3),rep("Experiment",3))
pdata <- cbind(SampleID,Type)
# 获取表达矩阵
exprSet <- exprs(gset)

# QC检测
opar <- par(no.readonly = T)
par(cex = 0.7)
n.sample <- ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample)
boxplot(exprSet,main="expression value",las=2)
par(opar)

# 数据校正
library(limma) 
exprSet <- normalizeBetweenArrays(exprSet)

# 自动判断是否需要log转换
ex <- expr
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)

if (LogC) { ex[which(ex <= 0)] <- NaN
exprSet <- log2(ex)
print("log2 transform finished")}else{print("log2 transform not needed")}

platformMap <- read.csv("platformMap.csv",header = T)
index <-  gset@annotation
package <- platformMap$bioc_package[grep(index,platformMap$gpl)] 
package_db <- paste0(package,".db")

# 下载并调用R包
BiocManager::install(package_db)
library(hugene10sttranscriptcluster.db)

# 导出探针ID与基因名的对照矩阵
probe2symbol_df1 <- toTable(get(paste0(package,"SYMBOL")))
colnames(probe2symbol_df1)[1] <- "ID"

# 调整两个数据框的探针名字和类型，方便合并
exprSet <- cbind(rownames(exprSet),exprSet) 
colnames(exprSet)[1] <- colnames(probe2symbol_df1)[1]
exprSet <- as.data.frame(exprSet)
exprSet[,1] <- as.numeric(unlist(exprSet[,1]))
probe2symbol_df1[,1] <- as.numeric(unlist(probe2symbol_df1[,1]))

library(dplyr)
exprSet <- exprSet %>% 
  inner_join(probe2symbol_df1,by="ID") %>% #合并探针的信息
  select(-ID) %>% #去掉多余信息
  select(symbol, everything()) #重新排列，

#求出平均数(rowMean),并创建新的列（mutate函数） 
  a <- exprSet[,c(grep("GSM", colnames(exprSet)))]
for (i in 1:ncol(a)) {
  a[,i] <- as.numeric(unlist(a[,i]))
}
a <- cbind(rowMeans(a),exprSet) 

  exprSet <- a %>%
  arrange(desc("rowMeans(a)")) %>% #把表达量的平均值按从大到小排序
  distinct(symbol,.keep_all = T) %>% # symbol留下第一个
  select(-"rowMeans(a)") %>% #反向选择去除rowMean这一列
  tibble::column_to_rownames(colnames(.)[1]) # 把第一列变成行名并删除

# 查看内参表达量来判定ID转化是否正确
# 常见的内参有GAPDH，ACTB等(QC检测时显示了总体数据的数值水平，内参值应该高于实验组值)
exprSet["GAPDH",]

group_list=as.character(pdata[,2])# 描述了样本的状态，需要根据它来进行样本分类

# 分组矩阵(design)
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)


# 差异比较矩阵（contrast）
# 比较矩阵的组之间用-连接
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较
# 如果分组时有三个及三个以上，进行两两对比时，比较矩阵内就需要注明对象与顺序了
# contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)

# limma差异化分析
# step1，lmFit()
for (i in 1:ncol(exprSet)) {
  exprSet[,i] <- as.numeric(unlist(exprSet[,i]))
}
fit <- lmFit(exprSet,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
# step2，eBayes()
fit2 <- eBayes(fit2)
# step3，topTable()
tempOutput = topTable(fit2, coef=1, n=Inf)
# 使用coef参数，这里设为1，也就是表示根据上面比较矩阵的第一个比较信息来提取结果，由于这个例子中只有两个分组，所以只会有一个比较信息
# n=Inf：展示的最大值设为无限大，可以不需要该参数
# 参数adjust="BH"表示对p值的校正方法，包括了："none", "BH", "BY" and "holm"，默认为BH，适用与多数情况，所以这边没有修改这个参数。
nrDEG = na.omit(tempOutput) # 删除NA值
# write.csv(nrDEG,"C:/Users/16748/Desktop/limma_notrend.results.csv",quote = F) 
# 输出差异分析矩阵
福旺旺
关注
  • 13
    点赞
  • 踩
  • 78
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
零基础入门转录组数据分析——数据处理（GEO数据库——芯片数据）
weixin_49878699的博客
01-07 3464
GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据，也就是说只要是目前已经发表的论文，论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。并且这个网站作为各种高通量实验数据的公共存储库。这些数据包括基于单通道和双通道微阵列的实验，检测mRNA，基因组DNA和蛋白质丰度，以及非阵列技术，如基因表达系列分析（SAGE），质谱蛋白质组学数据和高通量测序数据。
零基础入门转录组数据分析——数据处理（GEO数据库——高通量测序数据）
weixin_49878699的博客
02-11 1万+
GEO数据库中高通量数据处理（结合了官方和自己理解），包括：基因symbol转化，获取count，fpkm处理，设置分组信息表。
GEO数据挖掘-GEO背景知识+表达芯片分析思路
最新发布
weixin_57975238的博客
05-16 2501
1.图表分析 2.GEO背景介绍及分析思路 3.代码分析流程 4.复杂数据分析 提示：这是GEO数据挖掘的大概内容： 广义的基因有 6w+个，如何缩小范围到课题相关？ 1.数据从哪里来 2.有什么类型的数据可挖掘 基因表达芯片 转录组 单细胞 突变(
生信小白菜之GEO数据库简介及数据获取方式
yangyulu1998的博客
03-21 3509
看差异基因pheatmap包输入数据是数值型matrixdata.frame颜色变化是数值大小聚类树：层次聚类。
GEO数据库学习一(简介 数据下载 芯片知识)
weixin_73362123的博客
11-27 1万+
NCBI Gene Expression Omnibus（GEO）是一个存储各种高通量实验数据的公共数据库。这些数据包括基于单通道和双通道微阵列的实验，检测mRNA，基因组DNA和蛋白质丰度，以及非阵列技术，如基因表达系列分析（SAGE），质谱蛋白质组学数据和高通量测序数据。
GEO数据挖掘（一）基础介绍
doctor_yuxiang的博客
05-07 1万+
GEO数据挖掘，火山图，热图，主成分分析
GEO数据库的使用（一）
热门推荐
xiaobai1_1的博客
12-04 3万+
1、GEO数据库介绍 GEO全称GENE EXPRESSION OMNIBUS，由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。创建于2000年，收录世界各国研究机构提交的高通量基因表达数据。 GEO上有四类数据GSM, GSE, GDS, GPL 1.GSM是单个样本的实验数据 2.GDS是人工整理好的关于某个话题的GSM的集合，一个GDS中的GSM的平台是一样的 3.GSE是一个...
GEO数据库架构介绍
庐州月光的博客
12-20 2455
欢迎关注”生信修炼手册”!GEO是一个国际化的开源项目，允许研究者提交自己的数据到该数据库，在世界范围内公开共享自己的数据，网址如下https://www.ncbi.nlm.nih.go...
GEO数据库介绍
bioprogrammer
08-10 1万+
本博客内容将同步更新到个人微信公众号：生信玩家。欢迎大家关注~~~
学习笔记(02):零基础搞定Python数据分析与挖掘-数据分析的流程
01-07
立即学习:...M:modify（数据修正，便于下一步的分析和建模） M:model（数据建模，侧重于未知事件的预测） A: assess（模型评估，检验模型的稳定性以及实用性。） 作者：卢云伟
数据挖掘学习笔记（四）
01-07
数据预处理是数据挖掘过程的第一个主要步骤，了解数据才能为分析与挖掘做好预处理。 数据的属性： 数据值类型 数据的分布 图形表示形式 数据的相似性与相异性 数据的属性 1.数据对象 数据集由数据对象组成。 一个...
数据挖掘学习笔记
12-21
数据挖掘学习笔记 现有一组facebook数据，x1-x10为已知属性，y为要预测的变量，首先导包，和数据。 import numpy as np import pandas as pd from statsmodels.formula.api import ols import statsmodels.api as ...
数据挖掘与分析学习笔记数据
03-22
数据挖掘与分析学习笔记数据
数据挖掘学习笔记（三）
01-06
数据分析与数据挖掘的方法 1.频繁模式 频繁模式：数据中频繁出现的模式。 频繁项集：频繁在事务数据集中一起出现的商品集合。 例如，信用卡分析、患者就诊分析、购物车分析… 2.分类与回归 分类与标签预测是找出描述...
GEO芯片数据分析更新（补富集分析与WGCNA）
Annaaphq的博客
12-24 2838
更新下GEO芯片数据分析
GEO数据挖掘-1 (基因芯片)
qq_44048812的博客
03-15 1209
boxplot发现异常样品怎么办(取值范围差异很大)？如何看作者提供的表达矩阵是否能直接用来差异分析？数据/代码源自生信技能树课程。富集分析(KEGG、GO)找数据，找到GSE编号。
GEO数据库学习四(差异分析 可视化 GSEA)
weixin_73362123的博客
11-29 3820
定义：在某一特定基因子集中，多个基因共同参与的pathway通路，与标准背景相比，是否比例升高、具有统计显著性。，其中N：所有gene总数n：N中差异表达gene的总数M：N中属于某个GO term的gene个数k: n中属于某个GO term的gene个数p：表示差异表达gene富集到这个GO term上的可信程度当p
GEO芯片数据基本分析
Annaaphq的博客
12-23 4428
GEO数据库芯片基本分析——R语言
使用oligo软件包处理芯片数据
R语言与生物信息学
03-27 2万+
本博客介绍过 Affy芯片的处理方法 ，其中所使用的软件包有一定的局限性，无法读取和分析一些新版Affy芯片。本文介绍oligo软件包的处理方法以解决这些问题。oligo软件包并不是新出现的软件包，只因新类型芯片的不断推出，关注它的用户越来越多。而且，除了用于Affy芯片处理外，oligo软件包还可处理NimbleGen芯片。 oligo处理芯片的原理和其他方法相同，难点在最后一步：从探针到
opencv学习笔记——cv::mean()函数详解
06-12
cv::mean()函数是OpenCV中用于计算图像均值的函数，它的函数原型如下： ``` Scalar cv::mean(InputArray src, InputArray mask = noArray()); ``` 其中，src参数是输入的图像，mask参数是用于掩码的可选参数。 函数返回值是一个Scalar类型的值，包含了计算得到的均值。 使用cv::mean()函数的步骤如下： 1. 加载图像 ``` Mat img = imread("image.jpg", cv::IMREAD_GRAYSCALE); ``` 2. 计算均值 ``` Scalar mean = cv::mean(img); ``` 在这个例子中，我们计算了灰度图像的均值，并将结果存储在名为mean的Scalar变量中。 如果想要使用掩码来计算均值，可以使用以下代码： ``` Mat mask = imread("mask.jpg", cv::IMREAD_GRAYSCALE); Scalar mean = cv::mean(img, mask); ``` 在这个例子中，我们使用名为mask的掩码来计算灰度图像的均值。 需要注意的是，如果输入图像是多通道的，则返回的Scalar类型变量包含每个通道的均值。例如，如果输入图像是RGB图像，则返回的Scalar类型变量包含红、绿、蓝三个通道的均值。 希望这个解释能够帮助你理解cv::mean()函数的使用方法。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值