驴的单细胞数据基因ID如何转换？

生信技能树

发布于 2024-12-27 19:57:41

2940

大部分人研究的物种都是人、大鼠、小鼠，最近我们的生信入门群里面有个小伙伴学习了单细胞后，开始做自己的课题，物种比较少见，是个驴子：donkey，遇到的一个问题是驴的基因ID转换。

驴的基因ID如下：

head(rownames(ct))
[1] "ENSEASG00005011113" "ENSEASG00005020859" "ENSEASG00005021205" "ENSEASG00005021266" "ENSEASG00005021268"
[6] "ENSEASG00005021261"

学习过我们《转录组测序分析专题》课程的人，肯定一眼就看出来了这个ID来自数据库：Ensembl数据库。

其特征如下：ENS开头+物种缩写+Feature+11位唯一的数字。

那么，我们就需要去这个数据库下载这个物种对应的gtf文件进行ID与Symbol关系提取，而这个小技巧也是我们《转录组测序分析专题》中重点讲过的知识点：

参考基因组注释文件介绍：

每一列的含义：

第九列的具体含义：

课上专门针对这个文件的练习：R版本

rm(list=ls())
## 使用西湖大学的 Bioconductor镜像
# options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
# options("repos"=c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"))
# BiocManager::install("rtracklayer")
# install.packages("tidyverse")
library(tidyverse)
library(rtracklayer)

# 读取gtf文件
gtf <- rtracklayer::import('data/Homo_sapiens.GRCh38.104.chr.gtf.gz') 
class(gtf)
str(gtf)
head(gtf)

# 转为数据框，简单探索
gtf_df <- as.data.frame(gtf)
table(gtf_df$type)
table(gtf_df$gene_biotype)
colnames(gtf_df)

# 过滤基因行
gtf_df <- filter(gtf_df, type=="gene")
head(gtf_df)
as.data.frame(table(gtf_df$gene_biotype))


# 提取 gene id与entrez id
id2name <- gtf_df[, c("gene_id", "gene_name", "gene_biotype")]
head(id2name) # 并不是所有的gene_id都有对应的gene_name

# 去重
id2name <- unique(id2name)
# 保存基因id对应关系
save(id2name, file = "data/id2name.Rdata")
write.table(id2name, file = "data/id2name.xls",row.names = F,sep = "\t",quote = F)


# 也可以提取编码蛋白基因
id2name_mrna <- filter(id2name, gene_biotype=="protein_coding")
head(id2name_mrna)
# 保存基因id对应关系
save(id2name_mrna, file = "data/id2name_mrna.Rdata")
write.table(id2name_mrna, file = "data/id2name_mrna.xls",row.names = F,sep = "\t",quote = F)

bash编程代码版本

zless -S /home/t_rna/database/GRCh38.104/Homo_sapiens.GRCh38.104.chr.gtf.gz | awk -F '\t' '{if($3=="gene") {print$9}}' | awk ' BEGIN{print "gene_id\tgene_name\tgene_biotype"} !/^#/{ \
 id = name = type = "NA";
 for (i=1;i<=NF;i++) { 
   if($i~"gene_id") id=$(i+1)
   if($i~"gene_name") name=$(i+1)
   if($i~"gene_biotype") type=$(i+1)
  }
  print id"\t"name"\t"type } ' | \
  sed 's/[\";]//g' > GRCh38.gene_info.txt

head GRCh38.gene_info.txt

现在到了该融汇贯通的时候了，我们上课讲解的是人这个物种，现在换为驴：

首先先去下载驴的gtf文件：

在这个页面搜索：donkey https://www.ensembl.org/info/about/species.html

进入：https://ftp.ensembl.org/pub/release-113/gtf/equus_asinus/

wget -c https://ftp.ensembl.org/pub/release-113/gtf/equus_asinus/Equus_asinus.ASM1607732v2.113.gtf.gz

提取基因ID与symbol对应关系：

## id 转换
library(tidyverse)
library(rtracklayer)

# 读取gtf文件
gtf <- rtracklayer::import('Equus_asinus.ASM1607732v2.113.gtf.gz') 
# 转为数据框，简单探索
gtf_df <- as.data.frame(gtf)
table(gtf_df$type)
table(gtf_df$gene_biotype)
colnames(gtf_df)

# 过滤基因行
gtf_df <- filter(gtf_df, type=="gene")
head(gtf_df)
as.data.frame(table(gtf_df$gene_biotype))

# 提取 gene id与entrez id
id2name <- gtf_df[, c("gene_id", "gene_name", "gene_biotype")]
id2name <- unique(id2name)
head(id2name) # 并不是所有的gene_id都有对应的gene_name
#id2name <- na.omit(id2name)
loc <- which(is.na(id2name$gene_name))
id2name[loc,2] <- id2name[loc,1]
head(id2name)

对应关系：这里的基因id有很多没有symbol，没有symbol的我使用gene_id进行了替补。

             gene_id          gene_name   gene_biotype
1 ENSEASG00005011113 ENSEASG00005011113 protein_coding
2 ENSEASG00005020859            SLC41A2 protein_coding
3 ENSEASG00005021205 ENSEASG00005021205 protein_coding
4 ENSEASG00005021266               MYF5 protein_coding
5 ENSEASG00005021268               MYF6 protein_coding
6 ENSEASG00005021261              ACSS3 protein_coding

读取单细包数据：

ct <- Read10X('data/',gene.column = 1)
dim(ct)
head(rownames(ct))

comid <- intersect(rownames(ct), id2name$gene_id)
ct_symbol <- ct[comid, ]
ct_symbol[1:4,1:5]
names <- id2name[match(comid, id2name$gene_id), ]
rownames(ct_symbol) <- names$gene_name

# 创建seurat对象
sce.all <- CreateSeuratObject(counts = ct_symbol)
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])

结果如下：