jamesbang

V1

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)
counts
counts
annotation
annotation

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
V1

wx🔍: Grassssss 卷起来了