GEO数据挖掘-PCA、差异分析

已于 2024-05-18 23:55:08 修改
阅读量1k 收藏 16
点赞数 14
分类专栏： GEO数据挖掘 文章标签： 数据挖掘 windows GEO
于 2024-05-17 23:29:48 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_57975238/article/details/139012671
版权

From 生物技能树 GEO数据挖掘第二节

文章目录

探针注释

在这里插入图片描述

自主注释流程(了解)

在这里插入图片描述

rm(list = ls())  
load(file = "step2output.Rdata")
#输入数据：exp和Group

PCA图、top1000基因热图

探针注释

查看示例代码

添加链接描述

示例代码：
# The variable Species (index = 5) is removed
# before PCA analysis
iris.pca <- PCA(iris[,-5], graph = FALSE)

fviz_pca_ind(iris.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = iris$Species, # color by groups
             palette = c("#00AFBB", "#E7B800", "#FC4E07"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
             )

dat=as.data.frame(t(exp))
#根据sthda上的示例数据，更改自己的数据，需要转置后转为dataframe

在这里插入图片描述

#画PCA图，使用pca分析后的数据dat.pca
如果每组样本在四个以下，是不会有圈的，圈是置信区间，在统计学里，小于4个样本没法计算置信区间

top 1000 sd 热图

离散基因热图，top1000表达基因，只是看一下，不用放文章里

g = names(tail(sort(apply(exp,1,sd)),1000)) #day7-apply的思考题,apply(exp,1,sd)对exp的每一行，就是一个基因，求sd，sort排序，小-到排序后，取后1000个，再提取基因名字（向量名字）
n = exp[g,]
library(pheatmap)
annotation_col = data.frame(row.names = colnames(n),
                            Group = Group)
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row", #对数据进行转换，按行标准化，只保留行内差别，不保留行间差别，会把数据范围缩放到大概-5~5之间
         breaks = seq(-3,3,length.out = 100) #设置色带分布范围为-3~3之间，超出此范围的数字显示极限颜色
         ) 
# ?pheatmap，查看帮助文档

在这里插入图片描述

差异分析

1.得到差异基因结果表格，2.做了火山图，3.做热土图
在这里插入图片描述

rm(list = ls()) 
load(file = "step2output.Rdata")
#load上一步做完后得到用于差异分析的结果表格
#需要输入数据：exp、ids、group，使用limma包根据贝叶斯检验原理进行差异分析
library(limma)
design = model.matrix(~Group)
#model.matrix()根据分组信息生成一个模型矩阵,线性拟合函数需要用模型矩阵
fit = lmFit(exp,design)
#线性拟合函数，当你执行 fit = lmFit(exp, design) 这行代码时，你正在尝试拟合一个线性模型到表达数据 exp，其中：
exp：这是一个矩阵或数据框，包含了基因表达水平的测量值。每一行通常代表一个基因，每一列代表一个样本或实验条件。
design：这是一个设计矩阵，用于指定模型中每个样本的实验条件。它通常是一个因子向量或指示变量矩阵，用于定义模型中的各项，比如不同的处理组或时间点。
fit = eBayes(fit)
#贝叶斯检验
deg = topTable(fit,coef = 2,number = Inf)
#提取差异分析结果：coef = 2指design的第二列，number = infinity指提取全部的差异基因

得到一个probe_id和对应的logFC、pValue的表格，但还需要和symbol还有entrezid连接在一起

代码如下：
1.加probe_id列，把行名变成一列，防止行名丢失

library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#mutate新增一列，probe_id = rownames(deg)：这里创建了一个新的列 probe_id，其值是 deg 的行名。rownames() 函数用于获取数据框的行名，注意不能直接添加

2.加上探针注释
probe_id和基因symbol不是一对一的关系，是多对一或者一对多的关系，因此需要去重

  1. 一个探针对应多个基因：非特异性探针，在表格中可以看到，则需要直接去除
  2. 一个基因对应多个探针：相当于一个基因测了好多遍
    处理方法：随机去重/保留行和或行平均值最大的探针/取多个探针的平均值
ids = distinct(ids,symbol,.keep_all = T)
#此处是随机去重的方法，其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
#用inner_join取交集并把差异分析的结果deg和之前的id—symbol表格连接在一起
nrow(deg) 
## [1] 20824
#检查一下，如果行数为0就是你找的探针注释是错的。

#保留最大值
exp2 = exp[ids$probe_id,]
identical(ids$probe_id,rownames(exp2))
library(dplyr)
ids = ids %>% 
  mutate(exprowsum = rowSums(exp2)) %>% 
  arrange(desc(exprowsum)) %>% 
  select(-3) %>% 
  distinct(symbol,.keep_all = T)
nrow(ids)
# 拿这个ids去inner_join

#求平均值
rm(list = ls())
load("step2output.Rdata")
exp3 = exp[ids$probe_id,]
rownames(exp3) = ids$symbol
exp3[1:4,1:4]
exp4 = limma::avereps(exp3)

# 此时拿到的exp4已经是一个基因为行名的表达矩阵，直接差异分析，不再需要inner_join

3.加change列,标记上下调基因

logFC_t = 1
p_t = 0.05
#设置logFC和pValue的阈值
#使用ifelse两次或者casewhen判断down、up、stable并输出成新的一列change
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
## 
##   down stable     up 
##    579  19621    624

#思考：如何使用padj而非p值
把代码里的所有的P.Value替换成adj.P.Val
> colnames(deg)
[1] "logFC"     "AveExpr"   "t"         "P.Value"   "adj.P.Val" "B"        
[7] "probe_id"  "symbol"

火山图

library(ggplot2)
ggplot(data = deg, aes(x = logFC, y = -log10(P.Value))) +
  geom_point(alpha=0.4, size=3.5, aes(color=change)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +
  geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
  theme_bw()

在这里插入图片描述

#使用ggplot的geom_point画火山图，vline和hline画阈值的线

差异基因热图

# 表达矩阵行名替换为基因名，分为两步：
exp = exp[deg$probe_id,]
#按deg中的symbol列的内容在exp中按行取子集,把最终使用的探针取出来
rownames(exp) = deg$symbol
#把exp中行名改为deg的symbol列，此时已经是一一对应的,上述俩代码只能运行一次,运行一次直接就把探针表达矩阵转换为基因表达矩阵了
diff_gene = deg$symbol[deg$change !="stable"]
#取出change列不是stable的基因symbol
n = exp[diff_gene,]
#按有差异的基因symbol在exp中按行取子集，即为有差异的基因的logFC和pValue，赋值到数据框n中，用于画差异基因热图

library(pheatmap)
annotation_col = data.frame(group = Group)
rownames(annotation_col) = colnames(n) 
pheatmap(n,show_colnames =F,
         show_rownames = F,
         scale = "row",
         #cluster_cols = F, 
         #即不按照列聚类，此时按照表达矩阵的顺序聚类
         annotation_col=annotation_col,
         breaks = seq(-3,3,length.out = 100)
)

在这里插入图片描述
#如果差异基因的聚类分组还是错乱的，则加cluster_col = F
#如果加了还是错乱的，去看：小洁老师的语雀/分组聚类的热图
添加链接描述

#如果行名比较少，例如改成10行，就显示出基因
pheatmap(n[1:10,],show_colnames =F,
         #show_rownames = F,
         scale = "row",
         #cluster_cols = F, 
         annotation_col=annotation_col,
         breaks = seq(-3,3,length.out = 100)
)

在这里插入图片描述

转换行名的快捷函数：探针矩阵如何转换为基因矩阵

> library(tinyarray)
> exp[1:4,1:4]
          GSM175766 GSM175767 GSM175768 GSM175769
1007_s_at  8.045017  8.314098  8.342717  8.261483
1053_at    6.444243  6.330321  6.168972  6.422393
117_at     6.158540  5.805438  5.565754  6.082891
121_at     7.737116  7.640965  7.835118  7.631916
> exp2 = trans_array(exp,ids)
20824 rownames transformed after duplicate rows removed
> exp2[1:4,1:4]
      GSM175766 GSM175767 GSM175768 GSM175769
DDR1   8.045017  8.314098  8.342717  8.261483
RFC2   6.444243  6.330321  6.168972  6.422393
HSPA6  6.158540  5.805438  5.565754  6.082891
PAX8   7.737116  7.640965  7.835118  7.631916

转换id

symbol:常说的基因名
entrezid：富集分析指定用
两个并非一一对应，损失或增加部分基因属于正常，两者可以转换
在这里插入图片描述

加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
entrezid是富集分析指定用的，需要symbol转entrezid，然后inner_join
使用clusterProfiler中的bitr函数实现，另外数据库根据物种不同

library(clusterProfiler)
library(org.Hs.eg.db)
s2e = bitr(deg$symbol, 
           fromType = "SYMBOL",
           toType = "ENTREZID",
           OrgDb = org.Hs.eg.db)#人类,注意物种，不同物种R包不同，如果物种写错，**也不会报错**，所以要检查代码错了没

一部分基因没匹配上是正常的。<30%的失败都没事。

其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb nrow(deg)
添加链接描述

看剩下数量，如果只有几十说明有问题

deg = inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
#把差异基因的表和entrezid通过inner_join连接在一起，用于后面的富集分析

多了几行少了几行都正常，SYMBOL与ENTREZID不是一对一的

nrow(deg)  #检查
## [1] 20827
#再看看还有几行，然后保存
save(exp,Group,deg,logFC_t,p_t,file = "step4output.Rdata")

富集分析-KEGG数据库

KEGG(Kyoto Encyclopedia of Genes and Genomes)是系统分析基因功能、基因组 信息数据库，它有助于研究者把基因及表达信息作为一个整体网络进行研究，以“理解生物系统的高级功能和实用程序资源库”著称。

补充

两个文件夹里分析出来的差异基因如何改交集

# 首先保存第一个project里的差异基因
exp = exp[deg$probe_id,]
rownames(exp) = deg$symbol
diff_gene = deg$symbol[deg$change !="stable"]
diff_gene2 = diff_gene
save(diff_gene2,file = "lianxi_diff_gene.Rdata")

#取交集
exp = exp[deg$probe_id,]
rownames(exp) = deg$symbol
diff_gene = deg$symbol[deg$change !="stable"]
load("../pipeline/GEO自己练习-GPL6887/lianxi_diff_gene.Rdata")
intersect(diff_gene,diff_gene2)


boxplot(exp)  #查看exp的阈值
遗落凡尘的萤火-生信小白
关注
  • 14
    点赞
  • 踩
  • 16
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
数据分析与挖掘之实验 PCA主成分提取
06-28
现有数据‘Country-data.xlsx’。数据文件包含10列，其中第1列为【国家名称】，其余九列我们暂不考虑其具体含义，仅将其视为数据标签列X1~X9，其中数据均为数字类型。现要求提取该数据集X1~X9的主成分，并设置其累计贡献率达90%，并输出其特征向量和贡献率属性。 现有数据‘Country-data.xlsx’。数据文件包含10列，其中第1列为【国家名称】，其余九列我们暂不考虑其具体含义，仅将其视为数据标签列X1~X9，其中数据均为数字类型。现要求提取该数据集X1~X9的主成分，并设置其累计贡献率达90%，并输出其特征向量和贡献率属性。
GEO数据挖掘-GEO背景知识+表达芯片分析思路
weixin_57975238的博客
05-16 1368
1.图表分析 2.GEO背景介绍及分析思路 3.代码分析流程 4.复杂数据分析 提示：这是GEO数据挖掘的大概内容： 广义的基因有 6w+个，如何缩小范围到课题相关？ 1.数据从哪里来 2.有什么类型的数据可挖掘 基因表达芯片 转录组 单细胞 突变(
学习笔记Day8：GEO数据挖掘-基因表达芯片
m0_71952366的博客
03-15 1515
GEO数据库GEO工具：GEO2R，可以导入代码进入R修改。Series：用户提交给数据库一个完整的研究，包括其样本数据(GSM)，包含提供研究描述，包括对数据描述并总结分析(GSE)。GEO数据集筛选表达芯片数据：Expression profiling by array单细胞/普通转录组（高通量测序）：Expression profiling by high throughput sequencingGSE界面：GPL（平台）中看ID和Symbol Gene；
学习笔记Day9&10：GEO数据挖掘-基因表达芯片代码
m0_71952366的博客
03-16 1026
不要在同一个文件夹中处理两套数据。
GEO数据挖掘（一）基础介绍
doctor_yuxiang的博客
05-07 1万+
GEO数据挖掘，火山图，热图，主成分分析
典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集
悟道西方
03-24 7150
典型医学设计实验GEO数据分析 (step-by-step) - 数据获取到标准化介绍了实验的设计、数据获取、数据标准化和注释，下面是如何利用Limma和线性模型鉴定差异基因，并进行GO富集分析。 线性模型 为了分析发炎和未发炎组织的差异表达，我们需要构建一个线性模型。线性模式是实验数据分析的常用方法，适用于几乎任何复杂的实验设计。后面我们专门出文介绍，推荐Mike Love和Michael Ir...
GEO数据挖掘-从数据下载到富集分析-各种绘图（学习笔记）
热门推荐
bioprogrammer
01-16 4万+
############################### ##---------GEO-data1---------## ############################### rm(list = ls()) source(&quot;http://bioconductor.org/biocLite.R&quot;) #before install package biocLite('DOSE') #...
GEO数据挖掘-DAY1
urologist的博客
03-18 793
总而言之，使用降维的算法，自定义坐标轴——把多个指标转化为少数几个综合指标（即主成分），根据这些主成分进行样本聚类，进而实现样本组内/组间的差异可视化。boxplot(exp,las = 2) #看是否有异常样本（跟其他样本截然不同），还要看纵坐标有无log （纵坐标范围在2-4之间大概率是取重复了）view(pd) #一般看title这一列（也可能是在别的列），哪些是实验组，哪些是对照组，并找到提取的关键词。#(4)提取芯片平台编号，后面要根据它来找探针注释。
GEO数据挖掘全流程分析
Eric_blog的博客
02-29 2万+
声明：以下学习资料根据“生信技能树”网络系列免费教学材料整理而成，代码来自“生信技能树”校长jimmy的github。GEO数据库挖掘系列知识分享课程，于2016年首发于生信菜鸟团博客。配套教学视频在B站，特此声明。 前言：关于GEO数据 我们的目标是要从读懂文献到复刻文献实验，再到掌握GEO数据挖掘的能力。首先便是要广泛阅读，在读文献时，提炼脉络，读懂文献使用了哪个或哪些GSE数据集，对数据做了...
数据处理-主成分分析-数据降维-pca聚类-pca
09-04
主成分分析应用 1用主成分分析法求2018年1-12月小河闸断面主要污染物（主成分）和对应的主成分得分，下面列出详细计算过程。 2pca实现人脸检测 代码思路：先用PCA实现对人脸特征的提取并降维，再使用感知机实现多...
FaceRecognitionUsing-PCA-2D-PCA-And-2D-Square-PCA:用于识别人脸的Python中PCA2D-PCA2D（Square）-PCA的实现
05-08
使用PCA-2D-PCA和2D-Square-PCA进行人脸识别用于识别人脸的Python中PCA / 2D-PCA / 2D（Square）-PCA的实现： 单人图像集团形象识别视频中的人脸ORL数据集的准确性PCA（93.42％） 二维PCA（96.05％） 2D（平方）-PCA...
《数据挖掘与大数据分析》实验报告-数据预处理
04-21
本实验报告使用数据集选自机器学习存储库 UCI，数据集标题为心脏病数据库，数据采集自...实验报告内容包括数据清洗（缺失处理、异常值、噪声）、数据归约（维度约减（特征选择、PCA）、样本抽样）、各种距离计算
数据降维-PCA算法详解
04-13
pca降维
GEO数据挖掘
weixin_40215208的博客
04-21 1164
输入数据输入数据是**数值型矩阵/数据框**颜色的变化表示数值的大小实验目的：实验目的：通过基因表达量数据的差异分析和富集分析来解释生物学现象组别设计流程：不同处理 – 差异基因 – 找功能、找关联、找机制、同一通路有多个基因表达变化 – 解释差异，缩小基因范围流程：基因表达芯片的原理探针的表达量代表基因的表达量分析思路。
数据挖掘--引论
Frankabcdefgh的博客
06-05 729
数据挖掘是从大量数据中挖掘有趣模式和知识的过程。数据源包括数据库、数据仓库、Web、其他信息存储库或动态地流入系统的数据。数据清理（消除噪声和删除不一致数据）数据集成（多种数据源可以组合在一起）数据选择（从数据库中提取与分析任务相关的数据）数据变化（通过汇总或聚焦操作，把数据变换和统一成适合挖掘的形式）数据挖掘（基本步骤，使用智能方法提取数据模式）模式评估（根据某种兴趣度量度，识别知识的真正有趣模式）知识表示（使用可视化和知识表示技术，向用户提供挖掘的知识）
【机器学习&数据挖掘】基于ARIMA 自回归积分滑动平均模型的销售价格&库存分析报告 附完整python代码
Obliviate的博客
06-04 1550
​同时销售量后1000的sku品类占比中（不畅销产品）如上，精品类产品占比第一，达到66.7%，其次是香化类产品，占比11.90%，远远小于精品类产品，酒水类产品占比7.3%，有税商品免税其他商品和电子类产品分别占比6.40%、6.40%、1.3%，将数据按照毛利进行排序，毛利前1000和后1000的sku品类占比如下，
数据挖掘--数据仓库与联机分析处理
Frankabcdefgh的博客
06-07 350
（面集时非）面向主题的：围绕某一主题来构建集成的：图片文字杂糅在一起时变的：随时间变化的数据非易失的：硬盘存放，不易丢失。
数据挖掘--挖掘频繁模式、关联和相关性：基本概念和方法
最新发布
Frankabcdefgh的博客
06-07 517
闭频繁项集：一个集合他的超集(包含这个集合的集合)在数据库里面的数量和这个集合在这个数据库里面的数量不一样,这个集合就是闭项集。，完成这一操作后，将它和题目中给出的最小支持度作比较，从而得到频繁1象集！拿到候选数据后，我们先筛选出候选频繁1象集，并算出它们的。候选3项集→频繁3象集的过程。后续同理，我们可以得到所有关系的置信度，最后我们再拿。如果这个集合还是频繁的,那么他就是极大频繁项集。频繁项集：出现的次数超过最小支持度计数阈值。1正相关，一个出现另一个就出现。
gpio-pca953x.c驱动代码分析
03-27
gpio-pca953x.c是一个GPIO驱动，用于控制PCA953x系列的GPIO扩展器。在Linux内核中，GPIO驱动的实现通常是通过向内核注册一个gpio_chip结构体，并实现一些gpio_chip操作函数来实现。而gpio-pca953x.c就是通过实现gpio_chip操作函数来控制PCA953x系列的GPIO扩展器。 下面分析gpio-pca953x.c中的一些重要函数： 1. pca953x_gpio_direction_input() 该函数用于设置GPIO方向为输入。 ```c static int pca953x_gpio_direction_input(struct gpio_chip *chip, unsigned offset) { struct pca953x_chip *pc = gpiochip_get_data(chip); int ret; mutex_lock(&pc->lock); pca953x_gpio_config(pc, offset, PCA953X_DIR_INPUT); ret = pca953x_gpio_to_irq(chip, offset); mutex_unlock(&pc->lock); return ret; } ``` 该函数首先通过gpiochip_get_data()获取pca953x_chip结构体，然后调用pca953x_gpio_config()函数设置GPIO的方向为输入。最后调用pca953x_gpio_to_irq()函数获取GPIO的中断号。 2. pca953x_gpio_direction_output() 该函数用于设置GPIO方向为输出。 ```c static int pca953x_gpio_direction_output(struct gpio_chip *chip, unsigned offset, int value) { struct pca953x_chip *pc = gpiochip_get_data(chip); int ret; mutex_lock(&pc->lock); pca953x_gpio_config(pc, offset, PCA953X_DIR_OUTPUT); pca953x_gpio_set_value(chip, offset, value); ret = pca953x_gpio_to_irq(chip, offset); mutex_unlock(&pc->lock); return ret; } ``` 该函数首先通过gpiochip_get_data()获取pca953x_chip结构体，然后调用pca953x_gpio_config()函数设置GPIO的方向为输出，并调用pca953x_gpio_set_value()函数设置GPIO的值。最后调用pca953x_gpio_to_irq()函数获取GPIO的中断号。 3. pca953x_gpio_get() 该函数用于获取GPIO的值。 ```c static int pca953x_gpio_get(struct gpio_chip *chip, unsigned offset) { struct pca953x_chip *pc = gpiochip_get_data(chip); int ret; mutex_lock(&pc->lock); ret = (pca953x_gpio_read(pc) & (1 << offset)) ? 1 : 0; mutex_unlock(&pc->lock); return ret; } ``` 该函数首先通过gpiochip_get_data()获取pca953x_chip结构体，然后调用pca953x_gpio_read()函数读取PCA953x寄存器的值，并根据GPIO的偏移量计算出GPIO的值。 4. pca953x_gpio_set() 该函数用于设置GPIO的值。 ```c static void pca953x_gpio_set(struct gpio_chip *chip, unsigned offset, int value) { struct pca953x_chip *pc = gpiochip_get_data(chip); mutex_lock(&pc->lock); pca953x_gpio_set_value(chip, offset, value); mutex_unlock(&pc->lock); } ``` 该函数首先通过gpiochip_get_data()获取pca953x_chip结构体，然后调用pca953x_gpio_set_value()函数设置PCA953x寄存器的值。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值