适用背景
在单细胞分析中，必要的一步是对每个细胞进行注释，这相当于给每个细胞打个标签，或者说起个别名，这些信息都存在@meta.data信息里，所以不用特意去更改细胞名字，Seurat官网也有相关教程进行这些步骤。但是，如果想改基因名字该怎么改呢？
为什么要改基因名字？常用的基因名一般是字母+数字的基因Symbol形式，例如GAD1和FLT1等等，但有时候拿到的矩阵基因名是类似ENS*这种以基因ID的形式命名基因的，这对于后续分析来说不太方便，毕竟大多数研究是直接使用基因Symbol，因此转换后会更好分析和理解。另外，在做跨物种分析时，我们需要进行跨物种的数据整合分析，这个时候就需要进行同源基因转换，然后使用同源基因集去做整合聚类，因此也需要把基因名字进行重命名，例如小鼠抑制性神经元的经典marker基因是Gad1，而人与其同源的基因是GAD1，如果以人作为参考，则需要把小鼠的基因也改成GAD1，类似的，小鼠里的所有能与人同源的基因都需要做这种转换，之后才能同时把人和小鼠的表达矩阵进行整合分析。
怎么改Seurat的基因名？当然，最简单的方法是导出矩阵，对矩阵的基因名进行替换即可，然后用替换后的矩阵再创建Seurat对象，这是比较简单粗暴也是比较容易理解的办法。此外，另一个解决方法是把Seurat对象里所有涉及基因名的地方都进行基因名转换，本文就是基于此想法写了个函数。
注意：使用本篇的内容，建议使用DietSeurat函数尽量精简Seurat对象之后再改基因名，不然可能会有很多的bugs。
函数灵感来源
然而，官网目前并没有现成的函数对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)
更新20230519
RenameGenesSeurat <- function(obj,newnames,gene.use=NULL,de.assay="Spatial") {
# 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$*@counts, @data and @scale.data, @var.features,@reductions$pca@feature.loadings")
lassays <- Assays(obj)
#names(obj@assays)
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)
df1 <- df1[rownames(obj@reductions$pca@feature.loadings),]
all_genenames1 <- df1$v1
newnames1 <- df1$v2
rownames(obj@reductions$pca@feature.loadings) <- newnames1
}
try(obj[[de.assay]]@meta.features <- data.frame(row.names = rownames(obj[[de.assay]])))
return(obj)
}
小结与讨论
以上函数适用于已经完成聚类分析的Seurat对象，也就是含有pca，umap等降维信息的对象，如果是什么都没做，完全就是一个新的只载入矩阵的Seurat对象，建议还是直接对矩阵操作更方便快捷。
目前用这个RenameGenesSeurat_v2函数还没发现什么问题，可能还要遗漏需要改基因名的地方，欢迎评论区补充！
