Seurat查看不同亚群的高表达基因排行（AverageExpr

适用背景

虽然单细胞绝大多数时候都关注差异表达基因（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

小结与补充

这个函数大多数时候应该用不上，因为只看亚群的高表达基因的意义不大，毕竟高表达基因往往并不是亚群特异性基因，研究单细胞往往也只关注这些特异性基因，因此大家一般都只关注差异表达基因。