12.GEO数据集的R语言差异分析和代码解析—1.数据下载和整理
一、举例介绍
本节下载GSE1009数据集，使用limma包进行差异分析举例。
GSE1009
样本量：共6个样本，其中3个为糖尿病肾病(DN)肾小球样本，另3个为正常肾小球样本。
使用芯片：Affymetrix Human Genome U95 Version 2 Array。
平台：GPL8300。
二、差异分析举例：
###第一步：GEO数据下载
>setwd("D:\\Rfile")
#设置工作目录为D 盘的Rfile文件夹，大家可根据自己需要设置工作目录。
>rm(list = ls())
#清除所有变量
>options(stringsAsFactors=F)
#R导入数据时不自动转换为因子变量
##下面需要安装GEOquery包，用于下载GEO数据库的文件，已经安装GEOquery包不需要重复安装，直接调用GEOquery包。
#安装命令如下(本人已经安装了，所以只给出安装命令，去掉前面的注释符号即可安装)：
#if (!requireNamespace("BiocManager", quietly = TRUE))
#install.packages("BiocManager")
#BiocManager::install("GEOquery")
#bioconductor中的包的安装命令经常更新，如果大家怕安装命令已经更新了，现在的命令不能运行，可以百度“GEOquerybioconductor”进入主页，下拉到安装命令处，将网站上的最新安装命令复制即可，这就是最新的安装命令，其他bioconductor中的包的安装是一样的道理。具体如下：
复制上面红色这段即可。
##下面是用GEOquery包下载数据
>library(GEOquery)
#加载GEOquery包
>gset <- getGEO("GSE1009",destdir = ".",AnnotGPL = F,getGPL = F)
#下载GSE1009表达矩阵文件并赋值给gset变量,destdir = "."表示下载后的文件放在什么地方，默认为当前工作目录。AnnotGPL = F表示不下载注释文件，getGPL = F表示不下载平台文件。我们这里为了下载速度，就设置成不下载平台文件和注释文件。
#下载好后打开D盘,可看到Rfile中有一个“GSE1009_series_matrix.txt.gz”文件。
>save(gset,file = 'GSE1009.gset.Rdata')
#将gset(也就是下载后的GSE1009矩阵文件)保存为R文件，文件名为“GSE1009.gset”方便以后调用。
>gset[["GSE1009_series_matrix.txt.gz"]]@annotation
#查看GSE1009_series_matrix矩阵的平台文件
>gpl <- getGEO('GPL8300', destdir=".")
#下载这个平台文件。用于后面基因symbol注释。
总体先看看上面数据下载代码：
代码运行过程及结果：
运行代码后Rfile中的文件：
###第二步：数据整理（将探针ID给为基因symbol）
>a=gset[[1]]
#提取gset中第一个元素（包含基因表达矩阵和注释信息），并赋值给变量a
>dat=exprs(a)
#提取a中的表达矩阵并赋值给dat，这时的表达矩阵的基因名字为探针ID
>head(dat)
#展示表达矩阵的前6行，看下数据是否经过log转换，一般数据在20以内，是经过log转换的，若有成百上千的数据，表示这个数据还需要进一步log处理。
>ex <- dat
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
(qx[6]-qx[1] > 50 && qx[2] > 0) ||
(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
if (LogC) { ex[which(ex <= 0)] <- NA
dat <- log2(ex)
print("log2 transform finished")}else{print("log2 transform not needed")}
#也可以输入上面这段代码，如果数据已经经log后，显示log2 transform not needed；如果数据没有尽行log，需要log程序按照代码尽行log转换，完成后显示log2 transform finished。
再看看转换后的数据：
不是成千上百的啦。>pd=pData(a)
#查看a的临床信息，为后面选择用于分组的变量做准备
>View(pd)
#查看pd
本例中pd共有26个变量(未完全展示)，我们可以看到只有title下的字段可以分为正常、疾病，其余变量下的字段都是一样的，所以我们选择title字段用于后续分组。
Title变量中的记录为control 1a....Diabetes 1a.......，中间用空格分割，我们需要提取出Control和Diabetes作为分组，所以需要用字符分割包来实现。
##安装字符分割包stringr包,命令如下，已经安装过的可以直接调用，不用重复安装。
#install.packages("stringr")
>library(stringr)
#调用stringr包
>group_list=str_split(pd$title,' ',simplify = T)[,1]
#把pd中title变量的字段，按照空格分割，为了给大家展示，我们先运行分割代码，分割后如下图所示，所以我们选择分割后的第一列即为分组状态
>table(group_list)
#计数每一组个数
>gpl1<-Table(gpl)
>save(gpl1,file = 'gpl1.Rdata')
#我们将平台文件转为list格式，并赋值给gpl1，将gpl1保存为R文件，方便以后调用。
>colnames(Table(gpl))
#查看平台文件的列名，我们看到有ID和gene symbol,记住gene symbol在第几列，我们这里在第11列。
>View(gpl1)
#再了解一下平台文件的数据，当然这里大家可以直接选择gene symbol字段也行，主要了解一下symbol中基因symbol值是否只有一个。
我们这里的gene symbol字段中的symbol，有的基因就不止一个名称，///后面有重名，我们需要第一个名字，所以需要用字符分割包再处理一下,原理同上面处理title一样。
>gpl1$`Gene Symbol`=str_split(gpl1[,11],'///',simplify = T)[,1]
现在是不是变成唯一啦。
>probe2symbol_df<-gpl1[,c(1,11)]
#提取平台文件gpl1中的ID和gene symbol，并赋值给probe2symbol_df
>colnames(probe2symbol_df)=c('probe_id','symbol')
#将列名改为probe_id和symbol
#这两步是因为我懒，不想再调整代码，为了方便后面代码运行，我改的，大家也可以不改。
>length(unique(probe2symbol_df$symbol))
#查看symbol为唯一值的个数
>table(sort(table(probe2symbol_df$symbol)))
#查看每个基因出现n次的个数,我们可以看到，symbol出现一次的基因有7050个，出现2次有1493个。。。
>ids=probe2symbol_df[probe2symbol_df$symbol != '',]
#去掉没有注释symbol的探针（其实这里没有注释的探针数量即为上面出现次数最多的基因440，也就是说有440个探针没有symbol）
>ids=probe2symbol_df[probe2symbol_df$probe_id %in% rownames(dat),]
##%in%用于判断是否匹配，
#注意这里dat是gset表达矩阵的数据，这一步就是把平台文件的ID和矩阵中的ID匹配。
>dat=dat[ids$probe_id,]
#取表达矩阵中可以与探针名匹配的那些，去掉无法匹配的表达数据
>ids$mean=apply(dat,1,mean)
#ids新建mean这一列，列名为mean，同时对dat这个矩阵按行操作，取每一行（即每个样本在这个探针上的）的均值，将结果给到mean这一列的每一行，这里也可以改为中位值，median.
>ids=ids[order(ids$symbol,ids$mean,decreasing = T),]
#即先按symbol排序，相同的symbol再按照均值从大到小排列，方便后续保留第一个值。
>ids=ids[!duplicated(ids$symbol),]
#将symbol这一列取出重复项，'!'为否，即取出不重复的项，去除重复的gene
#取median的话，这样就保留每个基因最大表达量结果.最后得到n个基因。
>dat=dat[ids$probe_id,]
#新的ids取出probe_id这一列，将dat按照取出的这一列，用他们的每一行组成一个新的dat
>rownames(dat)=ids$symbol
#把ids的symbol这一列中的每一行给dat作为dat的行名
>View(dat)
现在就得到行名为基因名，列名为样本名字的矩阵啦。
>dat<-dat[-9173,]
但是我们看到矩阵最后一行，是没有symbol名字的，我们把他去掉，数字自己更改。
>save(dat,group_list,file = 'GSE1009.Rdata')
>write.csv(dat,file="GSE1009expressionmetrix_GSE.csv")
#最后我们把结果保存。下一节会讲解差异分析。
现在看看，整体代码如下：
运行后的Rfile中的文件如下：
整理好的数据如下：
