小丫画图出品,作者:haitao and jimmy


点击📎FigureYa59Plus_GEO2DEG.zip,下载FigureYaPlus压缩包,包含本文的代码、输入、输出文件。跟FigureYa59vocano放在一起,更方便查看和使用哦~


需求描述

从GEO数据库下载芯片数据,筛选差异表达基因,无缝对接FigureYa59volcano画火山图。


应用场景

适用于GEO数据库里的芯片数据


环境设置

加载R包

library(GEOquery)
library(limma)
library(ggplot2)
library(ggrepel)
library(ggthemes)

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor


输入数据的获得

如果你的数据已经整理成“easy_input_expr.csv”📎easy_input_expr.csv和“easy_input_pheno.csv”📎easy_input_pheno.txt的格式,就可以跳过这步,直接进入“输入数据预处理”;


如果已经整理成“easy_input_limma.csv”📎easy_input_limma.csv的格式,包含所有基因的差异表达分析结果,每行一个基因,第一列是基因名,后面包含log2foldchange、P.value或adj.Pvalue,就直接使用FigureYa59volcano画最美火山图。


从GEO下载表达矩阵和芯片注释文件

此处以GSE17215为例,用芯片平台自带的探针组和gene symbol对应关系:Affymetrix HT Human Genome U133A Array: GPL3921。如果芯片平台自带的注释信息不符合你的要求,需要重注释,请使用FigureYa58lncRNAreannotation,获得探针ID跟gene symbol的对应关系。


First, let’s start by opening an R session and creating a function to return the eset(ExpressionSet) object or the original list object downloaded by the getGEO() function in R


把下載的步驟寫成一個小的function getGEOdata,每次load function之後就可以簡單通過修改GSE ID拿到自己想要的data進行分析了。


getGEOdata <- function(x, GSE.id=FALSE){
  # Use the getGEO() function to download the GEO data for the id stored in x
  GSEdata <- getGEO(x, GSEMatrix=T, AnnotGPL=FALSE)
  # check the data: printing a summary of the expression values for the first 2 columns
  print(summary(exprs(GSEdata[[1]])[,1:2]))
  # Get the eset object
  eset <- GSEdata[[1]]
  # Save the objects generated for future use in the current working directory
  save(GSEdata, eset, file=paste(x, ".RData", sep=""))
  # check whether we want to return the list object we downloaded on GEO or
  # just the eset object with the getGSEobject argument
  if(GSE.id) return(GSEdata) else return(eset)
}

# Store the dataset ids in a vector GEO_dataset just in case you want to loop through several GEO ids
GEO_dataset <- c("GSE17215") #芯片数据在GEO数据库里的ID

# 获得表达矩阵
eset <- getGEOdata(GEO_dataset[1])
exprSet <- as.data.frame(exprs(eset))
# 提取GSM ID作为样品名
GSMID <- pData(eset)[,c(1,2,8,10,12)]
colnames(exprSet) <- GSMID$geo_accession

# 获得探針组跟gene symbol的对应关系
gpl <- getGEO(eset$platform_id[1])
dim(Table(gpl)) #查看一共有多少行探针组

# 到此為止,所需要的表达矩阵和芯片平台注释已經下載成功

整理芯片平台注释信息

此处使用芯片平台自带的探针组和gene symbol对应关系。


如果芯片平台自带的注释信息不符合你的要求,需要重注释,使用FigureYa58lncRNAreannotation,获得探针ID跟gene symbol的对应关系“gplTOgene.csv”文件,就跳过这步,直接进入“Transfer ID to Gene Symbol”。


colnames(Table(gpl)) #查看列名
#head(Table(gpl))

# 提取自己想要的信息,此处第1、11、12列,分别是行名ID和其对应的Gene Symbol和ENTREZ_GENE_ID
annotation <- Table(gpl)[,c(1,11,12)]
i <- sapply(annotation, is.factor) # Change factor to character
annotation[i] <- lapply(annotation[i], as.character)

# 有些探针组对应多個基因,使用“///”分割,當然你可以選擇保留,我的習慣是保留第一個,所以這裡對於其他的冗余信息進行clean
annotation[,2] <- data.frame(sapply(annotation[,2], function(x) unlist(strsplit(x,'///'))[1]),stringsAsFactors = F)[,1]
annotation[,3] <- data.frame(sapply(annotation[,3], function(x) unlist(strsplit(x,'///'))[1]),stringsAsFactors = F)[,1]
# 或者宁缺毋滥,删除对应多个基因的探针组,用下面两行代替上面两行
#annotation <- annotation[!grepl("///", annotation$`Gene Symbol`),]
#dim(annotation)

# 把整理好的探针组ID、gene symbol、ENTREZ ID对应关系保存到文件里
colnames(annotation) <- c("probe_id", "gsym", "entrez_id")
write.csv(annotation, "gplTOgene.csv", row.names = F)

把表达矩阵的探针组ID换成Gene Symbol

需要用到前面整理好的芯片自带的注释,或重注释获得的“gplTOgene.csv”文件。


annotation <- read.csv("gplTOgene.csv")
head(annotation)

exprSet$probe_id <- rownames(exprSet)
express <- merge(x = exprSet, y = annotation[,c(1:2)], by = "probe_id", all.x = T) 
#如果想轉換成ENTREZ_GENE_ID,就用下面这行
#express <- merge(x = exprSet, y = annotation[,c(1:3)], by = "probe_id", all.x = T) 
express$probe_id <- NULL

# 有些gene对应多个探针组,此处保留其中表达量最大值
dim(express)
rowMeans <- apply(express, 1, function(x) mean(as.numeric(x), na.rm = T))
express <- express[order(rowMeans, decreasing =T),]
express <- express[!duplicated(express[, dim(express)[2]]),]  #express的最后一列为gene name
express <- na.omit(express)
dim(express)

rownames(express) <- express[,dim(express)[2]]
express <- express[,-(dim(express)[2])]
head(express)
# 至此,表达矩阵(express)已构建好,保存到文件
write.csv(express, "easy_input_expr.csv", quote = F)

按照实验设计整理分组信息


# 获得样品信息
pdata <- pData(eset)
i <- sapply(pdata, is.factor)
pdata[i] <- lapply(pdata[i], as.character)

# 按照实验设计,构建分组矩阵
group_list <- subset(pdata, select = title) # Sample的分组信息
group_list$condition = rep(c("a","b"), each = 3)
group_list

# 至此,分组矩阵(design)已构建好
write.table(group_list, "easy_input_pheno.txt", quote = F, sep = "\t", row.names = F)

输入数据预处理

输入数据

- easy_input_expr.csv,基因表达矩阵,每行一个基因,以gene symbol命名,每列一个sample

- easy_input_pheno.txt,样品分组,第一列是sample ID,第二列是分组

均一化处理

如果你的数据已经进行过均一化处理,就直接运行下面这段的最后一行


express <- read.csv("easy_input_expr.csv", row.names = 1)

# Have to check whether this matrix is log2 based
# usually, log2 value should be in between 0 - 18
range(express)
# if yes no need to log transfer, if above this range, have to do log transfer.
express <- log2(express)

par(mfrow=c(2,1))
boxplot(express,outline = F) #均一化前

# Normalise data quantiles
library(preprocessCore)
express.norm <- normalize.quantiles(as.matrix(express))

colnames(express.norm) <- colnames(express)
rownames(express.norm) <- rownames(express)
boxplot(express.norm, outline=F) #均一化后

#不需要再做均一化处理,就只运行这一行:
#express.norm <- express

截屏2019-12-09下午9.50.52.png

离群sample的处理

如果不需要处理离群sample,就只运行下面这段的第一行。


group_list <- read.table("easy_input_pheno.txt", sep = "\t", header = T)
group_list

# Correlation between the samples
heatmap(cor(express.norm), scale = "none")
round(cor(express.norm), 2)

截屏2019-12-09下午9.50.59.png


# 通過correlation,可能会發現某个sample与同一組其他sample之间距離很遠,需要删除。
# 此处删除GSM431123
group_list <- group_list[c(-3),] #删除第3个sample
express.norm <- express.norm[,c(-3)] #删除第3列

# Correlation between the samples again
heatmap(cor(express.norm),scale="none")
round(cor(express.norm),2)

Differential expression analysis

Limma for differential expression gene analysis (DEG) (limma powers differential expression analyses)


design <- model.matrix(~ 0 + factor(group_list$condition))
colnames(design) <- levels(factor(group_list$condition))
rownames(design) <- colnames(express.norm)
design

# 构建差异比较矩阵
contrast.matrix <- makeContrasts(b-a, levels = design)
# 至此,差异表达矩阵已构建好

fit <- lmFit(express.norm,design)
fit2 <- contrasts.fit(fit,contrast.matrix)
fit2 <- eBayes(fit2)
# 得到两两差异表达的结果
# b vs. a
x <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P")

#把全部基因的limma分析结果保存到文件
write.csv(x, "easy_input_limma.csv", quote = F)

开始画图

用上面获得的📎easy_input_limma.csv作为输入文件,无缝对接FigureYa59volcano里的代码,画最美火山图。

classic模式

截屏2019-12-09下午8.08.07.png

advanced模式

截屏2019-12-09下午8.08.15.png

富集分析

用上面获得的📎easy_input_limma.csv作为输入文件,无缝对接FigureYa209batchEnrich里的代码,发现最适合你讲故事的注释数据