引言

❝
非模式生物单细胞亚群注释参考数据库局限，缺乏完整的基因注释信息和 marker 基因数据库。若基因组注释不完善，可能还需要自己重头组装新的完整的基因组来做后续的分析。
除了去收集已有文献的所提供的 marker 基因外，寻找同源基因是一个不错的选择，同源转换基于蛋白质序列的相似性和一一对应的比对关系，一般情况下，只有一部分的基因符合条件，物种分化时间越远，能保留的基因就越少，从而导致注释结果存在一定的不确定性和误差。
还有一些是基于模式物种已有注释好的单细胞进行标签转移（label transfer）映射到非模式物种的单细胞数据上面（Scanpy 的 scANVI 模型）。
今年 Nature Methods 发表了 SATURN，该模型利用深度学习技术和大型蛋白质语言模型，通过获取 cell embeddings 整合不同物种之间的单细胞转录组测序数据，从而实现跨物种单细胞数据注释。SATURN 通过整合 scRNA-seq 数据的 embeddings 和物种对应的 protein 序列的 embeddings 来构建宏基因空间（macrogene space）。通过不同物种的细胞数据在宏基因空间的相似性，SATURN 可以将不同物种的细胞数据进行聚类，从而解决了以往单纯只能靠有限一一对应的同源基因进行注释的问题。–动植物单细胞测序细胞类型注释还在用同源转换？! 快来试试AI新方法
❝
这里介绍一下使用 OrthoFinder 来寻找同源基因的方法。关于 OrthoFinder 输出文件的解释可以自行搜索查看，下面看看相关一些文章使用的这个方法案例：

安装 OrthoFinder

github 地址：https://github.com/davidemms/OrthoFinder

wget https://github.com/davidemms/OrthoFinder/releases/download/2.5.5/OrthoFinder.tar.gz
tar -zxvf OrthoFinder.tar.gz
测试

❝
OrthoFinder 需要输入蛋白质文件，这里我们下载鹅和人的蛋白质和注释文件作为示例：

Taihu_goose

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/040/182/565/GCF_040182565.1_Taihu_goose_T2T_genome/GCF_040182565.1_Taihu_goose_T2T_genome_protein.faa.gz

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/040/182/565/GCF_040182565.1_Taihu_goose_T2T_genome/GCF_040182565.1_Taihu_goose_T2T_genome_genomic.gtf.gz

human

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_protein.faa.gz

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_genomic.gtf.gz
寻找同源基因

首先需要使用 primary_transcript.py 处理一下冗余的蛋白序列：
python OrthoFinder/primary_transcript.py GCF_000001405.40_GRCh38.p14_protein.faa

Identified as NCBI file

Found 136282 accessions, 43125 genes, 0 unidentified transcripts

Wrote 43125 genes

/mnt/data/home/tycloud/junjun_proj/18.anser_project/primary_transcripts/GCF_000001405.40_GRCh38.p14_protein.faa

python OrthoFinder/primary_transcript.py GCF_040182565.1_Taihu_goose_T2T_genome_protein.faa

Identified as NCBI file

Found 53343 accessions, 22639 genes, 0 unidentified transcripts

Wrote 22639 genes

/mnt/data/home/tycloud/junjun_proj/18.anser_project/primary_transcripts/GCF_040182565.1_Taihu_goose_T2T_genome_protein.faa

然后把处理的文件放在primary_transcripts文件夹下，挂后台运行：
nohup ./OrthoFinder/orthofinder -f ./primary_transcripts -t 40 -a 20&
然后会在当前目录下生成 OrthoFinder/Results_日期/ 这些目录：

后续处理，首先读入注释文件，获取每个基因的 id 和名称对应关系：
setwd(“~/junjun_proj/18.anser_project/primary_transcripts/”)

library(rtracklayer)
library(tidyverse)

human <- import.gff(“../GCF_000001405.40_GRCh38.p14_genomic.gtf.gz”,format = “gtf”) %>%
data.frame()

get protein_id

human_id <- human %>%
select(protein_id,gene) %>% unique() %>%
na.omit()

goose <- import.gff(“../GCF_040182565.1_Taihu_goose_T2T_genome_genomic.gtf”,format = “gtf”) %>%
data.frame()

get protein_id

goose_id <- goose %>%
select(protein_id,gene) %>% unique() %>%
na.omit()
我们需要的文件是 Orthogroups/Orthogroups.tsv：

process orthofinder output

orth_data <- read.delim(“./OrthoFinder/Results_Dec19/Orthogroups/Orthogroups.tsv”) %>%
dplyr::filter(GCF_000001405.40_GRCh38.p14_protein != “”) %>%
dplyr::filter(GCF_040182565.1_Taihu_goose_T2T_genome_protein != “”)

❝
我们给这些蛋白 id 添加对应物种的基因名称：

x = 1

purrr::map_df(1:nrow(orth_data),function(x){
tmp <- orth_data[x,]

# add human symbol
hid <- unlist(strsplit(tmp$GCF_000001405.40_GRCh38.p14_protein,split = “, “)) human_id_tmp <- human_id %>% dplyr::filter(protein_id %in% hid)
human_id_tmp <- human_id_tmp[match(hid,human_id_tmp$protein_id),]

tmp$human_symbol <- paste(human_id_tmp$gene,collapse = “,”)

# add goose symbol
gid <- unlist(strsplit(tmp$GCF_040182565.1_Taihu_goose_T2T_genome_protein,split = “, “)) goose_id_tmp <- goose_id %>% dplyr::filter(protein_id %in% gid)
goose_id_tmp <- goose_id_tmp[match(gid,goose_id_tmp$protein_id),]

tmp$goose_symbol <- paste(goose_id_tmp$gene,collapse = “,”)

# output
return(tmp)
}) -> orth_data_anno

❝
可以看到有些一个蛋白 id 对应多个蛋白 id，也有多个蛋白 id 对应多个蛋白 id 的，有可能这些多个蛋白 id 都隶属与同一个基因，因此我们这里不筛选一对一的 id，我们保留那些多对一/一对多或者多对多那些多的蛋白 id 是同一个基因的行：

get symbol pairs

orth_data_anno_symbol <- orth_data_anno %>%
dplyr::select(human_symbol,goose_symbol)

extract one human symbol match with goose one symbol

x = 1

purrr::map_df(1:nrow(orth_data_anno_symbol),function(x){
tmp <- orth_data_anno_symbol[x,]

# check uniq symbol length
hidn <- length(unique(unlist(strsplit(tmp$human_symbol,split = “,”))))
gidn <- length(unique(unlist(strsplit(tmp$goose_symbol,split = “,”))))

cd <- hidn == 1 & gidn == 1

if(cd){
tmp$human_symbol <- unique(unlist(strsplit(tmp$human_symbol,split = “,”)))
tmp$goose_symbol <- unique(unlist(strsplit(tmp$goose_symbol,split = “,”)))

return(tmp)

}
}) -> uniq_orthologo

save

save(uniq_orthologo,file = “uniq_orthologo.rda”)

length(unique(goose_id$gene))

[1] 18126

length(unique(uniq_orthologo$goose_symbol))

[1] 13320

13320/18126

[1] 0.734856

❝
可以看到鹅的同源基因有大约 73%的在人里面可以找到对应的。 后面可以拿这些同源基因去做人的亚群注释，进而给自己的物种单细胞亚群添加细胞身份。
结尾

❝
路漫漫其修远兮,吾将上下而求索。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。