3种方法为你的稀有物种建立生物学知识数据库

生信技能树

发布于 2022-12-16 14:39:11

1.1K1

发布于 2022-12-16 14:39:11

文章被收录于专栏：生信技能树

一般来说，生命科学领域的数据分析方向，除了人和鼠之外的物种，都是稀有物种，大家很容易从目前有的测序数据的分布可以看出来，这里说的鼠是小鼠，并不是大鼠哈，如果是其他物种，比如拟南芥，水稻，斑马鱼，可能是稍微规范一点的稀有物种，猪马牛羊狗，甚至原核生物就更是难上加难了。

这里，我们以大鼠为例，演示3种方法为你的稀有物种建立生物学知识数据库，以单细胞数据分析环节的细胞通讯，CellChat为例，可以在：CellChat学习笔记【一】——通讯网络构建了解它的基础用法，**CellChat** 有一个专门的数据库，叫做CellChatDB，这个数据库是 CellChat 的作者们通过阅读大量文献，手动整理出来的“受体-配体”对，目前有人、鼠以及斑马鱼的版本。其中

人的叫做 CellChatDB.human，
鼠的叫做 CellChatDB.mouse，
斑马鱼的叫做 CellChatDB.zebrafish。

这样的话，如果大家做的是稀有物种的单细胞，想使用CellChat来做细胞通讯，就有点麻烦了。比如大鼠，我们演示3种方法为你的稀有物种建立生物学知识数据库。

首先看看CellChat自带的3个数据库

每个数据库文件在r里面的都是一个简单的list，所以很容易探索：

library(CellChat)


dbh <- CellChatDB.human 
lapply(dbh, dim)
dbm <- CellChatDB.mouse
lapply(dbm, dim)
dbz <- CellChatDB.zebrafish
lapply(dbz, dim)

cbind(
  do.call(rbind,lapply(dbh, dim)) ,
  do.call(rbind,lapply(dbm, dim)),
  do.call(rbind,lapply(dbz, dim))
)

如下所示：

         [,1] [,2]  [,3] [,4]  [,5] [,6]
interaction  1939   11  2019   11  2774   11
complex       157    4   155    4   267    2
cofactor       31   16    32   16    62    5
geneInfo    41787    6 45551    6 15778    2

可以看到人和鼠基本上没有差异，但是斑马鱼的数据库记录就有点迷惑了，它明明是基因数量少了一半，但是通讯数据库记录仍然是差不多，可能是因为细胞通讯是比较保守的行为，在多个物种都是可以互相转化的。

我们也可以肉眼去看看数据库里面的记录到底是什么：

dbh_i = dbh$interaction
tail(sort(table(dbh_i$pathway_name)))
dbh_p = dbh$complex
dbh_f = dbh$cofactor

dbm_i = dbm$interaction
tail(sort(table(dbm_i$pathway_name)))
dbm_p = dbm$complex
dbm_f = dbm$cofactor

dbz_i = dbz$interaction
tail(sort(table(dbz_i$pathway_name)))
dbz_p = dbz$complex
dbz_f = dbz$cofactor

可以看到，这个时候人和鼠基本上没有区别：

> tail(sort(table(dbh_i$pathway_name)))

     FGF      CCL      BMP  LAMININ COLLAGEN      WNT 
      60       66       76      143      187      320 

> tail(sort(table(dbm_i$pathway_name)))

     CCL      BMP    MHC-I  LAMININ COLLAGEN      WNT 
      68       76       88      143      198      320 

> tail(sort(table(dbz_i$pathway_name)))

     FGF    TIGIT  LAMININ COLLAGEN      BMP      WNT 
      81      104      144      170      214      504

就是斑马鱼的记录反而是多了一些。。。

那，我们选取其中一个肉眼看看吧，比如 SOMATOSTATIN，首先我们看看人和鼠有什么区别：


> dbm_i[dbm_i$pathway_name=='SOMATOSTATIN',1:4]
           interaction_name pathway_name ligand receptor
SST_SSTR1         SST_SSTR1 SOMATOSTATIN    Sst    Sstr1
SST_SSTR2         SST_SSTR2 SOMATOSTATIN    Sst    Sstr2
SST_SSTR3         SST_SSTR3 SOMATOSTATIN    Sst    Sstr3
SST_SSTR4         SST_SSTR4 SOMATOSTATIN    Sst    Sstr4
SST_SSTR5         SST_SSTR5 SOMATOSTATIN    Sst    Sstr5
CORT_SSTR1       CORT_SSTR1 SOMATOSTATIN   Cort    Sstr1
CORT_SSTR2       CORT_SSTR2 SOMATOSTATIN   Cort    Sstr2
CORT_SSTR3       CORT_SSTR3 SOMATOSTATIN   Cort    Sstr3
CORT_SSTR4       CORT_SSTR4 SOMATOSTATIN   Cort    Sstr4
CORT_SSTR5       CORT_SSTR5 SOMATOSTATIN   Cort    Sstr5

> dbm_i[dbm_i$pathway_name=='SOMATOSTATIN',1:4]
           interaction_name pathway_name ligand receptor
SST_SSTR1         SST_SSTR1 SOMATOSTATIN    Sst    Sstr1
SST_SSTR2         SST_SSTR2 SOMATOSTATIN    Sst    Sstr2
SST_SSTR3         SST_SSTR3 SOMATOSTATIN    Sst    Sstr3
SST_SSTR4         SST_SSTR4 SOMATOSTATIN    Sst    Sstr4
SST_SSTR5         SST_SSTR5 SOMATOSTATIN    Sst    Sstr5
CORT_SSTR1       CORT_SSTR1 SOMATOSTATIN   Cort    Sstr1
CORT_SSTR2       CORT_SSTR2 SOMATOSTATIN   Cort    Sstr2
CORT_SSTR3       CORT_SSTR3 SOMATOSTATIN   Cort    Sstr3
CORT_SSTR4       CORT_SSTR4 SOMATOSTATIN   Cort    Sstr4
CORT_SSTR5       CORT_SSTR5 SOMATOSTATIN   Cort    Sstr5

如我们所猜测的，其实人和鼠仅仅是基因名字大小写差异而已。。。。

反而是斑马鱼，我肉眼看了看，基因名字太诡异，这里不予考虑吧，听说是果蝇基因名字也很诡异，这样的很多上古时代的诺贝尔奖生物学家所在的领域的知识体系架构都及其的乱。。。

既然我们认识了数据库知识架构，现在来看我们的3个方案：

方案1：直接把你的大鼠单细胞表达量矩阵的基因名字修改为人类即可

前面我们提到的CellChat的对象构建，超级简单：

library(CellChat)
library(Seurat)
library(SeuratData)
data("pbmc3k.final")
ct <- GetAssayData(object = pbmc3k.final, slot = 'data')
meta <- pbmc3k.final@meta.data
# 这个 ct 是一个稀疏矩阵 
cellchat <- createCellChat(object = ct,
                           meta = meta,
                           group.by = 'seurat_annotations')

如果是大鼠单细胞表达量矩阵，有两个方案，首先是简单粗暴修改为大写基因就是人类基因名字啦；

rownames(ct) = toupper(rownames(ct) )

后续的数据分析直接当做是人类单细胞表达量矩阵即可，反正我们也不知道人类和小鼠大鼠的细胞通讯到底差异大不大，数据分析仅仅是启发，并不需要百分比精确。

另外一个方案是进行同源基因转化，比如 babelgene和homologene ，我们这里针对 CellChat里面的基因进行举例：

# install.packages("babelgene")
library(babelgene)
cg = unique(CellChatDB.human$interaction$receptor)
length(cg)
tmp = orthologs(genes =  cg ,
          species = "mouse")
tmp$new = toupper(tmp$symbol) 
tmp[tmp$human_symbol != tmp$new,c(1,5)]

可以看到确实是简单粗暴修改为大写基因会有误差；

> tmp[tmp$human_symbol != tmp$new,c(1,5)]
    human_symbol    symbol
13         AGTR1    Agtr1a
42         CD209    Cd209b
45         CD244    Cd244a
56          CD55     Cd55b
61          CD8B     Cd8b1
62         CD8B2     Cd8b1

但是有误差的基因占比少得可怜。

另外一个方法：


# install.packages('homologene')#R
library(homologene)
tmp <- human2mouse( cg )
tmp <- na.omit(tmp) #去除NA值
dim(tmp)
tmp$new = toupper(tmp$mouseGene) 
tmp[tmp$humanGene != tmp$new,]

也是可以看到误差：

> tmp[tmp$humanGene != tmp$new,]
    humanGene mouseGene humanID   mouseID       new
53       IL4R     Il4ra    3566     16190     IL4RA
71  TNFRSF10A Tnfrsf10b    8797     21933 TNFRSF10B
91     LILRB3      Pirb   11025     18733      PIRB
92     LILRB3   Gm14548   11025 100038909   GM14548
93     LILRB3     Pira2   11025     18725     PIRA2
94     LILRB3     Pira1   11025     18722     PIRA1
95     LILRB3   Gm15448   11025 100041146   GM15448
96     LILRB3    Lilra6   11025     18726    LILRA6
97     LILRB3     Pira6   11025     18729     PIRA6
109     AGTR1    Agtr1a     185     11607    AGTR1A
247      CD55      Daf2    1604     13137      DAF2

就很诡异， CD55 和 Daf2 ，看起来是八竿子打不着的啊，搞小鼠科学研究的那些人就不能低头吗，给人类科学让路，自己那一小撮人开会自己检讨一下自己的基因命名问题，以人类研究为准，身为科学家居然搞不清楚自己搞科学的目的，不应该是造福人类吗？