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

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中，线粒体基因高表达往往达标细胞状态不佳，在数据分析中应该剔除这类细胞，具体的操作我们在后面的教程中继续分享吧。🫵

最后祝大家早日不卷!~

scRNA-seq

🤩 scRNA-seq | 吐血整理的单细胞入门教程（ID转换）（六）

1写在前面

2用到的包