適用背景
在單細胞分析中，必要的一步是對每個細胞進行注釋扁誓，這相當于給每個細胞打個標簽防泵，或者說起個別名，這些信息都存在@meta.data信息里蝗敢，所以不用特意去更改細胞名字捷泞，Seurat官網(wǎng)也有相關教程進行這些步驟。但是寿谴，如果想改基因名字該怎么改呢锁右？
為什么要改基因名字？常用的基因名一般是字母+數(shù)字的基因Symbol形式讶泰，例如GAD1和FLT1等等咏瑟，但有時候拿到的矩陣基因名是類似ENS*這種以基因ID的形式命名基因的，這對于后續(xù)分析來說不太方便痪署，畢竟大多數(shù)研究是直接使用基因Symbol码泞，因此轉換后會更好分析和理解。另外狼犯，在做跨物種分析時余寥，我們需要進行跨物種的數(shù)據(jù)整合分析领铐，這個時候就需要進行同源基因轉換，然后使用同源基因集去做整合聚類宋舷，因此也需要把基因名字進行重命名绪撵，例如小鼠抑制性神經(jīng)元的經(jīng)典marker基因是Gad1，而人與其同源的基因是GAD1祝蝠，如果以人作為參考莲兢，則需要把小鼠的基因也改成GAD1，類似的续膳，小鼠里的所有能與人同源的基因都需要做這種轉換改艇，之后才能同時把人和小鼠的表達矩陣進行整合分析。
怎么改Seurat的基因名坟岔？當然谒兄，最簡單的方法是導出矩陣，對矩陣的基因名進行替換即可社付，然后用替換后的矩陣再創(chuàng)建Seurat對象承疲，這是比較簡單粗暴也是比較容易理解的辦法。此外鸥咖，另一個解決方法是把Seurat對象里所有涉及基因名的地方都進行基因名轉換燕鸽，本文就是基于此想法寫了個函數(shù)。
注意：使用本篇的內(nèi)容啼辣，建議使用DietSeurat函數(shù)盡量精簡Seurat對象之后再改基因名啊研，不然可能會有很多的bugs。
函數(shù)靈感來源
然而鸥拧，官網(wǎng)目前并沒有現(xiàn)成的函數(shù)對Seurat對象進行基因重命名党远。GitHub上的一個 issue 倒是有解決方案，然而只對obj@assays$RNA@counts, @data and @scale.data這些改了基因名富弦，實測在運行某些函數(shù)沟娱，例如FindVariableFeatures時就會報錯。
事實上腕柜，這個函數(shù)主要就是對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那一步會容易報錯砰蠢，這個函數(shù)雖然簡單也容易理解，但是就是不好用蛾找，改完之后的Seurat對象在運行很多函數(shù)時都會報錯娩脾。
升級版函數(shù) RenameGenesSeurat_v2
由于原來的函數(shù)太過簡單，局限性太多而且打毛，而且有不少bugs柿赊，所以基于它的版本我做了一些改進，然后封裝成函數(shù)RenameGenesSeurat_v2：
- obj幻枉，Seurat對象
- newnames碰声，新的基因名字集向量，與gene.use的基因集一一對應
- gene.use熬甫，使用的基因集向量胰挑，默認使用所有原始基因，如果給定基因集會對obj取基因子集椿肩，長度需與newnames一樣瞻颂，并且與newnames對應
- de.assay，默認的assay郑象，單細胞數(shù)據(jù)填RNA贡这，空間組數(shù)據(jù)填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)
}
這個函數(shù)不僅僅對單個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)
}
小結與討論
以上函數(shù)適用于已經(jīng)完成聚類分析的Seurat對象，也就是含有pca柜砾，umap等降維信息的對象湃望，如果是什么都沒做，完全就是一個新的只載入矩陣的Seurat對象痰驱，建議還是直接對矩陣操作更方便快捷喜爷。
目前用這個RenameGenesSeurat_v2函數(shù)還沒發(fā)現(xiàn)什么問題，可能還要遺漏需要改基因名的地方萄唇，歡迎評論區(qū)補充檩帐！
