2022-07-01 | 教你实战单细胞信息注释函数：R语言数据处理中的替换操作（附常用示例gsub函数）

适用背景

在R语言中，我们需要对字符串、向量和数据框等数据类型进行替换操作，有时候是因为需要更换别名，有时候是因为数据存在错误需要修正，有时候则是因为需要删除某些信息。本文将介绍常用的替换函数gsub的常用用法，但gsub也存在某些局限性，一般只能进行一次指定情况的操作。例如在单细胞数据分析的信息注释过程中，我们常常需要把无监督聚类得到的clusters注释成细胞类型，如果每一个clusters都写一行替换的代码就会显得相当冗余，因此可以封装成一个函数进行类似的处理就会简单一些。因此，本文后半部分将介绍批量替换写成函数的方法。

gsub函数

R语言中，最常用的替换函数是gsub，其用法也比较容易理解，一般只需传入三个参数，gsub(匹配内容，替换内容，操作对象）：

string='strings'
gsub('s','S',string)
[1] "StringS"

• 删掉某些字符串:

string='strings'
gsub('s','',string)
[1] "tring"

• 如果只想替换开头的s，则需要使用正则表达式，代码如下：

string='strings'
gsub('^s','S',string)
[1] "Strings"

• 只替换结尾的s则写成:

gsub('s$','S',string)

• 如果有多个匹配项，则使用管道符 |

string='strings'
gsub('^s|i','X',string)
[1] "XtrXngs"

• 把s及s后面的一个字符替换掉：

string='strings'
gsub('s..','X',string)
[1] "Xings"

• 把i及后面的字符串替换掉：

string='strings'
gsub('i.*','X',string)
[1] "strX"

• 如果字符串里有特殊的字符，而匹配内容又含有这种特殊字符，例如^，则用方括号[]括起来：

#不用方括号
string='^strings'
gsub('^','X',string)
[1] "X^strings"
#使用方括号
gsub("[.^]",'X',string)
[1] "Xstrings"
#不用方括号
> gsub(".^",'X',string)
[1] "^strings"

• 模糊匹配数字后替换：

string='0strings1'
gsub('[0-9]','X',string)
[1] "XstringsX"

• 模糊匹配字母后替换

string='0strings1'
gsub('[a-z]','X',string)
[1] "0XXXXXXX1"

单细胞注释简单方法

以Celltype列为例，以下all为Seurat对象，也可以是数据框。

• 可以用which函数，返回的是逻辑值

all$Celltype[which(all$seurat_clusters == "24")]<-"XCL+ NK"

• 用grep函数，返回索引值

all$Celltype[grep("24",all$seurat_clusters)]<-"XCL+ NK"

• 用grepl函数，与grep类似，grepl返回的是逻辑值，但效果与grep一样

all$Celltype[grepl("24",all$seurat_clusters)]<-"XCL+ NK"

以上方法都比较好理解，但是如果clusters比较多，写起来就很冗余，有时候clusters多达几十种，写起来不美观，可以写个函数封装一下。

单细胞注释函数

其实单细胞数据注释也可以看成是一种对数据框的替换操作，因此可以提取Seurat对象的meta.data信息运行以下函数。

函数一 anno_cell

参数简介：

• meta，输入的数据框，例如obj@meta.data
• anno.list，注释对应信息，例如anno.list=list("OLI"=c(1,13,11,3,23,18), "MIC"=c(6,22))
• ref.col，原始列的列名，默认是无监督聚类得到的聚类列seurat_clusters
• anno.col，注释列的列名，默认是Celltype，如果此列不存在则会自动建一列并标记为'unknown'，如果存在，则此列之前的信息会保留，但可能会被覆盖

anno_cell <- function(meta,anno.list=NULL,ref.col='seurat_clusters',anno.col='Celltype'){
colist <- unique(colnames(meta))
if (length(grep(anno.col,colist))==0) {
meta[,anno.col] <- 'unknown'
}
nlen <- length(anno.list)
for (i in 1:nlen) {
name <- names(anno.list)[i]
clust <- anno.list[[name]]
c1 <- which(meta[,ref.col] %in% clust)
if (length(c1)!=0) {
meta[,anno.col][c1] = name
}
}
return(meta)
}

使用示例（obj是Seurat对象）：

olic <- c(1,13,11,3,23,18)
micc <- c(6,22)
vlmc <- c(37,40)
nfol <- c(38)
opcc <- c(2)
unknownc <- c(39)
inc <- c(31,25,16,29)
exc <- c(14,12,7,32,34,28,27,24,4,0,35,8,10,30,33,21,36,19)
ast <- c(17,20,41,9,5,15,26)
obj@meta.data <- anno_cell(obj@meta.data,anno.list=list("OLI"=olic,
                                      "MIC"=micc,
                                      "ENDO"=vlmc,
                                      "NFOL"=nfol,
                                      "OPC"=opcc,
                                      "IN"=inc,
                                      "EX"=exc,
                                      "AST"=ast)

以上就完成注释了，可以通过Dimplot函数检查注释得是否正确

函数二 map_anno

有时候，注释结果不是根据无监督聚类得到的，而是需要从其它方法，或者亚群聚类得到信息来注释，简单来说就是需要把从数据框df1中信息转移到数据框df2中，但df1和df2的维度不相等，因此需要进行匹配注释。当然，Seurat自带AddMetaData函数，但由于维度问题，可能会报错。
参数简介如下：

• ref，参考数据框，含有注释信息
• que，待注释数据框
• ref.col，默认为Celltype，参考数据框注释信息所在列
• que.col，默认为Celltype，待注释的列名，如果此列不存在则会自动建一列并标记为'unknown'，如果存在，则此列之前的信息会保留，但可能会被覆盖

map_anno <- function(ref,que,ref.col='Celltype',que.col='Celltype'){
colist <- unique(colnames(que))
if (length(grep(que.col,colist))==0) {
que[,que.col] <- 'unknown'
}
c1 <- which(rownames(ref) %in% rownames(que))
if (length(c1)!=dim(ref)[1]) {
print('Note: some cells in ref are not in que!')
ref <- ref[c1,]
}


nlen <- length(unique(ref[,ref.col]))
for (i in 1:nlen) {
name <- unique(ref[,ref.col])[i]
c1 <- which(ref[,ref.col] %in% name)
cellist <- rownames(ref)[c1]
que[cellist,que.col] <- name
}
return(que)
}

使用示例：

df1 <- read.table('cell_info.xls',sep="\t",header=T,stringsAsFactor=F)
obj@meta.data <- map_anno(ref=df1,que=obj@meta.data)

更新函数三 seq_anno

最近拿到师兄给的注释表格是把seurat_clusters从0群到n群的注释表格，如下图：

因此上面的两个函数都不太方便，还需要进一步整理才能用上前面的两个函数，这个跟每个人的细胞注释信息习惯有关，因此写了这个新函数seq_anno，可以按seurat_clusters列的顺序进行注释。

• meta，数据框
• ref.col，参考列，默认为seurat_clusters，可以改成其它分辨率的结果，例如RNA_snn_res.1.5
• anno.col，注释列，默认为Celltype
• anno.list，从0群到n群注释信息的向量，例如c("BC","MC","TC")，其中0群为BC，1群为MC，2群为TC。

seq_anno <- function(meta,ref.col='seurat_clusters', anno.col='Celltype',anno.list=NULL) {
colist <- unique(colnames(meta))
meta$ref.col <- meta[,ref.col]
if (length(grep(anno.col,colist))==0) {
meta[,anno.col] <- 'unknown'
}
nmax <- max(meta$ref.col)
for (i in 0:nmax){
meta[,anno.col][grep(i,meta$ref.col)] <- anno.list[i+1]
}
return(meta)
}

使用示例：

meta为数据框，cell.anno为向量

> head(cell.anno)
[1] "CD16+Mono" "CD14+Mono" "CD16+Mono" "CD16+Mono" "CD1C+DC"   "CD1C+DC"
> meta <- seq_anno(meta,anno.list=cell.anno,ref.col='seurat_clusters', anno.col='Celltype')

可以简写为

df2 <- seq_anno(df1,anno.list=cell.anno)

映射关系替换

基因ID与基因名之间的替换是分析过程中容易遇到的问题，简单来说可以写个循环，但是基因有上万个时就比较麻烦了，不过已经有人写了现成的函数。

目标：df2是基因ID与基因名的映射关系，需要把df1的gene列的基因ID替换成基因名。

df1数据结构：

                      p_val avg_logFC pct.1 pct.2    p_val_adj cluster
AMEX60DD035333 1.276051e-76 0.2897204 0.811 0.836 3.139595e-72       0
AMEX60DD000806 4.985819e-61 0.5910303 0.739 0.894 1.226711e-56       0
AMEX60DD039553 1.624536e-44 0.2744919 0.175 0.469 3.997009e-40       0
AMEX60DD025106 1.054587e-41 0.2832646 0.145 0.396 2.594705e-37       0
AMEX60DD008073 2.048906e-40 0.2503263 0.136 0.375 5.041127e-36       0
AMEX60DD002836 4.554993e-40 0.2639632 0.168 0.430 1.120710e-35       0
                         gene       FC
AMEX60DD035333 AMEX60DD035333 1.336054
AMEX60DD000806 AMEX60DD000806 1.805848
AMEX60DD039553 AMEX60DD039553 1.315862
AMEX60DD025106 AMEX60DD025106 1.327456
AMEX60DD008073 AMEX60DD008073 1.284444
AMEX60DD002836 AMEX60DD002836 1.302080

df2数据结构：

      Axolotl_ID  hs_gene
1 AMEX60DD000002   ZNF569
2 AMEX60DD000006   ZNF141
3 AMEX60DD000016    FZD10
4 AMEX60DD000033   GLT1D1

• 函数一 stri_replace_all_regex

library(stringi)
df3 <- stri_replace_all_regex(df1$gene,df2$Axolotl_ID,df2$hs_gene,vectorize=F)

• 函数二 mgsub

library(mgsub)
df4 <- mgsub(df1$gene,df2$Axolotl_ID,df2$hs_gene)

两个函数用法一样，输入三个参数：待替换向量，映射关系的原始ID，映射关系的替换内容。实测stri_replace_all_regex比mgsub快，推荐使用stri_replace_all_regex。

---

根据映射关系替换数据框某一列

map_index <- function(df,index=NULL,content=NULL,mapcol=NULL,newcol='newcol'){
df1 <- data.frame(content)
rownames(df1) <- index
df2 <- df1[df[,mapcol],]
df[,'newcol'] <- df2
return(df)
}

• 使用示例

df3 <- map_index(df2,index=celt,content=major,mapcol='Celltype',newcol='newcol')

使用到的数据格式如下：

> head(df2)
        orig.ident nCount_Spatial nFeature_Spatial percent.mt nCount_SCT
BIN.214    Spatial           3096             1030 0.38759690      11238
BIN.371    Spatial           5609             1716 0.44571225      11438
BIN.374    Spatial          11045             2608 0.09959258      11533
BIN.375    Spatial           9026             2246 0.33237314      11646
BIN.378    Spatial           9520             2412 0.23109244      11541
BIN.379    Spatial          12353             3047 0.09714239      12151
        nFeature_SCT SCT_snn_res.0.5 seurat_clusters BIN.100 bin100.y bin100.x
BIN.214         1814               0               0 BIN.214        2       52
BIN.371         1853               4               4 BIN.371        3       47
BIN.374         2608               4               4 BIN.374        3       50
BIN.375         2248               4               4 BIN.375        3       51
BIN.378         2412               4               4 BIN.378        3       54
BIN.379         3047               2               2 BIN.379        3       55
        refined_pred ref.col major  sub Celltype Major
BIN.214           10      10   mDC  mDC      mDC    DC
BIN.371            5       5    LC LC_2       LC    LC
BIN.374            5       5    LC LC_2       LC    LC
BIN.375            5       5    LC LC_2       LC    LC
BIN.378            5       5    LC LC_2       LC    LC
BIN.379            5       5    LC LC_2       LC    LC
> celt
 [1] "CP"            "DMC"           "DVR"           "EG"
 [5] "LC"            "lDC"           "MC"            "mDC"
 [9] "PT"            "Septum_l"      "Septum_m"      "Striatum"
[13] "Striatum_CCK"  "Striatum_TAC4" "VLMC"
> major
 [1] "CP"     "HIP"    "DVR"    "VZ"     "LC"     "DC"     "HIP"    "DC"
 [9] "PT"     "Septum" "Septum" "STR"    "STR"    "STR"    "VLMC"

注意index和content是对应的映射关系向量，顺序不能乱，且index不能有重复的元素。

---

小结与讨论

替换操作不只有gsub函数这一种方法，我写的2个函数就没有使用到gsub，当然也可以写成gsub的脚本，但已经写好了就懒得写了。