• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

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

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

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

      未分类

      • 首页
      • 博客
      • 未分类
      • 用R提取代表转录本

      用R提取代表转录本

      • 发布者 weinfoadmin
      • 分类 未分类, 老俊俊的生信笔记
      • 日期 2021年9月10日
      • 评论 0评论

      感谢老俊俊的大力支持。我们会每日跟新,欢迎您关注老俊俊的生信笔记。

      用 R 提取代表性转录本

      前言:

      上上篇文章讲了从 GTF 文件的提取代表性的转录本 gtf 信息:提取代表转录本之 gencode,如果想获得代表性转录本序列的话,可以使用 gffread 软件取获得 cDNA 序列。

      也可以直接下载转录组序列,但是里面包含了同一个基因的多个转录本,我们可以筛选出代表性的转录本。筛选原则与上面的文章有一些不同,我们先看下面这张图:

      基于 gtf 注释文件挑选的最长转录本的长度是包含了内含子的长度,这样选择不能真实反应转录本的实际长度,在选择 CDS 最长的前提下,再选择转录本外显子的总长度最长更加合理。

      提取实操:

      1、数据准备

      首先去 gencode 数据库下载蛋白编码基因的转录组序列文件:

      # 下载
      $ wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.pc_transcripts.fa.gz
      # 查看有多少个转录本序列
      $ zless gencode.v38.pc_transcripts.fa.gz | grep '>' | wc -l
      106143
      # 简单查看前面几行
      $ zless gencode.v38.pc_transcripts.fa.gz | head
      >ENST00000641515.2|ENSG00000186092.7|OTTHUMG00000001094.4|OTTHUMT00000003223.4|OR4F5-201|OR4F5|2618|UTR5:1-60|CDS:61-1041|UTR3:1042-2618|
      CCCAGATCTCTTCAGTTTTTATGCCTCATTCTGTGAAAATTGCTGTAGTCTCTTCCAGTT
      ATGAAGAAGGTAACTGCAGAGGCTATTTCCTGGAATGAATCAACGAGTGAAACGAATAAC
      TCTATGGTGACTGAATTCATTTTTCTGGGTCTCTCTGATTCTCAGGAACTCCAGACCTTC
      CTATTTATGTTGTTTTTTGTATTCTATGGAGGAATCGTGTTTGGAAACCTTCTTATTGTC
      ATAACAGTGGTATCTGACTCCCACCTTCACTCTCCCATGTACTTCCTGCTAGCCAACCTC
      TCACTCATTGATCTGTCTCTGTCTTCAGTCACAGCCCCCAAGATGATTACTGACTTTTTC
      AGCCAGCGCAAAGTCATCTCTTTCAAGGGCTGCCTTGTTCAGATATTTCTCCTTCACTTC
      TTTGGTGGGAGTGAGATGGTGATCCTCATAGCCATGGGCTTTGACAGATATATAGCAATA
      TGCAAGCCCCTACACTACACTACAATTATGTGTGGCAACGCATGTGTCGGCATTATGGCT

      ENST00000641515.2|ENSG00000186092.7|OTTHUMG00000001094.4|OTTHUMT00000003223.4|OR4F5-201|OR4F5|2618|UTR5:1-60|CDS:61-1041|UTR3:1042-2618|

      每个序列的 ID 由这么一长串内容组成,分别是:

      • ENST00000641515.2 : ensembl 的转录本 id
      • ENSG00000186092.7 : ensembl 的基因 id
      • OTTHUMG00000001094.4 : havana 的基因 id
      • OTTHUMT00000003223.4 : havana 的转录本 id
      • OR4F5-201 : 转录本名称
      • OR4F5 : 基因名
      • 2618 : 转录本序列总长度
      • UTR5:1-60 : 转录本 5’UTR 的位置
      • CDS:61-1041 : 转录本 CDS 的位置
      • UTR3:1042-2618 : 转录本 3’UTR 的位置

      我还想重新命名每个转录本序列的 ID,只保留基因名、gene_id、transcript_id、transcript_name、CDS 起始位置和终止位置、转录本长度,类似于下面这种更加简洁一点:

      >OR4F5_ENSG00000186092.7_ENST00000641515.2_OR4F5-201_61_1041_2618
      >OR4F29_ENSG00000284733.2_ENST00000426406.4_OR4F29-201_1_939_939
      >OR4F16_ENSG00000284662.1_ENST00000332831.4_OR4F16-201_20_958_995

      把转录组序列名字输出保存到一个文件里,同时去除 havana 的基因 id 和转录本 id:

      # 按 | 字符分割
      $ zless gencode.v38.pc_transcripts.fa.gz |grep '>' |sed 's/|/t/g' |cut -f 1,2,5,6,7,8,9,10 > name.txt
      # 查看内容
      $ head -3 name.txt
      >ENST00000641515.2      ENSG00000186092.7       OR4F5-201       OR4F5   2618    UTR5:1-60       CDS:61-1041     UTR3:1042-2618
      >ENST00000426406.4      ENSG00000284733.2       OR4F29-201      OR4F29  939     CDS:1-939
      >ENST00000332831.4      ENSG00000284662.1       OR4F16-201      OR4F16  995     UTR5:1-19       CDS:20-958      UTR3:959-995

      2、改名操作

      导入 R 里:

      # 设置工作路径
      setwd('C:/Users/admin/Desktop/test/')
      # 读取id数据
      all_name <- read.delim('name.txt',header = F)
      # 添加列名
      colnames(all_name) <- c('tt_id','gene_id','tt_name','gene_name','tt_length','info1','info2','info3')
      # 查看行数
      nrow(all_name)
      [1] 106143

      可以看到最后 3 列相同的内容不在一列里,这是因为有些转录本只有 CDS 区或者 5’UTR 和 CDS 区,我们需要增加一列都是 CDS 区的位置信息:

      # 添加cds列
      library(stringr)
      # 新建储存cds信息的list
      cds_list <- list()
      # 循环
      for (i in 1:nrow(all_name)) {
          # 提取数据
          all_name_filter <- all_name[i,]
          # 分割info1
          info_1 <- unlist(strsplit(all_name_filter$info1,split = ':|-'))
          if (info_1[1] == 'UTR5') {
              # 如果第一个是UTR5,则取第二个元素
              cds_list[[i]] <- all_name_filter$info2
          }else if(info_1[1] == 'CDS') {
              # 如果第一个是CDS,则取第一个元素
              cds_list[[i]] <- all_name_filter$info1
          }else {
          }
      }

      # 查看个数
      length(cds_list)
      [1] 106143
      # 合并cds数据
      cds_res <- do.call('rbind',cds_list)
      # 储存在新数据里
      all_name_cds <- all_name
      all_name_cds$cds <- cds_res
      # 查看内容
      head(all_name_cds,4)
                     tt_id            gene_id    tt_name gene_name tt_length      info1        info2          info3
      1 >ENST00000641515.2  ENSG00000186092.7  OR4F5-201     OR4F5      2618  UTR5:1-60  CDS:61-1041 UTR3:1042-2618
      2 >ENST00000426406.4  ENSG00000284733.2 OR4F29-201    OR4F29       939  CDS:1-939
      3 >ENST00000332831.4  ENSG00000284662.1 OR4F16-201    OR4F16       995  UTR5:1-19   CDS:20-958   UTR3:959-995
      4 >ENST00000616016.5 ENSG00000187634.13 SAMD11-209    SAMD11      3465 UTR5:1-509 CDS:510-3044 UTR3:3045-3465
                 cds
      1  CDS:61-1041
      2    CDS:1-939
      3   CDS:20-958
      4 CDS:510-3044

      ###################################
      # 改名
      # 新建储存添加列名的list
      st_sp_list <- list()
      # 循环
      for (i in 1:nrow(all_name_cds)) {
          # 提取数据
          all_name_cds_filter <- all_name_cds[i,]
          # 分割info1
          info_cds <- unlist(strsplit(all_name_cds_filter$cds,split = ':|-'))
          # 去掉tt_id的 > 符号
          ttid <- unlist(strsplit(all_name_cds_filter$tt_id,split = '>'))[2]
          # 添加改名列
          all_name_cds_filter$new_name <- paste(all_name_cds_filter$gene_name,
                                                all_name_cds_filter$gene_id,
                                                ttid,
                                                all_name_cds_filter$tt_name,
                                                info_cds[2],
                                                info_cds[3],
                                                all_name_cds_filter$tt_length,
                                                sep = '_')
          # 储存结果
          st_sp_list[[i]] <- all_name_cds_filter
      }

      length(st_sp_list)
      [1] 106143
      # 合并数据
      st_sp_res <- do.call('rbind',st_sp_list)
      # 查看内容
      head(st_sp_res$new_name,4)
      [1] "OR4F5_ENSG00000186092.7_ENST00000641515.2_OR4F5-201_61_1041_2618"
      [2] "OR4F29_ENSG00000284733.2_ENST00000426406.4_OR4F29-201_1_939_939"
      [3] "OR4F16_ENSG00000284662.1_ENST00000332831.4_OR4F16-201_20_958_995"
      [4] "SAMD11_ENSG00000187634.13_ENST00000616016.5_SAMD11-209_510_3044_3465"

      这样我们就成功的改好了名字!接下来我们需要读入转录组序列,把它们重新命名。

      # 读取转录组序列fasta文件
      # 安装读取fasta文件的R包
      BiocManager::install('seqinr')
      library(seqinr)
      fasta <- read.fasta('gencode.v38.pc_transcripts.fa.gz',
                          seqtype = 'DNA',
                          forceDNAtolower = F)
      # fasta数据类型
      class(fasta)
      [1] "list"

      # 查看命名前的名称
      names(fasta)[1]
      [1] "ENST00000641515.2|ENSG00000186092.7|OTTHUMG00000001094.4|OTTHUMT00000003223.4|OR4F5-201|OR4F5|2618|UTR5:1-60|CDS:61-1041|UTR3:1042-2618|"

      # 命名
      names(fasta) <- st_sp_res$new_name

      # 查看命名后的名称
      names(fasta)[1]
      [1] "OR4F5_ENSG00000186092.7_ENST00000641515.2_OR4F5-201_61_1041_2618"

      # 保存命名后的数据
      write.fasta(sequences = fasta,
                  names = names(fasta),
                  file.out = 'name_changed.fa')

      3、提取代表性转录本

      主要思路,我们已经获得了转录本长度的信息,只要计算每个转录本的 CDS 区长度即可,然后按 CDS 列和转录本长度列降序排序取第一个就行了。

      library(tidyverse)
      # 提取代表性转录本

      # 新建储存筛选出代表性转录本的信息
      repsentive_tspt <- list()

      # gene
      gene <- unique(st_sp_res$gene_name)
      # 循环
      for (g in gene) {
          # 筛选gene
          ginfo <- st_sp_res %>% filter(gene_name == g)

          # 添加cds length
          # 切分cds列
          cds_list <- strsplit(ginfo$cds,split = ':|-')
          # 获取cds start和stop位置
          start <- sapply(cds_list, '[',2)
          stop <- sapply(cds_list, '[',3)
          # 计算赋值cds length
          ginfo$cds_length <- as.numeric(stop) - as.numeric(start)

          # 先按cds_length列降序,再按tt_length列降序
          ginfo_order <- ginfo %>% arrange(desc(cds_length,tt_length))

          # 筛选transcript
          tid <- ginfo$tt_id[1]
          list_name <- paste(g,tid,sep = '_')
          # 取第一行元素保存
          repsentive_tspt[[list_name]] <- ginfo_order[1,]

      }

      # 查看筛选结果数量
      length(repsentive_tspt)
      [1] 20336
      repsentive_results <- do.call('rbind',repsentive_tspt)
      # 查看有多少个基因
      length(unique(repsentive_results$gene_id))
      [1] 20336
      # 查看有多少个转录本
      length(unique(repsentive_results$tt_id))
      [1] 20336

      # 查看结果
      head(repsentive_results,4)
                                             tt_id            gene_id    tt_name gene_name tt_length      info1
      OR4F5_>ENST00000641515.2  >ENST00000641515.2  ENSG00000186092.7  OR4F5-201     OR4F5      2618  UTR5:1-60
      OR4F29_>ENST00000426406.4 >ENST00000426406.4  ENSG00000284733.2 OR4F29-201    OR4F29       939  CDS:1-939
      OR4F16_>ENST00000332831.4 >ENST00000332831.4  ENSG00000284662.1 OR4F16-201    OR4F16       995  UTR5:1-19
      SAMD11_>ENST00000616016.5 >ENST00000618323.5 ENSG00000187634.13 SAMD11-213    SAMD11      3468 UTR5:1-509
                                       info2          info3          cds
      OR4F5_>ENST00000641515.2   CDS:61-1041 UTR3:1042-2618  CDS:61-1041
      OR4F29_>ENST00000426406.4                                CDS:1-939
      OR4F16_>ENST00000332831.4   CDS:20-958   UTR3:959-995   CDS:20-958
      SAMD11_>ENST00000616016.5 CDS:510-3047 UTR3:3048-3468 CDS:510-3047
                                                                                            new_name cds_length
      OR4F5_>ENST00000641515.2      OR4F5_ENSG00000186092.7_ENST00000641515.2_OR4F5-201_61_1041_2618        980
      OR4F29_>ENST00000426406.4      OR4F29_ENSG00000284733.2_ENST00000426406.4_OR4F29-201_1_939_939        938
      OR4F16_>ENST00000332831.4     OR4F16_ENSG00000284662.1_ENST00000332831.4_OR4F16-201_20_958_995        938
      SAMD11_>ENST00000616016.5 SAMD11_ENSG00000187634.13_ENST00000618323.5_SAMD11-213_510_3047_3468       2537

      接下来该怎么从筛选的结果中取提取 fasta 数据呢?

      # 读取改名后的fasta
      named_fasta <- read.fasta('name_changed.fa',seqtype = 'DNA',forceDNAtolower = F)

      # 对代表性转录本循环
      tar_name <- repsentive_results$new_name

      # 建立储存筛选结果的list
      target_falist <- list()
      # 循环

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

      • 分享:
      作者头像
      weinfoadmin

      上一篇文章

      RNA-seq:Salmon 定量结果差异分析
      2021年9月10日

      下一篇文章

      画个CNS级别火山图!
      2021年9月10日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

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

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月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年
      在线支付 激活码

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