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

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

文章目录

  • 探针注释
  • 自主注释流程(了解)
  • PCA图、top1000基因热图
    • 探针注释
    • 查看示例代码
  • top 1000 sd 热图
    • 离散基因热图，top1000表达基因，只是看一下，不用放文章里
  • 差异分析
  • 火山图
  • 差异基因热图
  • 转换id
  • 富集分析-KEGG数据库
  • 补充
    • 两个文件夹里分析出来的差异基因如何改交集

探针注释

在这里插入图片描述

自主注释流程(了解)

在这里插入图片描述

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的阈值

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：/a/637182.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

安装mpi4py与dlio_profiler_py的总结

安装mpi4py与dlio_profiler_py的总结

安装mpi4py mpi4py是一个Python库&#xff0c;它提供了与MPI&#xff08;Message Passing Interface&#xff09;兼容的接口&#xff0c;使得Python程序能够利用MPI实现并行计算。mpi4py 的核心是基于MPI标准的C/C实现&#xff0c;它能够在高性能计算环境下进行高效的并行处理…
阅读更多...
网页版收银系统比安装板收银系统的四大优势

网页版收银系统比安装板收银系统的四大优势

在当今竞争激烈的零售市场中&#xff0c;高效的收银系统对于连锁实体店的管理至关重要。随着科技的不断发展&#xff0c;网页版收银系统成为越来越多零售企业的首选。网页版收银系统以其灵活性、可定制性和便利性&#xff0c;成为现代零售业的利器。本文将探讨网页版收银系统相…
阅读更多...
pycharm 关闭项目卡死

pycharm 关闭项目卡死

PyCharm2023.3.4 关闭一直卡在 closing projects 解决办法&#xff1a; 打开PyCharm&#xff0c; 选择 Help -> Find Action -> 输入 Registry -> 禁用ide.await.scope.completion
阅读更多...
MYSQL 集群

MYSQL 集群

1.集群目的:负载均衡 解决高并发 高可用HA 服务可用性 远程灾备 数据有效性 类型:M M-S M-S-S M-M M-M-S-S 原理:在主库把数据更改(DDL DML DCL&#xff09;记录到二进制日志中。 备库I/O线程将主库上的日志复制到自己的中继日志中。 备库SQL线程读取中继日志…
阅读更多...
51cto已购买的视频怎么下载到电脑上？

51cto已购买的视频怎么下载到电脑上？

在数字学习的浪潮中&#xff0c;51CTO已成为众多专业人士和爱好者的知识宝库。但购买了视频课程后&#xff0c;如何将其下载到电脑上以便离线学习呢&#xff1f;这不仅是技术问题&#xff0c;更是时间管理和学习效率的关键。本文将为您揭示简单而高效的步骤&#xff0c;无论您使…
阅读更多...
前端面试项目细节重难点（已工作｜做分享）

前端面试项目细节重难点（已工作｜做分享）

面试官提问&#xff1a;需求场景&#xff1a;页面上有一个单选框&#xff0c;有是否两个选项&#xff1a;当用户选择是&#xff0c;出现一个输入框&#xff0c;用户可以输入内容&#xff0c;给后端的保存接口传入参数radio和content这两个字段&#xff0c;值分别是用户选项和输…
阅读更多...
西门子WINCC8.0VBS脚本学习讲解

西门子WINCC8.0VBS脚本学习讲解

WinCC VBS脚本置位/复位/取反 二进制变量 "TAG1_BOOL1" 进行置位复位取反操作 步骤&#xff1a;按钮-->对象属性-->事件-->单击鼠标VBS动作填入代码如下: 对二进制变量进行复位 对二进制变量进行置位 对二进制变量进行取反 VBS脚本数学运算/读写批处理 …
阅读更多...
百度智能云参与信通院多项边缘计算标准编制，「大模型时代下云边端协同 AI 发展研讨会」成功召开

百度智能云参与信通院多项边缘计算标准编制，「大模型时代下云边端协同 AI 发展研讨会」成功召开

1 中国信通院联合业界制定、发布多项标准化成果&#xff0c;推动产业发展 大模型开启了 AI 原生时代&#xff0c;云边端协同 AI 构建了「集中式大规模训练」、「边缘分布式协同推理」新范式&#xff0c;有效降低推理时延和成本&#xff0c;提升数据安全和隐私性&#xff0c;也…
阅读更多...
安卓App封装全攻略：利用小猪APP分发提升应用发布效率

安卓App封装全攻略：利用小猪APP分发提升应用发布效率

在快速迭代的移动应用市场&#xff0c;高效且安全地分发安卓应用程序是开发者面临的一大挑战。安卓App封装技术&#xff0c;作为这一挑战的解决方案之一&#xff0c;不仅能够提升应用的安全性&#xff0c;还能简化分发流程。本文将深入探讨安卓App封装的核心概念&#xff0c;以…
阅读更多...
小型发电机不发电原因和解决方法

小型发电机不发电原因和解决方法

小型发电机不发电可能由多种原因造成&#xff0c;以下是一些常见原因及其解决方法&#xff1a; 1.电池电量不足&#xff1a;小型发电机通常需要电池来启动。如果电池电量不足&#xff0c;可能导致发电机无法启动。此时&#xff0c;您可以使用充电设备对电池进行充电&#xff0…
阅读更多...
Ue目标身上放准星

Ue目标身上放准星

阅读更多...
2025秋招深度学习基础面试题（一）

2025秋招深度学习基础面试题（一）

01. 卷积和BN如何融合提升推理速度 Conv和BN的融合：在网络的推理阶段，可以将BN层的运算融合到Conv层中，减少运算量，加速推理。本质上是修改了卷积核的参数，在不增加Conv层计算量的同时，略去了BN层的计算量。def fuse_conv_bn(conv, bn):std = (bn.running_var + bn.eps).…
阅读更多...
深入解读TuGraph计算引擎模型推理系统

深入解读TuGraph计算引擎模型推理系统

作者&#xff1a;李文凯 TuGraph计算引擎模型推理系统将基于迭代计算的图计算框架与模型推理系统相结合&#xff0c;推理系统可自定义推理依赖环境&#xff0c;图迭代计算与推理链路实现隔离。基于共享内存的跨进程通信方式&#xff0c;提高了推理数据交换效率&#xff0c;满足…
阅读更多...
VUE3好看的酒网站模板源码

VUE3好看的酒网站模板源码

文章目录 1.设计来源1.1 首页界面1.2 十大名酒界面1.3 名酒新闻界面1.4 联系我们界面1.5 在线留言界面 2.效果和结构2.1 动态效果2.2 代码结构 3.VUE框架系列源码4.源码下载 作者&#xff1a;xcLeigh 文章地址&#xff1a;https://blog.csdn.net/weixin_43151418/article/detai…
阅读更多...
2024/5/22 学习杂记

2024/5/22 学习杂记

为什么功率放大电路在模电中经常提到&#xff1f; 模拟信号&#xff1a;它是连续变化的电信号&#xff0c;它在时间上和幅度上都是连续的&#xff0c;能够代表信息的连续变化。大多数物理量为模拟信号&#xff0c;如&#xff1a;温度、压力、流量… 非电物理量通过传感器变换成…
阅读更多...
Nginx - 健康检查终极指南：探索Upstream Check模块

Nginx - 健康检查终极指南：探索Upstream Check模块

文章目录 概述upstream_check_module模块安装和配置指南模块安装步骤基本配置示例详细配置说明检查类型和参数常见问题及解决方案 SSL检查和DNS解析功能SSL检查配置示例和说明配置示例 DNS解析配置示例和说明配置示例 结合实际应用场景的高级配置示例综合SSL检查与DNS解析 总结…
阅读更多...
Discourse 中可能使用的 HMAC 算法 Java 实现

Discourse 中可能使用的 HMAC 算法 Java 实现

在 DiscourseConnect 中&#xff0c;对数据的签名使用的是 HMAC 算法。 实际使用的算法为 HmacSHA256。 Java 生成签名的方法很简单。 String hmac new HmacUtils(HmacAlgorithms.HMAC_SHA_256, "55619458534897682511405307018226").hmacHex(ssoPayload);HmacUti…
阅读更多...
robosuite导入自定义机器人

robosuite导入自定义机器人

目录 目的&#xff1a;案例一&#xff1a;成果展示具体步骤&#xff1a;URDF文件准备xml文件生成xml修改机器人构建 目的&#xff1a; 实现其他标准/非标准机器人的构建 案例一&#xff1a; 成果展示 添加机器人JAKA ZU 7 这个模型 具体步骤&#xff1a; URDF文件准备 从…
阅读更多...
在深度学习中常见的初始化操作

在深度学习中常见的初始化操作

目录 截断正态分布来初始化张量 逐行代码解释 相关理论解释 截断正态分布函数 截断正态分布的定义 截断正态分布的作用 计算截断点的作用 具体步骤 正态分布的累积分布函数&#xff08;CDF&#xff09; 正态分布的累积分布函数与误差函数的关系 示例计算 误差函数 应…
阅读更多...
切换分支报错：Untracked Files Prevent Checkout

切换分支报错：Untracked Files Prevent Checkout

切换分支报错&#xff1a;Untracked Files Prevent Checkout 分支切换 Untracked Files Prevent Checkout 新起的项目在切换master分支到工作分支时&#xff0c;出现下图的问题&#xff1a; Untracked Files Prevent Checkout Move or commit them before checkout 网上的解决…
阅读更多...
最新文章

Copyright @ 2022~2023