基因注释函数

2022-05-23  本文已影响0人  郭师傅

基本思想:
1、有些ensembl id里包含".",按目前说法需要去掉
2、用clusterprofiler包进行ID转换,然后向表达矩阵进行注释,定义成函数,以后更省事。

# 去掉ID中的"."
exprSet <- exprSet %>% 
  tidyr::separate(id,into = c("id"),sep="\\.") %>% 
  mutate(rowMean = rowMeans(.[,-1])) %>% 
  arrange(desc(rowMean)) %>% 
  distinct(id,.keep_all = T) %>% 
  select(-rowMean) %>% 
  tibble::column_to_rownames(var = "id")
  
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)

geneid1 <- bitr(rownames(exprSet),
               fromType = "ENSEMBL",
               toType = "SYMBOL",
               OrgDb = org.Hs.eg.db)

                          
# 注释函数
# exprSet:表达矩阵
# df:注释文件,两列,1,ID,2,symbol

annot <- function(exprSet,df){
  probe2symbol_df <- df 
  #以下为重复代码,无需修改
  colnames(probe2symbol_df) <- c("probe_id","symbol")
  exprSet <- as.data.frame(exprSet)                                                 #change express matrix to dataframe
  exprSet$probe_id <- rownames(exprSet)                                       #make a new column of probe_id by rowname
  exprSet$probe_id <- as.character(exprSet$probe_id)
  
  #match probe_id and gene symbol
  exprSet <- exprSet %>%
    inner_join(probe2symbol_df,by="probe_id") %>%                                 #合并探针的信息
    dplyr::select(-probe_id) %>%                                                                #去掉多余信息
    dplyr::select(symbol, everything()) %>%                                                #重新排列,
    mutate(rowMean = rowMeans(.[grep("GSM", names(.))])) %>%            #求出平均数(这边的.真的是画龙点睛),改成rowMax,最大值
    arrange(desc(rowMean)) %>%                                                             #把表达量的平均值按从大到小排序
    distinct(symbol,.keep_all = T) %>%                                                      # symbol留下第一个
    dplyr::select(-rowMean) %>%                                                              #反向选择去除rowMean这一列
    tibble::column_to_rownames(colnames(.)[1])                                        # 把第一列变成行名并删除
  return(exprSet)
}

exprSet <- annot(exprSet,geneid1)

save(exprSet,pdata,file = ".\\data\\exprSet_annotated.Rdata")
上一篇下一篇

猜你喜欢

热点阅读