day 7-8 GEO数据挖掘

rm(list = ls())
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) 
options(scipen = 20)#不要以科学计数法表示

#传统下载方式
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)#下载及/或读取文件
#网速太慢，下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下，放在工作目录里
#2.试试geoChina,只能下载2019年前的表达芯片数据
#library(AnnoProbe)
#eSet = geoChina("GSE7305") #选择性代替第8行
#研究一下这个eSet
class(eSet)
length(eSet)

eSet = eSet[[1]] #使得列表的第一个元素重新赋值给eSet
class(eSet)

#(1)提取表达矩阵exp
exp <- exprs(eSet)#获取对象的表达矩阵
#⭐第一个要检查的地方👇，表达矩阵行列数，正常是几万行，列数=样本数，几百行比较难分析出东西
#如果0行说明不是表达芯片或者是遇到特殊情况，不能用此流程分析
dim(exp)#检查行列的函数
#⭐二个要检查的地方👇
range(exp)#看数据范围决定是否需要log，是否有负值，异常值，如有负值，结合箱线图进一步判断
#⭐可能要修改的地方👇
exp = log2(exp+1) #需要log才log，如果不需要log要注释掉这一句，exp+1是为了避免出现负值和0
#⭐第三个要检查的地方👇
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,"Normal");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(hgu133plus2.db))BiocManager::install("hgu133plus2.db",ask = F,update = F)
library(hgu133plus2.db)
ls("package:hgu133plus2.db") #列出R包里都有啥
ids <- toTable(hgu133plus2SYMBOL) #把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的转换）》