前言
廣告時(shí)間：如果覺(jué)得大家覺(jué)得推文有用，可以看下我的個(gè)人簡(jiǎn)介。給個(gè)關(guān)注，謝謝~
生信分析有時(shí)候需要提取最長(zhǎng)“轉(zhuǎn)錄本”，但是有時(shí)候往往會(huì)遇到一些極品，轉(zhuǎn)錄本和基因命名上看起來(lái)毫無(wú)關(guān)系。一些可視化的軟件提取最長(zhǎng)轉(zhuǎn)錄本，但我還未見(jiàn)能實(shí)現(xiàn)這個(gè)功能的。別人寫(xiě)的各種語(yǔ)言版本的提取最長(zhǎng)轉(zhuǎn)錄本還看不懂，也不敢用。自己的R語(yǔ)言能力弱的可憐，不能從頭到位只用tidyverse。沒(méi)辦法，只能用別人造的輪子了。
這篇推文，我將不同的部分使用一種或者多種方法實(shí)現(xiàn)。其中篩選轉(zhuǎn)錄本部分，R語(yǔ)言幾個(gè)簡(jiǎn)單的函數(shù)實(shí)現(xiàn)普適性的最長(zhǎng)轉(zhuǎn)錄組提取。
示例數(shù)據(jù)為Ensembl的pep序列。數(shù)據(jù)位置：http://ftp.ensemblgenomes.org/pub/plants/release-51/fasta/actinidia_chinensis/pep/ 這里我在使用的時(shí)候?qū)⑵渲孛麨?code>Red5.fa。
提取最長(zhǎng)轉(zhuǎn)錄本思路
思路分為以下三個(gè)步驟：（每個(gè)步驟選一個(gè)方法即可）
將fasta序列讀取為數(shù)據(jù)框格式（三種方法，建議新手使用第二種）；
篩選最長(zhǎng)轉(zhuǎn)錄本，期間要挑出基因ID，按基因ID進(jìn)行分組，按序列長(zhǎng)度排序，然后選擇最長(zhǎng)轉(zhuǎn)錄本（此時(shí)仍為data.frame格式）（一種方法）；
將序列讀出為fasta格式（兩種方法，建議新手使用第二種）；
難點(diǎn)：第二部分：這部分中選取分割的內(nèi)容需要根據(jù)自己的數(shù)據(jù)修改(這里pep兩側(cè)就各緊鄰了一個(gè)空格)" pep |gene:| transcript:"最為重要，需要注意前后空格。所以，對(duì)于新手，建議下載我用的示例數(shù)據(jù)先進(jìn)行復(fù)現(xiàn)，然后再用于自己的數(shù)據(jù)提取。
注意：本文linux命令與R命令混合，以$開(kāi)頭的為linux命令行，其余為R代碼。
perl腳本使用今天的次推《perl腳本練習(xí)：fasta格式與tab格式序列之間相互轉(zhuǎn)換(2)(適用范圍更廣)》中的腳本。R包biozhuoer見(jiàn)我之前的推文介紹：R包biozhuoer——將fasta序列讀入為tibble格式
分步提取最長(zhǎng)轉(zhuǎn)錄本
第一步：將fasta序列讀取為數(shù)據(jù)框
（1）第一種方法
使用一個(gè)R包biozhuoer實(shí)現(xiàn)。
library(biozhuoer)
df <- read_fasta("./fasta/Red5.fa")
class(df)
(2) 第二種方法
使用seqkit fx2tab 功能將fasta序列轉(zhuǎn)換為tab分割I(lǐng)D與seq的內(nèi)容，其中ID部分無(wú)‘>’，但是最后有一個(gè)單獨(dú)的制表符，R讀入后會(huì)有個(gè)空白行，這個(gè)問(wèn)題在讀入的同時(shí)會(huì)選擇相應(yīng)列。
seqkit有win系統(tǒng)版本，但是我還沒(méi)用過(guò)。
$ seqkit fx2tab -o Red5.fa.tsv Red5.fa #當(dāng)然，你也可以在R里使用system()執(zhí)行。
library(readr)
df <- read_delim("./fasta/Red5.fa.tsv", delim = "\t", col_names = FALSE,
col_select = c(X1,X2) #去除最后一個(gè)tab造成的空白
)
names(df) <- c("name", "seq")
class(df)
(3) 第三種方法
用我自己寫(xiě)的蹩腳perl 腳本 將fasta序列轉(zhuǎn)換為tab格式。
$ perl fa2tab_v2.pl -fa2tab Red5.fa Red5.fa.tsv
用R讀入為data.frame。
library(readr)
#用基礎(chǔ)函數(shù)read.table()無(wú)法將全部數(shù)據(jù)讀入。
df <- read_delim("./fasta/Red5.fa.tsv", delim = "\t", col_names = FALSE)
names(df) <- c("name", "seq")
class(df)
第二步：篩選最長(zhǎng)轉(zhuǎn)錄本（這里只寫(xiě)一種方法）
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其實(shí)可以合并為一個(gè)。主要是拆分寫(xiě)自己看著更能一目了然。
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")
第三步：保存篩選的最長(zhǎng)轉(zhuǎn)錄本
（1）第一種：使用biozhuoer包實(shí)現(xiàn)
library(biozhuoer)
#write_fasta()中的width參數(shù)是個(gè)假參數(shù)
write_fasta(longest_pep,"./fasta/Red5_longest_trans_pep.fa")
(2) 第二種：使用seqkit tab2fx功能實(shí)現(xiàn)
先保存為tab分割的序列，然后用seqkit tab2fx 功能轉(zhuǎn)化為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轉(zhuǎn)化為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
總結(jié)
看著眼花繚亂，仔細(xì)看看條理還是有的。選擇自己覺(jué)得能實(shí)現(xiàn)的方法提取轉(zhuǎn)錄本就可以了。
提取最長(zhǎng)轉(zhuǎn)錄本的第二步是關(guān)鍵，稍微學(xué)下這幾個(gè)函數(shù)，注意分割的時(shí)候是否帶空格就可以了。祝大家提取自己的最長(zhǎng)轉(zhuǎn)錄本成功。
