今天做了这样一个分析还是蛮有意思的 记录一下吧
下载的表达矩阵是这个样子,转换以后多个转录本对应一个symbol
采用aggregate加和,参考基因的FPKM等于各转录本FPKM之和 - 简书以及ensembl_id转换与gene symbol基因名去重复的两种方法 - 知乎
# GSE177029
library(BiocManager)
if (!require("dplyr"))BiocManager::install("dplyr")
if (!require("tidyverse"))BiocManager::install("tidyverse")
if (!require("org.Hs.eg.db"))BiocManager::install("org.Hs.eg.db") #加载人类转录本
if (!require("clusterProfiler"))BiocManager::install("clusterProfiler")#加载转换工具
rm(list=ls())
exp <- fread("GEO177029xianyu/GSE177029_Raw_gene_counts_matrix.txt",data.table = F)
head(exp)[1:6,1:6]
###把转录本转为symbol
ID<- bitr(exp[,1],
fromType = "ENSEMBLTRANS",
toType=c("SYMBOL"),
OrgDb = org.Hs.eg.db)
###合并转录后的结果
exp_ID <- merge(exp,
ID,
by.x="transcript_id",
by.y="ENSEMBLTRANS",
all=F)
exp_ID <- exp_ID[,-1]
colnames(exp_ID)
# [1] "FC-1_FPKM" "FC-2_FPKM" "FC-3_FPKM" "FSLE-1_FPKM"
# [5] "FSLE-2_FPKM" "FSLE-3_FPKM" "MC-1_FPKM" "MC-2_FPKM"
# [9] "MC-3_FPKM" "MC-4_FPKM" "MC-5_FPKM" "MSLE-1_FPKM"
# [13] "MSLE-2_FPKM" "MSLE-3_FPKM" "MSLE-4_FPKM" "MSLE-5_FPKM"
# [17] "FC-1_count" "FC-2_count" "FC-3_count" "FSLE-1_count"
# [21] "FSLE-2_count" "FSLE-3_count" "MC-1_count" "MC-2_count"
# [25] "MC-3_count" "MC-4_count" "MC-5_count" "MSLE-1_count"
# [29] "MSLE-2_count" "MSLE-3_count" "MSLE-4_count" "MSLE-5_count"
# [33] "SYMBOL"
###error因为有重复的名字。这是由于多个转录本可能对应一个基因造成的
rownames(exp_ID) <- exp_ID$SYMBOL
###继续运行
###可以看到有如下基因都是测到多个转录本
head(sort(table(exp_ID$SYMBOL),decreasing = T),100)
# MYO15B STRADA CORO7 MYO19 RPL17 ACY1
# 35 35 27 26 23 22
###数据分为count和FPKM两个 分开
###count
exp_count <- exp_ID[,17:33]
exp_count_gene <- aggregate(exp_count[,1:16], by=list(exp_count$SYMBOL), FUN=sum)
exp_count_gene <- column_to_rownames(exp_count_gene ,'Group.1')
###FPKM
exp_fpkm <- exp_ID[,c(1:16,33)]
###合并所有重复的gene symbol
###基因的FPKM等于各转录本FPKM之和
###主要思路为利用aggregate函数根据symbol列中的相同基因合并基因表达矩阵
###使用aggregate根据symbol列中的相同基因进行合并
exp_fpkm_gene <- aggregate(exp_fpkm[,1:16], by=list(exp_fpkm$SYMBOL), FUN=sum)
exp_fpkm_gene <- column_to_rownames(exp_fpkm_gene ,'Group.1')
save(exp_fpkm_gene,exp_count_gene,file = "GEO177029xianyu/exp.rda")