R语言通路基因ID提取表达table

生信喵

背景

用基因ID提取基因表达量，输出csv文件，作为下一步分析的输入文件。

应用场景

场景一：GO富集分析获得了GO term里的多个基因，需要获得这些基因的表达矩阵。

场景二：手里有已分类的基因列表，想提取相应的表达矩阵。

下一步，根据实际需要，用下面这些代码进行可视化：

• 绘制heatmap，为单个GO term对应的基因进行聚类分析。
• 已分类的heatmap绘制多个GO term对应基因的heatmap
• 批量绘制单个基因的box plot
• 批量绘制表达谱曲线

输入文件

需要两种文件：

• 基因列表文件，not_easy_input.txt或 GO*.txt
• 表达矩阵文件，not_easy_input_expr.txt或 easy_input_expr.txt。

基因表达矩阵

每行一个基因，每列一个sample

exprSet<-read.table("not_easy_input_expr.txt",as.is = T)
exprSet[1:3,1:5]

如果你不需要排序，可以跳过下面这段，直接进入“基因列表”。

下面按照表达量由高到低给基因排序：

library(dplyr)
library(tidyr)
library(tibble)
exprSet <- exprSet %>% 
  #把行名变成一列，命名为symbol
  rownames_to_column(var = "symbol")%>% 
  #新建一列rowMeans，每一行求平均值并填充进去
  dplyr::mutate(rowMeans =rowMeans(.[grep("TCGA", names(.))])) %>%
  #按照求得的平均值把其他列从高到低排序
  dplyr::arrange(desc(rowMeans)) %>% 
  #把symbol这一列变成行名
  tibble::column_to_rownames(var = "symbol")%>% 
  #去掉平均值那一列
  dplyr::select(-rowMeans)

write.table(exprSet[,1:10],"easy_input_expr.txt",quote = F)

easy_input_expr.txt是按基因表达量从高到低排好序的表达矩阵。

基因列表

如果你的基因列表已经整理成 GO*.txt那样，就可以跳过这步，直接进入“开始提取”。

此处的输入文件是GO富集分析的输出文件：not_easy_input_GO.txt，每行一个GO term，GO terms对应的基因位于第8列，并且以“/”号分隔。

以下代码同样适用于其他富集分析结果，如果文件里基因名以“, ”分隔，只需把下面代码里的 /改为 , 。

ego_BP_df<- read.table("not_easy_input_GO.txt",sep = "\t")
ego_BP_df[1:3,]

#输出排名靠前的GO term里的基因
#此处输出前两个
for (i in 1:2){
  #此处用"/"分隔基因列表，这取决于你的基因ID之间用的是哪个分隔符
  genelist <- unlist(strsplit(as.character(ego_BP_df[i,]$geneID),"/")) 
  write.table(genelist,paste0(rownames(ego_BP_df)[i],".txt"),row.names = F,col.names = F,quote = F)
}

#或者提取特定的几个GO term
GO_id <- c("GO:0060333","GO:0050900") #把你想要提取的GO term放在这里
for (i in 1:length(GO_id)){
  index <- grep(GO_id[i],ego_BP_df$ID) #在总GO term列表中的位置
  genelist <- unlist(strsplit(as.character(ego_BP_df[index,]$geneID),"/"))
  write.table(genelist,paste0(GO_id[i],".txt"),row.names = F,col.names = F,quote = F)
}

到这里，基因名就被整理成 GO*.txt文件的格式：

基因列表保存在多个文件里，每个文件包含一个GO term的基因ID。

基因名呈一列，每行一个基因名。

开始提取

读入基因表达矩阵

exprSet<-read.table("easy_input_expr.txt",as.is = T)
exprSet[1:3,1:5]

从基因列表文件逐一提取表达矩阵，保存到文件

#按照文件名的规律，读取基因列表文件
#即使你只有一个基因ID列表文件，也可以这样操作
fnames<-Sys.glob("GO*.txt") 

for (i in 1:length(fnames)){
  genelist<-read.table(fnames[i])
  #按表达矩阵中基因的顺序排列
  genelist <- rownames(exprSet)[rownames(exprSet) %in% genelist$V1]
  genelist_expr <- exprSet[genelist,]
  write.csv(genelist_expr,paste0(unlist(strsplit(fnames[i],".txt"))[1],".csv"),quote = F)
}

这里输出的 GO*.csv文件里就是基因表达矩阵。

通常情况下，到这里就结束了。除非：

你要用一条pheatmap命令画出多个GO term里的基因，就要继续运行下面的代码：

重复出现的基因名处理

例如你要画已分类heatmap，一步画出所有GO term里的基因表达谱heatmap。

会遇到报错提示“基因名不唯一”。

那是因为很多基因在多个GO term间重复出现，因此，用下面代码在基因名后面加上数字，以区分来源于不同GO term的同一基因。

fnames<-Sys.glob("GO*.txt")

for (i in 1:length(fnames)){
  genelist<-read.table(fnames[i])
  #按照表达矩阵中的位置排序
  genelist <- rownames(exprSet)[rownames(exprSet) %in% genelist$V1]
  genelist_expr <- exprSet[genelist,]
  #在基因名后面加上“.数字”
  rownames(genelist_expr)<-paste0(rownames(genelist_expr),".",i) 
  write.csv(genelist_expr,paste0(unlist(strsplit(fnames[i],".txt"))[1],".csv"),quote = F)
}

生信喵

文中所用数据可以关注公众号“生信喵实验柴”

发送关键词“20250413”获取