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

news2024/11/20 1:49:50

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

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

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

无人机集群路径规划：遗传算法求解无人机集群路径规划，提供MATLAB代码

无人机集群路径规划：遗传算法求解无人机集群路径规划，提供MATLAB代码

一、单个无人机路径规划模型介绍 无人机三维路径规划是指在三维空间中为无人机规划一条合理的飞行路径&#xff0c;使其能够安全、高效地完成任务。路径规划是无人机自主飞行的关键技术之一&#xff0c;它可以通过算法和模型来确定无人机的航迹&#xff0c;以避开障碍物、优化…
阅读更多...
基于springboot的毕业设计系统的开发源码

基于springboot的毕业设计系统的开发源码

风定落花生&#xff0c;歌声逐流水&#xff0c;大家好我是风歌&#xff0c;混迹在java圈的辛苦码农。今天要和大家聊的是一款基于springboot的毕业设计系统的开发。项目源码以及部署相关请联系风歌&#xff0c;文末附上联系信息 。 项目简介&#xff1a; 毕业设计系统能够实现…
阅读更多...
微软：最新ChatGPT-4o模型，可在 Azure OpenAI上使用

微软：最新ChatGPT-4o模型，可在 Azure OpenAI上使用

北京时间5月14日凌晨&#xff0c;OpenAI 一场不到 30 分钟的发布会&#xff0c;正式发布了 GPT-4o&#xff0c;视频语音交互丝滑到吓人&#xff0c;还即将免费可用&#xff01; GPT-4o&#xff0c;其中的「o」代表「omni」&#xff08;即全面、全能的意思&#xff09;&#xff…
阅读更多...
和程序员de 相处之道

和程序员de 相处之道

1、不要"一遇到问题就去找程序员" 通常&#xff0c;技术问题通过阅读使用说明就可以解决。比如你刚买了一个新的播放器&#xff0c;想要把它连接到你的电视&#xff0c;你只需要找到使用手册里关于如何连接接口的那一页即可。 错误信息通常会被很清晰地列出来。通过…
阅读更多...
需不需要选数据结构和算法的课程？

需不需要选数据结构和算法的课程？

网络工程是互联网方向&#xff1f; 不过如果你想走机器学习和大数据方向的话&#xff0c;数据结构以及算法应该说是必修了吧。刚好我有一些资料&#xff0c;是我根据网友给的问题精心整理了一份「数据结构的资料从专业入门到高级教程」&#xff0c; 点个关注在评论区回复“88…
阅读更多...
XShell-连接-Centos 7

XShell-连接-Centos 7

XShell 连接Centos 7 一.准备 安装XShell XShell下载地址&#xff1a; 在虚拟机上安装Centos 7&#xff0c;具体操作自行学习 二.Centos 7的准备 1.网络适配器修改为NAT 2.获取IP 输入命令&#xff1a; ip addr我的Centos 7对外IP为192.168.174.129 三.XShell连接Cento…
阅读更多...
全志T527芯片详解【一】：计算性能

全志T527芯片详解【一】：计算性能

高效架构 性能稳健 内置8*ARM Cortex-A55&#xff0c;8核主频可运行至1.8GHz&#xff0c;提供更稳健强劲的处理能力 Octa-core ARM Cortex-A55 in a DynamlQ big.LITTLE configuration, up to 1.8 GHz32KB L1 I-cache and 32KB L1 D-cache per A55 coreOptional 64KB L2 cach…
阅读更多...
Linux 如何用上次的checkpoint文件dist_train.sh 接着训练【mmdetection】

Linux 如何用上次的checkpoint文件dist_train.sh 接着训练【mmdetection】

在Linux环境下&#xff0c;如果你想要用上一次的checkpoint文件继续训练&#xff0c;你可以在你的dist_train.sh脚本中设置--resume_from参数。这个参数指定了checkpoint文件的路径&#xff0c;训练会从该文件的状态继续进行。 例如&#xff0c;如果你的checkpoint文件名为las…
阅读更多...
海山数据库（He3DB）线程池方案详解

海山数据库（He3DB）线程池方案详解

前言 对于应用开发人员来说肯定听说过连接池&#xff0c;却不一定听说过线程池&#xff0c;虽然二者都是池化的概念&#xff0c;但还是有所不同的&#xff1a; 连接池面向的是数据库连接&#xff0c;是针对数据库Client侧的优化。连接池可将数据库连接数固定在一定范围内&#…
阅读更多...
Unity3D读取Excel表格写入Excel表格

Unity3D读取Excel表格写入Excel表格

系列文章目录 unity工具 文章目录 系列文章目录&#x1f449;前言&#x1f449;一、读取Excel表格&#x1f449;二、写入Excel表格&#x1f449;三、Fileinfo和Directoryinfo的操作&#x1f449;四、壁纸分享&#x1f449;总结 &#x1f449;前言 有时候难免会遇到读取文件写…
阅读更多...
基于STM32看Cotex-M内核中断向量表重定向

基于STM32看Cotex-M内核中断向量表重定向

文章目录 ARM处理器的启动流程嵌套中断向量控制器NVIC向量表和向量表重定位IAP过程中的重定向没有VTOR寄存器的情况 中断向量表重定向在芯片IAP(或OTA)升级过程中&#xff0c;是不可避免的问题。 在Cotex-M3/M4内核中有向量表偏移寄存器VTOR 在Cortex-M0内核中没有向量表偏移寄…
阅读更多...
什么是住宅IP代理？为什么需要家庭 IP 代理

什么是住宅IP代理？为什么需要家庭 IP 代理

家庭代理 IP 允许您选择特定位置&#xff08;国家、城市或移动运营商&#xff09;并作为代理上网该区域的真实用户。住宅代理 IP 可以定义为保护用户免受一般网络流量影响的中介。它们在隐藏您的 IP 地址的同时充当缓冲区。住宅代理 IP 是服务提供商分配给用户的替代 IP 地址。…
阅读更多...
Windows部署Jar包到系统服务（Service）

Windows部署Jar包到系统服务（Service）

使用WinSW工具 1、工具下载地址&#xff1a;https://github.com/winsw/winsw/releases 选择最新版本下载 根据机器32位或者64位分别下载exe&#xff0c;再下载sample-minimal.xml文件 2、修改文件名 将两个文件名称修改为服务名&#xff0c;如&#xff1a; test.exe 和 test…
阅读更多...
Nodejs+Socket.io+Web端完成聊天

Nodejs+Socket.io+Web端完成聊天

前言 源码获取:nodeexpresssocket.ioweb: 聊天demo (gitee.com) 目录结构 后端依赖 启动方式 前端是html正常启动 后端是node app.js 后端app.js核心代码 const express require(express) const app express() var http require(http).Server(app) var io require(so…
阅读更多...
Ping32和IPguard两款文件透明加密软件进行实测分享

Ping32和IPguard两款文件透明加密软件进行实测分享

透明加密软件在现代企业的数据安全保护中扮演着至关重要的角色。它们通过自动加密和解密文件&#xff0c;确保数据在存储、传输和使用过程中的安全性&#xff0c;从而防止敏感信息泄露。本文将介绍几种常见的透明加密软件&#xff0c;并对Ping32和IPguard两款加密软件进行实测分…
阅读更多...
压铸模具适合采用3D打印随形水路吗

压铸模具适合采用3D打印随形水路吗

3D打印随形水路是一种基于3D打印技术的新型模具冷却水路设计。这种设计方式可以很好地贴合产品形状&#xff0c;有效提升产品良率和冷却效率&#xff0c;3D打印随形水路最初多应用在注塑模具上&#xff0c;而随着3D打印技术的发展和新材料的不断丰富&#xff0c;压铸模具在逐步…
阅读更多...
DOS学习-目录与文件应用操作经典案例-more

DOS学习-目录与文件应用操作经典案例-more

新书上架~&#x1f447;全国包邮奥~ python实用小工具开发教程http://pythontoolsteach.com/3 欢迎关注我&#x1f446;&#xff0c;收藏下次不迷路┗|&#xff40;O′|┛ 嗷~~ 目录 一.前言 二.使用 三.案例 一.前言 DOS系统的more命令是一个用于查看文本文件内容的工具。…
阅读更多...
一个典型的分布式缓存系统是什么样的？no.32

一个典型的分布式缓存系统是什么样的？no.32

分布式 Redis 服务 由于本课程聚焦于缓存&#xff0c;接下来&#xff0c;我将以微博内的 分布式 Redis 服务系统为例&#xff0c;介绍一个典型的分布式缓存系统的组成。 微博的 Redis 服务内部也称为 RedisService。RedisService 的整体架构如图所示。主要分为Proxy、存储、集…
阅读更多...
日处理100吨污水处理设备安装需要多久

日处理100吨污水处理设备安装需要多久

日处理100吨污水处理设备的安装时间取决于多种因素&#xff0c;包括设备的复杂性、安装地点的条件、所需的基础设施建设、以及安装团队的经验和效率等。以下是一个大致的安装时间框架和相关的考虑因素&#xff1a; 前期准备&#xff1a; 现场勘查和设计&#xff1a;1-2周&#…
阅读更多...
【科研小小白】Faster R-CNN论文阅读笔记+与FAST RCNN区别对比+个人补充知识

【科研小小白】Faster R-CNN论文阅读笔记+与FAST RCNN区别对比+个人补充知识

论文阅读笔记 网络结构 整个Faster R-CNN可以分为三部分&#xff1a; **backbone&#xff1a;**共享基础卷积层&#xff0c;用于提取整张图片的特征。例如VGG16&#xff0c;或Resnet101&#xff0c;去除其中的全连接层&#xff0c;只留下卷基层&#xff0c;输出下采样后的特征…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023