Learn R GEO

2023-04-02 00:54| 来源: 网络整理| 查看: 265

主要内容

•画图通用，仿制数据的思维通用，富集分析基本通用

•GEO数据库的背景知识

•GEO表达芯片的原理

•GEO表达芯片特有的下载方式

•表达芯片的差异分析（就几句代码）

•表达芯片的复杂分析

•主要学思维和方法，后面重点学习转录组的具体分析代码

图表介绍1.图表介绍1.热图·输入数据是数值型矩阵/数据框；·颜色变化表示数值大小；·热图上面横横竖竖是聚类树，为了展示数值的变化方向；·图例，根据输入的数值大小范围自动生成的颜色变化关系·相关性热图只有一半具有意义，画一半就好，但是专门的R包·差异基因热图纵坐标是样本普通热图相关性热图差异表达基因热图2.散点图3.箱线图比较组间的大小关系，以分组为单位

·输入数据是一个连续型向量和一个有重复值的离散型向量—横坐标；

·上下五条线的意思中间的又黑又粗的—中位数；上下两条线是最大值和最小值；方框的上下两条线是75%和25%（四分位数）；在外面的点-离群点

箱线图的五条线单个基因在两组之间的表达量差异散点图、箱线图4.火山图·根据logFC(横坐标)和 P value(纵坐标)可以画火山图多基因差异分析·Foldchange(FC): 处理组平均值/对照组平均值·logFoldchange(FC): Foldchange取值log2 上面标中的7.24实际上真正的表达量为2的7.24次方，是已经取过log2的数前n个样本想加除以n，后n个样本想加除以，相减(一定是处理组-对照组)log2FC·logFC>0,treat>control,基因表达量上升；上调基因：logFC>1,p解释差异，缩小基因范围数据库介绍GEO GEO网页工具GEO2R 给代码需修改GEO网页数据库介绍基因表达芯片的原理，探针的表达量代表基因的表达量，不是基因本身的表达量，所以需要将探针id转换为样本基因，他们之间存在关系，需要分组信息基因表达芯片原理分析思路代码分析流程#数据下载>rm(list = ls()) >library(GEOquery) #先去网页确定是否是表达芯片数据，不是的话不能用本流程。 >gse_number = "GSE56649" >eSet class(eSet) [1] "list" > length(eSet) [1] 1 > eSet = eSet[[1]] > class(eSet) [1] "ExpressionSet" #哪种数据类型 attr(,"package") # 属性之包的名字 [1] "Biobase" #Biobase包规定的格式"ExpressionSet"(1)提取表达矩阵exp> exp dim(exp) [1] 54675 22 > exp[1:4,1:4] #检查矩阵是否正常，如果是空的就会报错，空的和有负值的、有异常值的矩阵需要处理原始数据。是否取过loglog 是否有负值 GSM1366348 GSM1366349 GSM1366350 GSM1366351 1007_s_at 279.156 202.866 406.818 232.668 1053_at 487.700 409.018 393.810 397.127 117_at 666.869 388.589 704.633 953.481 121_at 240.646 361.198 305.229 374.757 > #如果表达矩阵为空，大多数是转录组数据，不能用这个流程(后面另讲)。 > #自行判断是否需要log 可以画boxplot() > exp = log2(exp+1) 正常log在0-20之间 > boxplot(exp)正常表达矩阵未取log的表达矩阵有异常样本的表达矩阵处理异常样本1.删掉异常样本2.把异常的拉回正常的exp=limma::normalizeBetweenArrays(exp)关于表达矩阵中有负值1.去过log，有负值-正常2.没取过loglog，有负值-错误数据3.有一半负值-做了标准化（2，3一般弃用除非处理原始数据）### (2)提取临床信息 >pd p = identical(rownames(pd),colnames(exp));p#判断信息是否以一对应 >if(!p) exp = exp[,match(rownames(pd),colnames(exp))] #分组信息来自临床信息，分组信息需要与表达矩阵列名一一对应 #临床信息需要与表达矩阵一一对应(4)提取芯片平台编号>gpl_number save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata") #保存 gse_number(原本的编号),pd(临床信息),exp(表达矩阵),gpl_number(芯片编号)图1Group(实验分组)和ids(探针注释)# 从临床样本中获得实验分组（在表格中慢慢找，代码如何实现看下） rm(list = ls()) load(file = "step1output.Rdata") library(stringr) # 标准流程代码是二分组，多分组数据的分析后面另讲 # 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if # 1.Group---- # 第一种方法，有现成的可以用来分组的列 Group = pd$`disease state:ch1` #pd$`cell type:ch1`注意这个``不能删 }else if(F){ # 第二种方法，自己生成，但需要看清哪是RA哪是control Group = c(rep("RA",times=13), rep("control",times=9)) Group = rep(c("RA","control"),times = c(13,9)) #简写 }else if(T){ # 第三种方法，使用字符串处理的函数获取分组 Group=ifelse(str_detect(pd$source_name_ch1,"control"), "control", "RA") #felse(str_detect()) 两个函数连用是用来分组的神器 #str_detect(pd$source_name_ch1,"control"）-source_name_ch1这一列是否含有control #ifelse(str_detect(pd$source_name_ch1,"control"),"control","RA") 是control就填control不是填RA #三种代码完整方法 if(F){ # 1.Group---- # 第一种方法，有现成的可以用来分组的列 Group = pd$`disease state:ch1` }else if(F){ # 第二种方法，自己生成 Group = c(rep("RA",times=13), rep("control",times=9)) Group = rep(c("RA","control"),times = c(13,9)) }else if(T){ # 第三种方法，使用字符串处理的函数获取分组 Group=ifelse(str_detect(pd$source_name_ch1,"control"), "control", "RA") } # 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后；因子正文与levels不对应时会产生NA Group = factor(Group,levels = c("control","RA")) Group #Group是一个有重复值的向量是分类型数据，适合用因子的形式 #factor直接转换并自动生成levels (control和RA)，顺序以字母排序为准 #levels顺序有意义，在第一个位置的水平是参考水平 #参考水平将在做差异分析时，被设为对照组 #所以需要控制levels的顺序 #levels = c("control","RA") 写了按照写的顺序，control位参考水平三种方法只运行最后一种了，因为有if探针注释

注释来源:不是所有的GPL都可以找到注释

Biocoductor的注释包GPL的表格文件解析 3.官网下载对应产品的注释表格 4.自主注释 https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA Biocoductor的注释包GPL的表格文件自主注释探针注释的获取首先是捷径方法-生信技能树写好的代码>library(tinyarray) #tinyarry这个包是生信技能树写的，find_anno()是其中的一个函数 >find_anno(gpl_number) #打出找注释的代码 [1] "`library(hgu133plus2.db);ids

【本文地址】

公司简介

联系我们