前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >🤩 scRNA-seq | 吐血整理的单细胞入门教程(ID转换)(六)

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

作者头像
生信漫卷
发布2022-10-31 17:20:37
1.1K0
发布2022-10-31 17:20:37
举报
文章被收录于专栏:R语言及实用科研软件

1写在前面

当我们拿到表达矩阵后,需要对ID进行一个转换,转换为大家可以看懂的gene symbol。 本期我们介绍一下如何转换,以及其中的一个大坑线粒体基因!!!😤

2用到的包

代码语言:javascript
复制
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 读入数据

代码语言:javascript
复制
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使用。

代码语言:javascript
复制
# 注意assays必须是matrix
umi <- SingleCellExperiment(
  assays = list(counts = as.matrix(counts)),
  colData = annotation)

我们把ERCC特征提取出来存到另一个地方,其实也用不到。

代码语言:javascript
复制
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

代码语言:javascript
复制
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 看看多少转换成功了吧

代码语言:javascript
复制
rowData(umi)$SYMBOL <- gene_names
table(is.na(gene_names))

4.3 去除转换失败的基因

代码语言:javascript
复制
umi <- umi[! is.na(rowData(umi)$SYMBOL),]

4.4 看一下有没有线粒体基因

很遗憾我们的data里并没有线粒体基因,但这显然是不正确的。🤒

代码语言:javascript
复制
grep("^MT-",rowData(umi)$SYMBOL,value = T)

4.5 看一下别的线粒体基因

这里我们看下不含MT-的线粒基因有没有,如ATP8,又叫MT-ATP8

代码语言:javascript
复制
grep("ATP8",rowData(umi)$SYMBOL,value = T)

hhhhhhhh,神奇了,居然是有的!!!!!解决方案往下看吧~


4.6 解决方案

问题就出在org.Hs.eg.db身上了。😂 这里推荐小伙伴们选择EnsDb.Hsapiens.v86进行转换。🫠 这里我们可以发现有13线粒体上的编码基因。

代码语言:javascript
复制
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转换 !~

代码语言:javascript
复制
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),]

我们再看一下是否有线粒体基因吧~

代码语言:javascript
复制
grep("^MT-",rowData(umi)$SYMBOL,value = T)

Perfect! 13线粒体基因。🥰🥳

5线粒体基因很重要吗?

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


最后祝大家早日不卷!~


本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2022-10-07,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信漫卷 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1写在前面
  • 2用到的包
  • 3示例数据
    • 3.1 读入数据
      • 3.2 创建SingleCellExperiment对象
      • 4ID转Symbol
        • 4.1 ID2SYMBOL
          • 4.2 看看多少转换成功了吧
            • 4.3 去除转换失败的基因
              • 4.4 看一下有没有线粒体基因
                • 4.5 看一下别的线粒体基因
                  • 4.6 解决方案
                  • 5线粒体基因很重要吗?
                  领券
                  问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档