
jamesbang
2022/10/07阅读:55主题:雁栖湖
🤩 scRNA-seq | 吐血整理的单细胞入门教程(ID转换)(六)
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)


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)

我们把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))

4.3 去除转换失败的基因
umi <- umi[! is.na(rowData(umi)$SYMBOL),]
4.4 看一下有没有线粒体基因
很遗憾我们的data
里并没有线粒体基因
,但这显然是不正确的。🤒
grep("^MT-",rowData(umi)$SYMBOL,value = T)

4.5 看一下别的线粒体基因
这里我们看下不含MT-
的线粒基因有没有,如ATP8
,又叫MT-ATP8
。
grep("ATP8",rowData(umi)$SYMBOL,value = T)

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)

用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)

Perfect! 13
个线粒体基因。🥰🥳
5线粒体基因很重要吗?
是的,非常重要,在scRNA-seq
中,线粒体基因高表达往往达标细胞状态不佳,在数据分析中应该剔除这类细胞,具体的操作我们在后面的教程中继续分享吧。🫵

需要示例数据的小伙伴,回复
scRNAseq
获取吧!点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰
作者介绍

jamesbang
wx🔍: Grassssss 卷起来了