1写在前面
当我们拿到表达矩阵后,需要对ID
进行一个转换,转换为大家可以看懂的gene symbol
。
本期我们介绍一下如何转换,以及其中的一个大坑,线粒体基因!!!😤
2用到的包
rm(list = ls())
library(tidyverse)
library(scater)
library(SingleCellExperiment)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(EnsDb.Hsapiens.v86)
3示例数据
这里我们用一下之前介绍的counts
文件和annotation
文件,然后通过SingleCellExperiment
创建SingleCellExperiment
格式的文件。
3.1 读入数据
counts <- read.delim("./molecules.txt", row.names = 1)
annotation <- read.delim("./annotation.txt", stringsAsFactors = T)
![counts](https://img-blog.csdnimg.cn/img_convert/793a9f57c3a144a5222a7d5948215e34.png)
![annotation](https://img-blog.csdnimg.cn/img_convert/84a1f415b847611054a60c085f404c45.png)
3.2 创建SingleCellExperiment对象
这个dataset
不仅使用了unique molecular identifiers
(UMIs
),还使用了ERCC spike-ins
。 但是!现在大部分droplet
为基础的protocol
已经不使用ERCC
了,只在一些低通量的方法中作为control
使用。
# 注意assays必须是matrix
umi <- SingleCellExperiment(
assays = list(counts = as.matrix(counts)),
colData = annotation)
![alt](https://img-blog.csdnimg.cn/img_convert/434871b34f04983596d210476809eced.png)
我们把ERCC
特征提取出来存到另一个地方,其实也用不到。
altExp(umi,"ERCC") <- umi[grep("^ERCC-",rownames(umi)), ]
umi <- umi[grep("^ERCC-",rownames(umi),invert = T), ]
4ID转Symbol
4.1 ID2SYMBOL
我们现在把scRNA-seq
后注释的ENSEMBL
转为SYMBOL
。
gene_names <- mapIds(org.Hs.eg.db,
keys = rownames(umi), # the keys to select records for from the database.
keytype="ENSEMBL",
columns="SYMBOL", # the columns or kinds of things that can be retrieved from the database.
column="SYMBOL" # the column to search on (for mapIds).
)
4.2 看看多少转换成功了吧
rowData(umi)$SYMBOL <- gene_names
table(is.na(gene_names))
![alt](https://img-blog.csdnimg.cn/img_convert/3ed4ddbbe2d636c8f8df5d4a1a6d887d.png)
4.3 去除转换失败的基因
umi <- umi[! is.na(rowData(umi)$SYMBOL),]
4.4 看一下有没有线粒体基因
很遗憾我们的data
里并没有线粒体基因
,但这显然是不正确的。🤒
grep("^MT-",rowData(umi)$SYMBOL,value = T)
![alt](https://img-blog.csdnimg.cn/img_convert/e323a155ceb270ecb51dcfa718361b3b.png)
4.5 看一下别的线粒体基因
这里我们看下不含MT-
的线粒基因有没有,如ATP8
,又叫MT-ATP8
。
grep("ATP8",rowData(umi)$SYMBOL,value = T)
![alt](https://img-blog.csdnimg.cn/img_convert/22a96b70f453bf9d112fe9f62c1b068d.png)
hhhhhhhh,神奇了,居然是有的!!!!!解决方案往下看吧~
4.6 解决方案
问题就出在org.Hs.eg.db
身上了。😂
这里推荐小伙伴们选择EnsDb.Hsapiens.v86
进行转换。🫠
这里我们可以发现有13个线粒体上的编码基因。
ensdb_genes <- genes(EnsDb.Hsapiens.v86)
MT_names <- ensdb_genes[seqnames(ensdb_genes) == "MT"]$gene_id
is_mito <- rownames(umi) %in% MT_names
table(is_mito)
![alt](https://img-blog.csdnimg.cn/img_convert/b727b747b1e37fbbe0556fa3bcbc5dcf.png)
用EnsDb.Hsapiens.v86
进行ID转换 !~
gene_names <- mapIds(EnsDb.Hsapiens.v86,
keys = rownames(umi), # the keys to select records for from the database.
keytype="GENEID",
column="SYMBOL" # the column to search on (for mapIds).
)
rowData(umi)$SYMBOL <- gene_names
table(is.na(gene_names))
umi <- umi[! is.na(rowData(umi)$SYMBOL),]
我们再看一下是否有线粒体基因吧~
grep("^MT-",rowData(umi)$SYMBOL,value = T)
![alt](https://img-blog.csdnimg.cn/img_convert/40f79d56955c9033940b457910d279a9.png)
Perfect! 13
个线粒体基因。🥰🥳
5线粒体基因很重要吗?
是的,非常重要,在scRNA-seq
中,线粒体基因高表达往往达标细胞状态不佳,在数据分析中应该剔除这类细胞,具体的操作我们在后面的教程中继续分享吧。🫵
![](https://img-blog.csdnimg.cn/img_convert/39d6736961c7c83eb8662411f12aef28.png)
需要示例数据的小伙伴,回复
scRNAseq
获取吧!点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰
本文由 mdnice 多平台发布