0. 准备
setwd("")
rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
library(limma)
library(tidyverse)
1. 数据下载
# 下载Series Matrix File(s)
gse <- "GSE75214"
if(!file.exists(paste0(gse,"_eSet.Rdata"))) {
GEO_file <- getGEO(gse, # 需要下载的series
destdir = ".", # 文件下载位置，"."表示当前目录
getGPL = T) # 是否下载GPL注释文件
save(GEO_file, file = paste0(gse,"_eSet.Rdata")) # 将下载下来的文件保存为R语言可以处理的格式
}
load(paste0(gse,"_eSet.Rdata"))
# 如果用getGEO函数下载不下来，就用下面代码中的函数。
# library(AnnoProbe)
# if(!file.exists(paste0(gse,"_eSet.Rdata"))) {
# gset <- geoChina(gse) # 下载表达矩阵和临床数据
# ids <- idmap(gpl = "GPL6244", type = "soft", mirror = "tercent")
# save(gset, ids, file = paste0(gse,"_eSet.Rdata")) # 将下载下来的文件保存为R语言可以处理的格式
# }
GEO_file：
下载完成之后记得比较下载的文件大小和网页上的文件大小是否一致或相差不大。
若相差太大，则“删除本地文件”之后，再用R下载，或将网址复制到浏览器下载。
2. 提取数据集中需要的部分
# 提取数据集中需要的部分 GEO_file[[1]] # 提取GEO_file中第一个数据，有的数据有两个平台测量的数据，会有[[]] exp <- exprs(GEO_file[[1]]) # 提取数据中的样本基因表达矩阵 plate <- fData(GEO_file[[1]]) #提取数据中的平台信息 clinical <- pData(GEO_file[[1]]) # 提取数据中的样本临床信息（比如：年龄、性别、是否存活）
AssayData：数据集包括33252个特征（探针），194个样本
Annotation：测序平台为GPL6244
exp：基因表达矩阵
plate：平台信息
有的探针会对应多个基因，有的基因会对应多个探针，有的探针不会对应基因。需要后续处理。
clinical：临床信息
包括分组信息，等等。
检查exp表达矩阵
# 看一下数据分布 exp <- as.matrix(exp) range(exp) #用range函数查看最大值和最小值。若有负值，则数据无法做差异分析。 boxplot(exp) # 作箱线图，查看数据分布情况，中位数基本在一条直线上为好。 #exp = log2(exp+1) # 若数据没有归一化（即log），则使用此行代码（“+1”是为了防止exp = 0）。 #normalizeBetweenArrays(exp) # 若中位数不在一条直线上，使用此行代码转换数据，使中位数在一条直线上。该函数在limma包中。或直接去掉中位数异常的数据。
boxplot：
3.基因注释
用探针去注释基因。但是会出现“一个探针对应多个基因”、“探针不对应基因”、“多个探针对应一个基因”的情况。
有两种方法去注释，一种是用已有的R包，另一种是用原始的平台数据。
方法1. 使用BioconductorR包直接注释
R包的作者已经删除了一个探针对应多个基因和不对应基因的情况。所以只需要考虑多个探针对应一个基因的情况。
得到的ID可能会比exp少，因为有些探针目前无法匹配到基因。
gpl <- GEO_file[[1]]@annotation # “@”符号类似于数据框中的“$”，用于在列表中取对象，输入“@”后使用“Tab”查看对象名称。
得到GPL编号后，在网站 http://www.bio-info-trainee.com/1399.html 中查找GPL对应的BioconductorR包。如下：
gpl organism bioc_package 1 GPL32 Mus musculus mgu74a 2 GPL33 Mus musculus mgu74b 3 GPL34 Mus musculus mgu74c 6 GPL74 Homo sapiens hcg110 7 GPL75 Mus musculus mu11ksuba 8 GPL76 Mus musculus mu11ksubb 9 GPL77 Mus musculus mu19ksuba 10 GPL78 Mus musculus mu19ksubb 11 GPL79 Mus musculus mu19ksubc 12 GPL80 Homo sapiens hu6800 13 GPL81 Mus musculus mgu74av2 14 GPL82 Mus musculus mgu74bv2 15 GPL83 Mus musculus mgu74cv2 16 GPL85 Rattus norvegicus rgu34a 17 GPL86 Rattus norvegicus rgu34b 18 GPL87 Rattus norvegicus rgu34c 19 GPL88 Rattus norvegicus rnu34 20 GPL89 Rattus norvegicus rtu34 22 GPL91 Homo sapiens hgu95av2 23 GPL92 Homo sapiens hgu95b 24 GPL93 Homo sapiens hgu95c 25 GPL94 Homo sapiens hgu95d 26 GPL95 Homo sapiens hgu95e 27 GPL96 Homo sapiens hgu133a 28 GPL97 Homo sapiens hgu133b 29 GPL98 Homo sapiens hu35ksuba 30 GPL99 Homo sapiens hu35ksubb 31 GPL100 Homo sapiens hu35ksubc 32 GPL101 Homo sapiens hu35ksubd 36 GPL201 Homo sapiens hgfocus 37 GPL339 Mus musculus moe430a 38 GPL340 Mus musculus mouse4302 39 GPL341 Rattus norvegicus rae230a 40 GPL342 Rattus norvegicus rae230b 41 GPL570 Homo sapiens hgu133plus2 42 GPL571 Homo sapiens hgu133a2 43 GPL886 Homo sapiens hgug4111a 44 GPL887 Homo sapiens hgug4110b 45 GPL1261 Mus musculus mouse430a2 49 GPL1352 Homo sapiens u133x3p 50 GPL1355 Rattus norvegicus rat2302 51 GPL1708 Homo sapiens hgug4112a 54 GPL2891 Homo sapiens h20kcod 55 GPL2898 Rattus norvegicus adme16cod 60 GPL3921 Homo sapiens hthgu133a 63 GPL4191 Homo sapiens h10kcod 64 GPL5689 Homo sapiens hgug4100a 65 GPL6097 Homo sapiens illuminaHumanv1 66 GPL6102 Homo sapiens illuminaHumanv2 67 GPL6244 Homo sapiens hugene10sttranscriptcluster 68 GPL6947 Homo sapiens illuminaHumanv3 69 GPL8300 Homo sapiens hgu95av2 70 GPL8490 Homo sapiens IlluminaHumanMethylation27k 71 GPL10558 Homo sapiens illuminaHumanv4 72 GPL11532 Homo sapiens hugene11sttranscriptcluster 73 GPL13497 Homo sapiens HsAgilentDesign026652 74 GPL13534 Homo sapiens IlluminaHumanMethylation450k 75 GPL13667 Homo sapiens hgu219 76 GPL15380 Homo sapiens GGHumanMethCancerPanelv1 77 GPL15396 Homo sapiens hthgu133b 78 GPL17897 Homo sapiens hthgu133a
重试
# 第一种方式：使用BioconductorR包直接注释
if(!require(hugene10sttranscriptcluster.db)) # 判断式安装R包（前缀 + .db）。在RStudio中修改代码时最好使用"Ctrl + F"查找并全部替换。
BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db") # 查看R包里有哪些函数/数据。
ID <- toTable(hugene10sttranscriptclusterSYMBOL) # 提取R包中的SYMBOL数据，即芯片平台的探针ID和对应的基因。R包的作者已经删除了一个探针对应多个基因和不对应基因的情况。
colnames(ID)[1] <- "probe_id"
exp <- as.data.frame(exp)
exp$probe_id <- rownames(exp) # 将探针ID添加到表达矩阵的新的一列
exp <- merge(exp, ID, by = "probe_id") # merge函数将exp的探针id与芯片平台的探针id相匹配。merge()将二者的共同列置于第一列，其余列依次置于之后。
exp[exp == ""] <- NA # 将空白赋值NA
exp <- na.omit(exp) # 删除gene_symbol缺失行的数据
exp <- as.data.frame(exp)
write.csv(exp, "exp.bioc.csv")
exp <- read.csv("exp.bioc.csv", row.names = 1)
exp[exp < 0] <- 0 # 防止后续重复基因取平均值时得到0
ID：芯片平台的探针ID和对应的基因。
exp：增加了基因symbol的表达矩阵
基因注释完成后，还需要处理重复基因，即多个探针对应一个基因的情况，下面的3种方式选择一种就可以。
# 几种处理重复基因的方式，下面的3种方式选择一种就可以
table(duplicated(exp$symbol)) # 看一下有多少重复，TRUE为重复数。duplicated()第二次出现的基因返回的结果为FALSE。
# 第一种，取重复基因的平均值
exp1 <- limma::avereps(exp[,-c(1,ncol(exp))], ID = exp$symbol) # 将重复基因替换为平均值
exp1 <- as.data.frame(exp1)
rownames(exp1) <- exp1$symbol # 加上行名
exp1 <- exp1[, -c(1,ncol(exp1))] # 剔除多余列
write.csv(exp1, "exp.bioc_average.csv") # 保存为csv格式
# 第二种，取重复基因中的最大值
exp2 <- aggregate(exp,
by = list(exp$symbol), # 将数据按symbol拆分为多个子集
FUN = max, # 计算每个子集的最大值，
na.rm = TRUE) # 去除NA
rownames(exp2) <- exp2$Group.1 # 加上行名
exp2 <- exp2[, -c(1,2,ncol(exp2))] # 剔除多余列
write.csv(exp2, "exp.bioc_max.csv") # 保存为csv格式
# 第三种，取重复基因中的最小值
exp3 <- aggregate(exp,
by = list(exp$symbol), # 将数据按symbol拆分为多个子集
FUN = max, # 计算每个子集的最小值，
na.rm = TRUE) # 去除NA
rownames(exp3) <- exp3$Group.1 # 加上行名
exp3 <- exp3[, -c(1,2,ncol(exp3))] # 剔除多余列
write.csv(exp3, "exp.bioc_min.csv") # 保存为csv格式
exp1：
exp2：
exp3：
方法2. 利用下载到的平台数据进行注释
这种方法会出现“一个探针对应多个基因”、“探针不对应基因”、“多个探针对应一个基因”的情况。需要一一处理。
对于“一个探针对应多个基因”的问题，有两种方式处理，可以只提取一个探针对应的多个基因中的第一个基因，也可以将探针对应的所有基因都提取出来。
“探针不对应基因”和“多个探针对应一个基因”的处理方式与“方法1. 使用BioconductorR包直接注释”相同。
方式1. 提取多个基因中的第一个
ID <- data.frame(ID_REF = plate$ID, Gene_Symbol = plate$gene_assignment) # 将平台文件的ID列和基因SYMBOL列取出
ID：出现了“一个探针对应多个基因”、“探针不对应基因”、“多个探针对应一个基因”的情况。
对于上图中一个探针对应多个基因的情况，Gene_Symbol的数据格式是：基因之间以'///'区分，基因说明用"//"区分。
每个基因说明包含5个部分：
①基因ID（包括Ensembl ID，refseq ID，具体区别详见 https://www.jianshu.com/p/a57c78bceb1a ）；
②基因symbol（Gene Official Name，是大家更愿意接受和理解的一种基因名）；
③官方全名（Official Full Name）；④染色体位置（Location）；⑤基因ID（Gene ID）
提取多个基因中的第一个，将其余基因舍去：
# 提取多个基因中的第一个
ID$Gene_Symbol <- sapply(ID$Gene_Symbol,
function(x)
unlist(strsplit(x, "//"))[2]) # [2]表示取第一个基因的Symbol
处理后的ID（和上面的ID对比）：
这样就解决了一个探针对应多个基因的问题。下面解决探针不对应基因的问题：
exp <- as.data.frame(exp)
exp$ID_REF <- rownames(exp) # 将探针ID添加到表达矩阵的新的一列
exp <- merge(exp, ID, by = "ID_REF") # merge函数将exp的探针id与芯片平台的探针id相匹配。merge()将二者的共同列置于第一列，其余列依次置于之后。
exp[, grep("Gene_Symbol", colnames(exp))] <- trimws(exp[, grep("Gene_Symbol", colnames(exp))]) # 去除数据头尾空格。grep("Gene_Symbol", colnames(exp))用于查找Gene_Symbol在exp的第几列。trimws用于从字符串中删除首尾空格。
#exp$Gene_Symbol <- trimws(exp$Gene_Symbol) # 功能同上，去除数据头尾空格
exp[exp == ""] <- NA # 将空白赋值NA
exp <- na.omit(exp) # 删除Gene_Symbol缺失的数据，解决探针不对应基因的问题
exp <- as.data.frame(exp)
write.csv(exp, "exp.first.csv")
exp <- read.csv("exp.first.csv", row.names = 1)
exp[exp < 0] <- 0
exp：
这样就解决了“一个探针对应多个基因”和“探针不对应基因”的问题。
接下来处理重复基因，与“方法1. 使用BioconductorR包直接注释”中的处理方式相同：
# 几种处理重复基因的方式，下面的3种方式选择一种就可以
table(duplicated(exp$Gene_Symbol)) # 看一下有多少重复，TRUE为重复数。duplicated()第二次出现的基因返回的结果为FALSE。
# 第一种，取重复基因的平均值
exp1 <- limma::avereps(exp[,-c(1,ncol(exp))], ID = exp$Gene_Symbol) # 将重复基因替换为平均值
exp1 <- as.data.frame(exp1)
rownames(exp1) <- exp1$Gene_Symbol # 加上行名
exp1 <- exp1[, -c(1,ncol(exp1))] # 剔除多余列
write.csv(exp1, "exp.first_average.csv") # 保存为csv格式
# 第二种，取重复基因中的最大值
exp2 <- aggregate(exp,
by = list(exp$Gene_Symbol), # 将数据按Gene_Symbol拆分为多个子集
FUN = max, # 计算每个子集的最大值，
na.rm = TRUE) # 去除NA
rownames(exp2) <- exp2$Group.1 # 加上行名
exp2 <- exp2[, -c(1,2,ncol(exp2))] # 剔除多余列
write.csv(exp2, "exp.first_max.csv") # 保存为csv格式
# 第三种，取重复基因中的最小值
exp3 <- aggregate(exp,
by = list(exp$Gene_Symbol), # 将数据按Gene_Symbol拆分为多个子集
FUN = max, # 计算每个子集的最小值，
na.rm = TRUE) # 去除NA
rownames(exp3) <- exp3$Group.1 # 加上行名
exp3 <- exp3[, -c(1,2,ncol(exp3))] # 剔除多余列
write.csv(exp3, "exp.first_min.csv") # 保存为csv格式
exp1：
exp2：
exp3：
方式2. 将每一个基因都对应出来
ID <- data.frame(ID_REF = plate$ID, Gene_Symbol = plate$gene_assignment) # 将平台文件的ID列和基因SYMBOL列取出
ID：出现了“一个探针对应多个基因”、“探针不对应基因”、“多个探针对应一个基因”的情况。
对于上图中一个探针对应多个基因的情况，Gene_Symbol的数据格式是：基因之间以'///'区分，基因说明用"//"区分。
每个基因说明包含5个部分：
①基因ID（包括Ensembl ID，refseq ID，具体区别详见 https://www.jianshu.com/p/a57c78bceb1a ）；
②基因symbol（Gene Official Name，是大家更愿意接受和理解的一种基因名）；
③官方全名（Official Full Name）；④染色体位置（Location）；⑤基因ID（Gene ID）
将每一个基因都对应出来：
### 将每一个基因都对应出来
# 将探针ID和基因symbol粘到一起，中间用"..."分隔
x <- tibble(unlist(apply(ID, 1,function(x){ # 对ID的每一行都执行function函数
paste(x[1], # 探针ID
if(str_detect(as.character(x[2]), "/")) # 当ID第二列内容有"/"时，执行下列操作，否则，不执行操作。
str_split(str_split(x[2],"///",simplify = TRUE), # 先用"///"将基因分开
"//", # 再用"//"将基因内部的信息分开
simplify = TRUE)[,2], # 取第二项基因symbol
sep = "...")})))
colnames(x) <- "ABC" # 修改x的列名
file <- separate(x, ABC, c("ID_REF", "Gene_Symbol"), sep = "\\...") # 将探针ID和基因symbol分开
x：
file（和上面的ID对比）：一个探针对应的多个基因都被提取出来了，但其中有些基因是重复的，后续再进行处理。
这样就把一个探针对应的多个基因都提取出来了。下面解决探针不对应基因的问题：
exp <- as.data.frame(exp)
exp$ID_REF <- rownames(exp) # 将探针ID添加到表达矩阵的新的一列
exp <- merge(exp, file, by = "ID_REF") # merge函数将exp的探针id与芯片平台的探针id相匹配。merge()将二者的共同列置于第一列，其余列依次置于之后。
exp[, grep("Gene_Symbol", colnames(exp))] <- trimws(exp[, grep("Gene_Symbol", colnames(exp))]) # 去除数据头尾空格。grep("Gene_Symbol", colnames(exp))用于查找Gene_Symbol在exp的第几列。trimws用于从字符串中删除首尾空格。
#exp$Gene_Symbol <- trimws(exp$Gene_Symbol) # 功能同上，去除数据头尾空格
exp[exp == ""] <- NA # 将空白赋值NA
exp <- na.omit(exp) # 删除Gene_Symbol缺失的数据
exp <- as.data.frame(exp)
write.csv(exp, "exp.all.csv")
exp <- read.csv("exp.all.csv", row.names = 1)
exp[exp < 0] <- 0
exp：
这样就将探针对应的多个基因都提取出来，并解决了“探针不对应基因”的问题。
接下来处理重复基因，与“方法1. 使用BioconductorR包直接注释”中的处理方式相同：
# 几种处理重复基因的方式，下面的3种方式选择一种就可以
table(duplicated(exp$Gene_Symbol)) # 看一下有多少重复，TRUE为重复数。duplicated()第二次出现的基因返回的结果为FALSE。
# 第一种，取重复基因的平均值
exp1 <- limma::avereps(exp[,-c(1,ncol(exp))], ID = exp$Gene_Symbol) # 将重复基因替换为平均值
exp1 <- as.data.frame(exp1)
rownames(exp1) <- exp1$Gene_Symbol # 加上行名
exp1 <- exp1[, -c(1,ncol(exp1))] # 剔除多余列
write.csv(exp1, "exp.all_average.csv") # 保存为csv格式
# 第二种，取重复基因中的最大值
exp2 <- aggregate(exp,
by = list(exp$Gene_Symbol), # 将数据按Gene_Symbol拆分为多个子集
FUN = max, # 计算每个子集的最大值，
na.rm = TRUE) # 去除NA
rownames(exp2) <- exp2$Group.1 # 加上行名
exp2 <- exp2[, -c(1,2,ncol(exp2))] # 剔除多余列
write.csv(exp2, "exp.all_max.csv") # 保存为csv格式
# 第三种，取重复基因中的最小值
exp3 <- aggregate(exp,
by = list(exp$Gene_Symbol), # 将数据按Gene_Symbol拆分为多个子集
FUN = max, # 计算每个子集的最小值，
na.rm = TRUE) # 去除NA
rownames(exp3) <- exp3$Group.1 # 加上行名
exp3 <- exp3[, -c(1,2,ncol(exp3))] # 剔除多余列
write.csv(exp3, "exp.all_min.csv") # 保存为csv格式
exp1：
exp2：
exp3：
两种基因注释方法的对比
方法1. 使用BioconductorR包直接注释
这种方式得到的基因数目最少，可能是由于R包比较旧，有些新的探针没有加上去。如果对数据要求比较高，可以用方法2。
方法2. 利用下载到的平台数据进行注释
方式1. 提取多个基因中的第一个
这种方式得到的基因数目较方式2少，是因为每个探针只取了第一个对应的基因，所以探针对应的第二个基因之后的基因是没有的。但是也不会有太大
方式2. 将每一个基因都对应出来
这种方式得到的基因数目较方式1多，但由于探针对应的基因都被提取出来，所以特异性可能不高。这种情况应该是探针对应的多个基因有重叠序列。
4.其他数据下载方式
前面的数据下载方式使用代码下载了表达矩阵和平台数据，持此之外，还可以直接从官网下载。
1. 官网下载表达矩阵和平台信息
NCBI官网（https://www.ncbi.nlm.nih.gov/）或GEO官网（https://www.ncbi.nlm.nih.gov/geo/）搜索 GSE75214：
下载表达量数据：点击 Series Matrix File(s) 进行下载。下载之后解压到工作目录：GSE75214_series_matrix.txt
下载平台数据：点击 Platforms (1) GPL6244 
进行下载：GPL6244-17930.txt
然后用代码加载数据：
# 将数据加载好
exp <- read.table("GSE75214_series_matrix.txt", header = TRUE, sep = "\t", dec = ".",
comment.char = "!", na.strings = c("NA"), fill = TRUE) # 读取基因列表矩阵
plate <- read.table("GPL6244-17930.txt", header = TRUE, quote = "", sep = "\t", dec = ".",
comment.char = "#", na.strings = c("NA"), fill = TRUE) # 读取平台文件（GPL）
exp：
plate：
后续分析与前面是一样的。只不过，下载的文件是矩阵格式，无法保存为eSet格式的数据，也就无法用函数获得其临床数据。
2. 官网下载原始数据
下载原始数据：点击（http）进行下载：GSE75214_RAW.tar
将下载的tar文件解压到文件夹 GSE75214_RAW 下。
下载平台数据：点击 Platforms (1) GPL6244 进行下载：GPL6244-17930.txt
library(tidyverse)
library(limma)
library(GEOquery)
library(affy)
# 读取表达数据
setwd("D:/BioInformationAnalyze/GEOdata/GEO_UC/my_analysis/1.Data/1.TestSet/GSE75214/GSE75214_RAW")
exp <- data.frame(t(data.frame(rma(ReadAffy())))) # rma：标准化处理。ReadAffy：读取CEL格式的文件。
# 看一下处理后的数据分布
exp <- as.matrix(exp)
boxplot(exp)
exp <- as.data.frame(exp)
# 读取平台数据
setwd("D:/BioInformationAnalyze/GEOdata/GEO_UC/my_analysis/1.Data/1.TestSet/GSE75214")
plate <- read.table("GPL6244-17930.txt", header = TRUE,
quote = "", sep = "\t", dec = ".",
comment.char = "#", na.strings = c("NA"), fill = TRUE) # 读取平台文件（GPL）
ID <- data.frame(ID_REF = plate$ID, Gene_Symbol = plate$gene_assignment) # 将平台文件的ID列和基因SYMBOL列取出
ID$ID_REF <- paste("X", ID$ID_REF,sep = "") # 如果exp中的ID_REF和ID中的ID_REF不一致时，需要设置为一致。视情况而定。
exp$ID_REF <- rownames(exp)
rawdata还没下载完，下载完后再说。
翻译失败: Failed to fetch
