• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    alphafold2培训

    alphafold2培训

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      alphafold2培训

      alphafold2培训

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      生信星球

      • 首页
      • 博客
      • 生信星球
      • ID转换失败,爬起来继续high

      ID转换失败,爬起来继续high

      • 发布者 weinfoeditor
      • 分类 生信星球
      • 日期 2020年3月18日
      测试开头

      【历史文章】那天是生信星球陪你的第501天


         大神一句话,菜鸟跑半年。我不是大神,但我可以缩短你走弯路的半年~

         就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~

         这里有豆豆和花花的学习历程,从新手到进阶,生信路上有你有我!

      0 前言

      拿到一个接近60000 x 6000的原始表达矩阵,行是基因,列是样本。并且行名都是Ensembl ID,还带着版本号。

      我的第一直觉就是先把ID转换成Symbol,保证后续操作的易读性

      下面👇是我对这个原始矩阵的操作

      1 第一步 对原始表达矩阵进行过滤

      很简单粗暴,直接去掉那些在所有样本不表达的基因

      > dim(origin_count)
      [1] 58347  6000
      # 过滤后的矩阵df
      df <- origin_count[rowSums(origin_count)>0,]
      > dim(df)
      [1] 43523  6295

      2 第二步 修改原始行名的Ensembl ID

      这里要先说一下Ensembl 的ID

      2.1 我们常见的Ensembl ID比如:

      ENSG00000279928.1
      ENSG00000279457.2
      详细说明:https://asia.ensembl.org/Help/Faq?id=488

      • 其中,ENS是固定字符,表示它是属于Ensembl数据库的。默认物种是人,如果是小鼠就要用ENSMUS开头,关于物种代码:http://www.ensembl.org/info/genome/stable_ids/index.html

      • G表示:这个ID指向一个基因;E指向Exon;FM指向 protein family;GT指向gene tree;P指向protein;R指向regulary feature;T指向transcript

      • 后面11位数字部分如00000279928 表示基因真正的编号

      • 小数点后不一定每个都有,表示这个ID在数据库中做了几次变更,比如.1做了1次变更,在分析时需要去除

      至于怎么去掉小数点后面的信息也很简单,之前花花介绍过(去掉ensembl id的version信息)

      2.2 小心两个坑

      🤡 第一坑

      不过这里,仅仅是简单去掉还不够,因为我们的目的是让新的Ensembl ID当做行名,但是去掉小数点后会出现一个问题:重复值 ,而行名是不允许重复值出现的

      因此,需要保证新的Ensembl ID中如果出现重复值,那么就只保留第一个

      🤡 第二坑

      之前原始的Ensembl ID是和行名的长度一致的,现在得到了去重复的Ensembl ID,那么由于长度不等,因此不能一一对应到行名上。

      要根据之前Ensembl ID的重复情况,去掉原始数据框对应的行(不过不要担心,重复的值没有多少,所以也去不掉多少行,对结果没有太大影响;并且既然是重复的行,它们的表达量应该也是类似的,所以取其中一个即可)

      2.3 那么整体的操作如下:

      head(rownames(df),3)
      # "ENSG00000000003.14" "ENSG00000000005.5"  "ENSG00000000419.12"

      # 去掉小数点后的gene name(gn)
      gn <- stringr::str_split(rownames(df),"\.",simplify = T)[,1]
      # 根据gn的重复情况,df也进行对应的去除
      df <- df[!duplicated(gn),]
      # 此时df的行数才和gn相等
      rownames(df) <- gn[!duplicated(gn)]

      > dim(df) # 看到相比之前(43523行),过滤掉了39个重复行
      [1] 43484  6295
      > head(rownames(df),3)
      [1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419"

      3 第三步 基因ID转换

      3.1 使用org.Hs.eg.db + mapIds

      suppressMessages(library(org.Hs.eg.db))
      symb <- mapIds(org.Hs.eg.db, keys=rownames(df), 
                     keytype="ENSEMBL", column="SYMBOL")

      > length(rownames(df))
      [1] 43484
      > length(na.omit(symb))
      [1] 22965

      可以看到基本包含了一半的NA值

      3.2 使用org.Hs.eg.db + bitr

      其实结果一样,因为都是使用的org.Hs.eg.db数据库,这里只是列一下方法

      library(clusterProfiler)
      library(org.Hs.eg.db)
      gene_tr <- bitr(rownames(df), fromType = "ENSEMBL",
                      toType = c("SYMBOL"),
                      OrgDb = org.Hs.eg.db)
      # 47.19% of input gene IDs are fail to map...

      > length(na.omit(symb_2$SYMBOL))
      [1] 23196

      可以看到比mapIds能多得到一些基因

      疑问

      这么多的NA基因,难道真的不存在对应的Symbol基因名吗?这个问题值得思考🤔

      很多时候,使用了两种工具,得到差不多的结果,可能就认为事已至此了。但是换种角度思考问题,如果问题不是出在工具呢?如果是数据库呢?

      好,接下来就看看没有对应Symbol的那些Ensembl ID
      not_mapped_ensembl <- names(symb[is.na(symb)])

      > length(not_mapped_ensembl)
      [1] 20519

      > head(not_mapped_ensembl,3)
      [1] "ENSG00000002079" "ENSG00000018607" "ENSG00000020219"

      那么随便挑一个(例如第一个)直接复制、搜索,看看到底有没有基因名

      搜索的结果让人大吃一惊,竟然它真的有基因名,而且还在HGNC、Ensembl、Genecards等数据库有收录
      ID转换失败,爬起来继续high
      image-20191209164711883
      ID转换失败,爬起来继续high
      image-20191209164858412

      看到Ensembl 列出的该基因的信息很全,不由得想到,基因ID转换不止这两种工具,还有一种:Ensembl自家的biomaRt,尝试一下

      3.3 使用hsapiens_gene_ensembl + biomaRt

      看它的帮助文档就很清楚它怎么操作,过程也很简单

      library("biomaRt")
      # 用useMart函数链接到人类的数据库
      ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
      # 除以以外还有很多:使用 listDatasets(ensembl) 查看

      # 依然是要列出查询的gene ID,以及要转换的格式
      # 【biomaRt强大之处就在于它的转换格式种类非常之丰富】
      attributes <- listAttributes(ensembl)
      View(attributes) # 查看转换格式

      # 列出要转换的ID(比如前面匹配失败的3个基因)
      value <- head(not_mapped_ensembl,3)
      # 列出当前格式和要比对的格式(当前是ensembl_gene_id;我要比对hgnc_symbol和chromosome_name)
      attr <- c("ensembl_gene_id","hgnc_symbol","chromosome_name")

      # 最后运行就好了(其中filters指定当前基因 ID类型)
      ids <- getBM(attributes = attr,
                   filters = "ensembl_gene_id",
                   values = value,
                   mart = ensembl)

      最后看一下结果,而且这个结果扩展性极强:

      > ids
        ensembl_gene_id  hgnc_symbol chromosome_name
      1 ENSG00000002079        MYH16               7
      2 ENSG00000018607       ZNF806               2
      3 ENSG00000020219      CCT8L1P               7
      biomaRt的小问题

      当提供的基因数量太多时,容易断线;它的检索速度非常慢!【不过用来应急还是可以的】

      心得

      当使用某个或某几个工具得到的结果感觉不满意时,不要只考虑工具的优劣,还要思考它们在后台帮我们做了什么。

      比如这里基因ID的转换,靠的其实就是背后的数据库,工具虽有差异但没那么大(比如bitr比mapIds多了几个基因)。有时一条路走不通,换个工具“背后的”数据库未尝不是个好的解决办法!


      点击底部的“阅读原文”,获得更好的阅读体验哦😻

      初学生信,很荣幸带你迈出第一步

      🤓生信星球 🌎~ 一个不拽术语、通俗易懂的生信知识平台

      测试结尾

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      作者头像
      weinfoeditor

      上一篇文章

      新手如何理解管道符号
      2020年3月18日

      下一篇文章

      maftools|TCGA肿瘤突变数据的汇总,分析和可视化
      2020年3月18日

      你可能也喜欢

      8-1651673488
      生信零基础入门学习小组长期报名中(2022仍继续)
      7 4月, 2022
      2-1651673738
      简化版的ROC曲线
      21 2月, 2022
      8-1651674718
      支持向量机模型
      19 11月, 2021

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      白介素-17受体信号的自主激活,维持炎症并促进疾病进展
      048月2023
      MCT4依赖的乳酸分泌抑制LKB1缺陷肺腺癌的抗肿瘤免疫
      187月2023
      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?