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

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

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

驴的基因ID如下：

代码语言：javascript代码运行次数：0运行复制

head(rownames(ct))
[1] "ESEASG0000501111" "ESEASG00005020859" "ESEASG00005021205" "ESEASG00005021266" "ESEASG00005021268"
[6] "ESEASG00005021261"

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

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

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

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

每一列的含义：

第九列的具体含义：

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

代码语言：javascript代码运行次数：0运行复制

rm(list=ls())
## 使用西湖大学的 Bioconductor镜像
# opti(BioC_mirror=";)
# opti("repos"=c(CRA="/"))
# BiocManager::install("rtracklayer")
# install.packages("tidyverse")
library(tidyverse)
library(rtracklayer)

# 读取gtf文件
gtf <- rtracklayer::import('data/Homo_sapiens.GRCh8.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")
(id2name, file = "data/id2name.xls", = 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")
(id2name_mrna, file = "data/id2name_mrna.xls", = F,sep = "\t",quote = F)

bash编程代码版本

代码语言：javascript代码运行次数：0运行复制

zless -S /home/t_rna/database/GRCh8.104/Homo_sapiens.GRCh8.gtf.gz | awk -F '\t' '{if($=="gene") {print$9}}' | awk ' BEGI{print "gene_id\tgene_name\tgene_biotype"} !/^#/{ \
 id = name = type = "A";
 for (i=1;i<=F;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' > GRCh8.gene_

head GRCh8.gene_

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

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

在这个页面搜索：donkey .html

进入：/

代码语言：javascript代码运行次数：0运行复制

wget -c .ASM160772v2.11.gtf.gz

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

代码语言：javascript代码运行次数：0运行复制

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

# 读取gtf文件
gtf <- rtracklayer::import('Equus_asinus.ASM160772v2.11.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 <- (id2name)
loc <- which((id2name$gene_name))
id2name[loc,2] <- id2name[loc,1]
head(id2name)

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

代码语言：javascript代码运行次数：0运行复制

             gene_id          gene_name   gene_biotype
1 ESEASG0000501111 ESEASG0000501111 protein_coding
2 ESEASG00005020859            SLC41A2 protein_coding
 ESEASG00005021205 ESEASG00005021205 protein_coding
4 ESEASG00005021266               MYF5 protein_coding
5 ESEASG00005021268               MYF6 protein_coding
6 ESEASG00005021261              ACSS protein_coding

读取单细包数据：

代码语言：javascript代码运行次数：0运行复制

ct <- Read10X('data/', = 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$RA$counts[1:10, 1:2])

结果如下：

后面学员应该就可以进行愉快的分析了~

本文参与 腾讯云自媒体同步曝光计划，分享自。原始发表：2024-12-27，如有侵权请联系 cloudcommunity@tencent 删除data编程入门数据数据库