分步提取最长转录本 | 适用性更广 | 懂点代码就会

2021-12-07 本文已影响0人学生信的大叔

前言

广告时间：如果觉得大家觉得推文有用，可以看下我的个人简介。给个关注，谢谢~

生信分析有时候需要提取最长“转录本”，但是有时候往往会遇到一些极品，转录本和基因命名上看起来毫无关系。一些可视化的软件提取最长转录本，但我还未见能实现这个功能的。别人写的各种语言版本的提取最长转录本还看不懂，也不敢用。自己的R语言能力弱的可怜，不能从头到位只用tidyverse。没办法，只能用别人造的轮子了。

这篇推文，我将不同的部分使用一种或者多种方法实现。其中筛选转录本部分，R语言几个简单的函数实现普适性的最长转录组提取。

示例数据为Ensembl的pep序列。数据位置：http://ftp.ensemblgenomes.org/pub/plants/release-51/fasta/actinidia_chinensis/pep/ 这里我在使用的时候将其重命名为Red5.fa。

提取最长转录本思路

思路分为以下三个步骤：（每个步骤选一个方法即可）

将fasta序列读取为数据框格式（三种方法，建议新手使用第二种）；
筛选最长转录本，期间要挑出基因ID，按基因ID进行分组，按序列长度排序，然后选择最长转录本（此时仍为data.frame格式）（一种方法）；
将序列读出为fasta格式（两种方法，建议新手使用第二种）；

难点：第二部分：这部分中选取分割的内容需要根据自己的数据修改(这里pep两侧就各紧邻了一个空格)" pep |gene:| transcript:"最为重要，需要注意前后空格。所以，对于新手，建议下载我用的示例数据先进行复现，然后再用于自己的数据提取。

注意：本文linux命令与R命令混合，以$开头的为linux命令行，其余为R代码。
perl脚本使用今天的次推《perl脚本练习：fasta格式与tab格式序列之间相互转换(2)(适用范围更广)》中的脚本。R包biozhuoer见我之前的推文介绍：R包biozhuoer——将fasta序列读入为tibble格式

分步提取最长转录本

第一步：将fasta序列读取为数据框

（1）第一种方法
使用一个R包biozhuoer实现。

library(biozhuoer)
df <- read_fasta("./fasta/Red5.fa")
class(df)

(2) 第二种方法
使用seqkit fx2tab 功能将fasta序列转换为tab分割ID与seq的内容，其中ID部分无‘>’，但是最后有一个单独的制表符，R读入后会有个空白行，这个问题在读入的同时会选择相应列。
seqkit有win系统版本，但是我还没用过。

$ seqkit fx2tab -o Red5.fa.tsv Red5.fa #当然，你也可以在R里使用system()执行。


library(readr)
df <- read_delim("./fasta/Red5.fa.tsv", delim = "\t", col_names = FALSE, 
                 col_select = c(X1,X2) #去除最后一个tab造成的空白
                 )
names(df) <- c("name", "seq")
class(df)

(3) 第三种方法
用我自己写的蹩脚perl 脚本将fasta序列转换为tab格式。

$ perl fa2tab_v2.pl -fa2tab Red5.fa Red5.fa.tsv

用R读入为data.frame。

library(readr)
#用基础函数read.table()无法将全部数据读入。
df <- read_delim("./fasta/Red5.fa.tsv", delim = "\t", col_names = FALSE)
names(df) <- c("name", "seq")
class(df)

第二步：筛选最长转录本（这里只写一种方法）

library(tidyr)
library(dplyr)
library(stringr)
pep <- separate(df,name,c("pep_id","other1","gene_id","other2")," pep |gene:| transcript:") %>%
  mutate(length = str_length(pep$seq)) %>% 
  select(c("pep_id","gene_id","seq","length"))

#这里定义的pep和longest_pep其实可以合并为一个。主要是拆分写自己看着更能一目了然。
longest_pep <- pep %>% group_by(gene_id) %>% 
  arrange(desc(length),by_group = TRUE) %>% 
  slice_head(n=1) %>% 
  select(c("gene_id", "seq"))
names(longest_pep) <- c("name", "seq")

第三步：保存筛选的最长转录本

（1）第一种：使用biozhuoer包实现

library(biozhuoer)
#write_fasta()中的width参数是个假参数
write_fasta(longest_pep,"./fasta/Red5_longest_trans_pep.fa")

(2) 第二种：使用seqkit tab2fx功能实现
先保存为tab分割的序列，然后用seqkit tab2fx 功能转化为fasta格式序列。

library(readr)
write_delim(longest_pep, "./fasta/Red5.fa.tsv", delim = "\t",col_names = F)

$ seqkit tab2fx -o Red5_longest_trans_pep.fa Red5.fa.tsv

(3) 第三种：使用自己的蹩脚perl脚本
先保存为tab分割的序列，然后用自己的蹩脚脚本fa2tab_v2.pl转化为fasta格式序列。

library(readr)
write_delim(longest_pep, "./fasta/Red5.fa.tsv", delim = "\t",col_names = F)

$ perl fa2tab_v2.pl -tab2fa Red5.fa.tsv Red5_longest_trans_pep.fa

总结

看着眼花缭乱，仔细看看条理还是有的。选择自己觉得能实现的方法提取转录本就可以了。
提取最长转录本的第二步是关键，稍微学下这几个函数，注意分割的时候是否带空格就可以了。祝大家提取自己的最长转录本成功。