GEO基础

GEO分析也是数据挖掘的一部分，目的在于减少基因数量，锁定核心基因。数据挖掘数据来源有不分方向的GEO和NHANES，及肿瘤专属的TCGA,ICGC,CCLE,SEER。类型可包括基因表达芯片，转录组，单细胞，突变，甲基化及拷贝数变异等等。

作图包括热图，箱线图点图，火山图等，箱线图的含义：

1.logFC

2.主成分分析

3.GEO表达芯片

分析表达芯片的主要目的是获得差异基因，后续再进行功能富集等分析。基因表达芯片的原理在于利用芯片探针与mRNA的核酸杂交荧光机制，探索基因的表达量。

主要有以下标志：

4.代码

表达矩阵的验证

表达矩阵的正常范围为0~20

1.处理异常样本：
第一个办法：删掉异常样本
第二个办法：exp = limma::normalizeBetweenArrays(exp)

2.关于表达矩阵里的负值

取过log，有少量的负值，4<中位数<15 —正常
没取过log，有负值— 错误数据

表达矩阵处理代码

 rm(list = ls()) #打破下载时间的限制,改前60秒，改后10w秒 options(timeout = 100000) options(scipen = 20)#不要以科学计数法表示 #传统下载方式 library(GEOquery) eSet = getGEO("GSE7305", destdir = '.', getGPL = F)#getGEO完成下载和读取 #网速太慢，下不下来怎么办 #1.从网页上下载/发链接让别人帮忙下，放在工作目录里 #2.试试geoChina,只能下载2019年前的表达芯片数据 #library(AnnoProbe) #eSet = geoChina("GSE7305") #选择性代替第8行 #研究一下这个eSet class(eSet) length(eSet) eSet = eSet[[1]] #把唯一的元素取出 class(eSet) #(1)提取表达矩阵exp exp <- exprs(eSet) #⭐第一个要检查的地方👇，表达矩阵行列数，正常是几万行，列数=样本数， #如果0行说明不是表达芯片或者是遇到特殊情况，不能用此流程分析 dim(exp) #⭐二个要检查的地方👇 range(exp)#看数据范围决定是否需要log，是否有负值，异常值，如有负值，结合箱线图进一步判断 #⭐可能要修改的地方👇 exp = log2(exp+1) #需要log才log，不需要log要注释掉这一句 #⭐第三个要检查的地方👇 boxplot(exp,las = 2) #看是否有异常样本 #(2)提取临床信息 pd <- pData(eSet) #⭐多分组中提取两分组的代码示例，二分组不需要 if(F){ #因为现在这个例子不是多分组，所以编造一列做示例。 pd$fake = paste0(rep(c("a","b","c","d"),each = 5),1:5) k1 = str_detect(pd$fake,"b");table(k1) k2 = str_detect(pd$fake,"c");table(k2) pd = pd[k1|k2,] } #(3)让exp列名与pd的行名顺序完全一致 p = identical(rownames(pd),colnames(exp));p if(!p) { s = intersect(rownames(pd),colnames(exp)) exp = exp[,s] pd = pd[s,] } #(4)提取芯片平台编号，后面要根据它来找探针注释 gpl_number <- eSet@annotation;gpl_number save(pd,exp,gpl_number,file = "step1output.Rdata") # 原始数据处理的代码，按需学习 # https://mp.weixin.qq.com/s/0g8XkhXM3PndtPd-BUiVgw 

探针注释：

 # Group(实验分组)和ids(探针注释) rm(list = ls()) load(file = "step1output.Rdata") # 1.Group---- library(stringr) # 标准流程代码是二分组，多分组数据的分析后面另讲 #⭐要修改的地方：分组信息,必须学会ifelse和str_detect k = str_detect(pd$title,"Control");table(k) #不在title就在pd的其他列 Group = ifelse(k,"Normal","Disease") # 需要把Group转换成因子，并设置参考水平，指定levels #⭐要修改的地方，对照组在前，处理组在后 Group = factor(Group,levels = c("Normal","Disease")) Group #⭐检查自己得到的分组是否正确 data.frame(pd$title,Group) #2.探针注释的获取----------------- #四种方法，方法1里找不到就从方法2找，以此类推。 #方法1 BioconductorR包(最常用) #⭐要操作的地方 library(tinyarray) gpl_number #首先看看编号是多少 #View(pkg_all) #然后在pkg_all里搜索gpl编号，找到对应的R包前缀(第二列)，没搜到就是没有R包,再看方法2。 #也可以用代码直接得到对应的R包前缀： pkg_all[pkg_all$gpl==gpl_number,2] #⭐要操作的地方 #用上面找到的前缀替换下面所有的hgu133plus2，共5处 if(!require(hugene10sttranscriptcluster.db))BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F) library(hugene10sttranscriptcluster.db) ls("package:hugene10sttranscriptcluster.db") #列出R包里都有啥 ids <- toTable(hugene10sttranscriptclusterSYMBOL) #把R包里的注释表格变成数据框 # 方法2 下载并读取GPL网页的表格文件，按列取子集 #⭐要操作的地方 library(tinyarray) get_gpl_txt(gpl_number) #获取表格文件的下载链接 # 接下来是复制网址去浏览器下载、放在工作目录下、读取、提取探针id和基因symbol(没有现成的需要拆分和转换)，不同文件代码不统一，等看同学们的例子。 # 注意:最终的数据ids只能有两列，第一列列名是probe_id,第二列列名是symbol,且都是字符型，否则后面代码要报错咯。 # 方法3 官网下载注释文件并读取 # 方法4 自主注释，了解一下 #https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA save(exp,Group,ids,file = "step2output.Rdata") #比较复杂的探针注释参考资料 #资料1：拆分取列https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/sv262capcgg9o8s5 #资料2：多种id的转换https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/pn0s1mmsaxocfynb?singleDoc# 《又一个有点难的探针注释（多种id的转换）》