查看原文
其他

提取 ensembl 的 gtf 文件中最长转录本信息

JunJunLab 老俊俊的生信笔记 2022-08-15

1、引言

一个基因有多个转录本,有时候我们只想保留一个转录本,我们可以提取出转录本序列最长的那个最为代表,当然也是很多文章里经常干的,写个循环把每个基因的转录本序列加和,然后取最大值,最后保留对应信息即可。

2、下载 ensembl 的注释文件

可以用迅雷或者去 ensembl 官网下载:

# 设置工作目录
setwd('c:/Users/admin/Desktop/')

# 下载gtf文件
download.file(url = 'http://ftp.ensembl.org/pub/release-104/gtf/homo_sapiens/Homo_sapiens.GRCh38.104.gtf.gz',
              destfile = 'Homo_sapiens.GRCh38.104.gtf.gz',
              method = 'libcurl',
              mode = 'wb')

3、导入 GTF 文件

# 加载R包
library(rtracklayer)
library(tidyverse)

# 导入 gtf 文件
gtf <- import('./Homo_sapiens.GRCh38.104.gtf.gz',format = 'gtf') %>%
  as.data.frame()

# 查看基因数量
length(na.omit(unique(gtf$gene_name)))
# [1] 39397

# 查看转录本数量
length(na.omit(unique(gtf$transcript_id)))
# [1] 236920

236920/39398
# [1] 6.013503

可以看到平均每个基因至少有 6 个转录本。

4、泛式函数版

拿 200 个基因测试:

微信扫一扫付费阅读本文

可试读22%

微信扫一扫付费阅读本文

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存