用R提取代表转录本
感谢老俊俊的大力支持。我们会每日跟新，欢迎您关注老俊俊的生信笔记。
用 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”，“生信星球”的支持！
你可能也喜欢
白介素-17受体信号的自主激活，维持炎症并促进疾病进展
4 8月, 2023
MCT4依赖的乳酸分泌抑制LKB1缺陷肺腺癌的抗肿瘤免疫
18 7月, 2023
lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包，干货满满
30 1月, 2023
