如何给Seurat对象的基因重命名？【基因名转换】2022-07

适用背景

在单细胞分析中，必要的一步是对每个细胞进行注释，这相当于给每个细胞打个标签，或者说起个别名，这些信息都存在@meta.data信息里，所以不用特意去更改细胞名字，Seurat官网也有相关教程进行这些步骤。但是，如果想改基因名字该怎么改呢？

为什么要改基因名字？常用的基因名一般是字母+数字的基因Symbol形式，例如GAD1和FLT1等等，但有时候拿到的矩阵基因名是类似ENS*这种以基因ID的形式命名基因的，这对于后续分析来说不太方便，毕竟大多数研究是直接使用基因Symbol，因此转换后会更好分析和理解。另外，在做跨物种分析时，我们需要进行跨物种的数据整合分析，这个时候就需要进行同源基因转换，然后使用同源基因集去做整合聚类，因此也需要把基因名字进行重命名，例如小鼠抑制性神经元的经典marker基因是Gad1，而人与其同源的基因是GAD1，如果以人作为参考，则需要把小鼠的基因也改成GAD1，类似的，小鼠里的所有能与人同源的基因都需要做这种转换，之后才能同时把人和小鼠的表达矩阵进行整合分析。

怎么改Seurat的基因名？当然，最简单的方法是导出矩阵，对矩阵的基因名进行替换即可，然后用替换后的矩阵再创建Seurat对象，这是比较简单粗暴也是比较容易理解的办法。此外，另一个解决方法是把Seurat对象里所有涉及基因名的地方都进行基因名转换，本文就是基于此想法写了个函数。

函数灵感来源

然而，官网目前并没有现成的函数对Seurat对象进行基因重命名。GitHub上的一个 issue 倒是有解决方案，然而只对obj@assays$RNA@counts, @data and @scale.data这些改了基因名，实测在运行某些函数，例如FindVariableFeatures时就会报错。

事实上，这个函数主要就是对Seurat对象里的各种需要用到基因名的地方进行了基因转换，以下是原始代码：

RenameGenesSeurat <- function(obj = ls.Seurat[[i]], newnames = HGNC.updated[[i]]$Suggested.Symbol) { # Replace gene names in different slots of a Seurat object. Run this before integration. Run this before integration. It only changes obj@assays$RNA@counts, @data and @scale.data.
  print("Run this before integration. It only changes obj@assays$RNA@counts, @data and @scale.data.")
  RNA <- obj@assays$RNA

  if (nrow(RNA) == length(newnames)) {
    if (length(RNA@counts)) RNA@counts@Dimnames[[1]]            <- newnames
    if (length(RNA@data)) RNA@data@Dimnames[[1]]                <- newnames
    if (length(RNA@scale.data)) RNA@scale.data@Dimnames[[1]]    <- newnames
  } else {"Unequal gene sets: nrow(RNA) != nrow(newnames)"}
  obj@assays$RNA <- RNA
  return(obj)
}

但是这个代码只改变了RNA的assays，而且如果修改的assay不是RNA而是integrated的话，修改scale.data那一步会容易报错，这个函数虽然简单也容易理解，但是就是不好用，改完之后的Seurat对象在运行很多函数时都会报错。

升级版函数 RenameGenesSeurat_v2

由于原来的函数太过简单，局限性太多而且，而且有不少bugs，所以基于它的版本我做了一些改进，然后封装成函数RenameGenesSeurat_v2：

• obj，Seurat对象
• newnames，新的基因名字集向量，与gene.use的基因集一一对应
• gene.use，使用的基因集向量，默认使用所有原始基因，如果给定基因集会对obj取基因子集，长度需与newnames一样，并且与newnames对应
• de.assay，默认的assay，单细胞数据填RNA，空间组数据填Spatial

RenameGenesSeurat_v2 <- function(obj,newnames,gene.use=NULL,de.assay="RNA") {
print("Run this before integration. It only changes obj@assays$*@counts, @data and @scale.data, @var.features,@reductions$pca@feature.loadings")
lassays <- Assays(obj)
assay.use <- obj@reductions$pca@assay.used
DefaultAssay(obj) <- de.assay
if (is.null(gene.use)) {
all_genenames <- rownames(obj)
}else{
all_genenames <- gene.use
obj <- subset(obj,features=gene.use)
}

order_name <- function(v1,v2,ref){
v2 <- make.names(v2,unique=T)
df1 <- data.frame(v1,v2)
rownames(df1) <- df1$v1
df1 <- df1[ref,]
return(df1)
}

df1 <- order_name(v1=all_genenames,v2=newnames,ref=rownames(obj))
all_genenames <- df1$v1
newnames <- df1$v2

if ('SCT' %in% lassays) {
if ('SCTModel.list' %in%  slotNames(obj@assays$SCT)) {
obj@assays$SCT@SCTModel.list$model1@feature.attributes <- obj@assays$SCT@SCTModel.list$model1@feature.attributes[all_genenames,]
rownames(obj@assays$SCT@SCTModel.list$model1@feature.attributes) <- newnames
}
}
change_assay <- function(a1=de.assay,obj,newnames=NULL,all_genenames=NULL){
RNA <- obj@assays[a1][[1]]
  if (nrow(RNA) == length(newnames)) {
    if (length(RNA@counts)) RNA@counts@Dimnames[[1]]            <- newnames
    if (length(RNA@data)) RNA@data@Dimnames[[1]]                <- newnames
    if (length(RNA@var.features)) {
        df1 <- order_name(v1=all_genenames,v2=newnames,ref=RNA@var.features)
        all_genenames1 <- df1$v1
        newnames1 <- df1$v2
        RNA@var.features            <- newnames1
    }
    if (length(RNA@scale.data)){
        df1 <- order_name(v1=all_genenames,v2=newnames,ref=rownames(RNA@scale.data))
        all_genenames1 <- df1$v1
        newnames1 <- df1$v2
        rownames(RNA@scale.data)    <- newnames1
    }

  } else {"Unequal gene sets: nrow(RNA) != nrow(newnames)"}
  obj@assays[a1][[1]] <- RNA
  return(obj)
}

for (a in lassays) {
DefaultAssay(obj) <- a
df1 <- order_name(v1=all_genenames,v2=newnames,ref=rownames(obj))
all_genenames1 <- df1$v1
newnames1 <- df1$v2
obj <- change_assay(obj=obj,a1=a,newnames=newnames1,all_genenames=all_genenames1)
}

hvg <- VariableFeatures(obj,assay=assay.use)
if (length(obj@reductions$pca)){
    df1 <- order_name(v1=all_genenames,v2=newnames,ref=hvg)
    all_genenames1 <- df1$v1
    newnames1 <- df1$v2
    rownames(obj@reductions$pca@feature.loadings) <- newnames1
}

return(obj)
}

这个函数不仅仅对单个assay进行基因名转换，而是会对所有的assay都进行转换，另外还会对obj@reductions$pca@feature.loadings进行转换，这个处理会使得FindVariableFeatures正常运行。另外对于做了SCTransform的含有SCT assay的Seurat对象，还转换了obj@assays$SCT@SCTModel.list$model1@feature.attributes。

• 使用示例

ingene <- read.table('gene_human.xls',sep='\t',header=T,stringsAsFactor=F)
ob3 <- RenameGenesSeurat(obj=ob3,newnames=ingene$hs_gene,gene.use=ingene$Axolotl_ID)

小结与讨论

以上函数适用于已经完成聚类分析的Seurat对象，也就是含有pca，umap等降维信息的对象，如果是什么都没做，完全就是一个新的只载入矩阵的Seurat对象，建议还是直接对矩阵操作更方便快捷。
目前用这个RenameGenesSeurat_v2函数还没发现什么问题，可能还要遗漏需要改基因名的地方，欢迎评论区补充！