适用背景
虽然单细胞绝大多数时候都关注差异表达基因（DEG），但在分析中途会发现有些奇奇怪怪的基因会在所有细胞中都高表达，所以这篇文章才写了下面这个函数ave_top，用于选取每个亚群的高表达基因。
这个函数主要使用了Seurat自带的AverageExpression函数对每个亚群计算其平均表达量，然后排序后选取高表达的top基因，用起来也比较简单。参数介绍：
- obj Seurat对象
- assays 默认为‘RNA’，一般也推荐使用原始的表达矩阵assay
- group.by 分组信息，默认为'seurat_clusters'
library(Seurat)
library(dplyr)
ave_top <- function(obj,assays='RNA',group.by ='seurat_clusters',ntop=20) {
ave <- AverageExpression(obj, assays = assays, group.by = group.by)
df1 <- data.frame(ave[[assays]])
col_sort <- function(col,df,ntop=20) {
df <- df[col]
df$gene <- rownames(df)
df$col <- col
colnames(df)<- c("score","gene","cluster")
rownames(df) <- NULL
df <- arrange(df,desc(score))
df <- df[1:ntop,]
}
lt <- lapply(colnames(df1),col_sort,df=df1,ntop=ntop)
df2 <- Reduce(rbind,lt)
return(df2)
}
函数返回的的是一个数据框，里面记录了每个cluster的前n个top基因平均表达量。
score gene cluster
1 584.010949081858 transcript/58277.mrna1 X0
2 374.43862301027 transcript/58511.mrna1 X0
3 118.511179761744 transcript/42113.mrna1 X0
4 104.396406838368 transcript/52106.mrna1 X0
5 88.9929905206066 transcript/55348.mrna2 X0
6 46.2031101527192 RBFOX1 X0
7 42.0222256206988 IL1RAPL1 X0
8 39.0578761072863 transcript/48150.mrna1 X0
9 38.8097305563072 unknown-transcript-1 X0
10 32.4727604332302 NRXN1 X0
小结与补充
这个函数大多数时候应该用不上，因为只看亚群的高表达基因的意义不大，毕竟高表达基因往往并不是亚群特异性基因，研究单细胞往往也只关注这些特异性基因，因此大家一般都只关注差异表达基因。
