Day08 生信马拉松-GEO数据挖掘 （上）

options("repos"="https://mirrors.ustc.edu.cn/CRAN/") #设置CRAN镜像
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") #设置Bioconductor镜像

#CRAN中需要安装的包
cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'cowplot',
                   'patchwork',
                   'basetheme',
                   'paletteer',
                   'AnnoProbe',
                   'ggthemes',
                   'VennDiagram',
                   'tinyarray')  

#Bioconductor中需要安装的包
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot")

for (pkg in cran_packages){
  if (! require(pkg,character.only=T,quietly = T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
} #character.only=T是为了精准的识别包名的字符串


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T,quietly = T) ) {
    BiocManager::install(pkg,ask = F,update = F) #不询问不更新，减少麻烦，确保代码顺利运行
    require(pkg,character.only=T) 
  }
}

#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}
#重复运行，如果没有问题控制台则不再出warning信息

#没有任何提示就是成功了，如果有warning xx包不存在，用library检查一下。
#library报错，就单独安装。

###前期准备工作###
rm(list = ls()) #清空环境变量
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) 
options(scipen = 20) #不要以科学计数法表示 scipen的数值代表阈值，默认1从10万开始，每加1阈值扩大10倍

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

eSet = eSet[[1]] #去除list外壳
class(eSet) #由Biobase包里面的“ExpressionSet”对象 

exp <- exprs(eSet) #Biobase中特定提取子集的函数
dim(exp) #看行、列数量 若出现异常“0”，GEO确认原始信息
range(exp) #看数据范围决定是否需要log（有特别大的，超过50的），是否有负值，异常值
# exp = log2(exp+1) #需要log才log，不需要就在这句代码前加#号 “+1”防止负数和“0”值出现
boxplot(exp,las = 2) #看是否有异常样本

pd <- pData(eSet)

p = identical(rownames(pd),colnames(exp));p #鉴定是否一致
if(!p) {
  s = intersect(rownames(pd),colnames(exp)) #取交集找出相同部分
  exp = exp[,s] #提取exp的列名
  pd = pd[s,] #提取pd的行名
}

###如果只有两个分组不需要此段###
k = pd$source_name_ch1 %in% c("Ctrl in adherent culture",
                          "DPN treatment in adherent culture")
#%in%中“”内的内容必须与pd中对应列完全一致
pd = pd[k,]
exp = exp[,k]

identical(rownames(pd),colnames(exp)) #再次确认临床数据的行基因数据的列名是否一致
dim(exp)
dim(pd) 
table(pd$source_name_ch1) #查看提取后的分组名称                                        

gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")