其他
提取 ensembl 的 gtf 文件中最长转录本信息
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%
微信扫一扫付费阅读本文