写一个 R 版本的 ORFfinder

引言


预测 开放阅读框 的软件和一下 非编码 RNA 翻译微肽 的数据库都有一些。根据提供的序列去寻找可能的 ORF(Open Reading frame) 。比如 NCBI 的 ORFfinder,我们需要满足最基本符合翻译的原则。
地址:https://www.ncbi.nlm.nih.gov/orffinder/
只需要粘贴 FASTA 序列即可:

页面端只能处理一定长度的序列,而且不能多序列处理,页面提供了一个本地安装版:


github 上淘宝了一下,找到一个 python 版本的,但是也只能单次处理一条序列。
地址:https://github.com/Chokyotager/ORFFinder

安装

自己造个:

install.packages(“devtools”)

devtools::install_github(“junjunlab/ORFfinderR”)

or

remotes::install_github(“junjunlab/ORFfinderR”)
python 转 R


把 python 版的 ORFFinder 代码改一下,增加一个批量处理并输出的结果。
放到 R 里面调用并进行分析返回结果。

batch get orf

def batch_get_orf(sequence=””,output_file=””, **kwargs):
from Bio import SeqIO

orf_lst = []

for sequence in SeqIO.parse(sequence, "fasta"):
    orf = getORFProteins(sequence, **kwargs)
    orf_lst.append(orf)

flatern = []
for i in orf_lst:
    for j in i:
        flatern.append(j)

out_file = open(output_file, 'w')
for i in flatern:
    out_file.write("\t".join([str(i["id"]),str(i["t_len"]),str(i["start"]),str(i["end"]),str(i["frame"]),str(i["sense"]),str(i["length"]),str(i["protein"])]) + '\n')
out_file.close()

测试的两条序列:

非常简单的使用,同时返回蛋白序列:
library(ORFfinderR)

df <- get_orf(sequence = “../test.fa”,remove_nested = “TRUE”,output_file = “pkg_test.txt”)

和 API 调用是一样的结果:

可视化一下:
df <- get_orf(sequence = “../test.fa”,output_file = “pkg_test.txt”)

orf_plot(data = df)

寻找 ORF


这里我们还是自己手动去写一个函数,基本原理是分别从给定序列的正链和负链寻找用户给定的起始密码子(ATG)和终止密码子(TAA,TAG,TGA)所有位置,然后筛选不是 3 的倍数的序列即可。
一句代码搞定:
tt <- find_orf(seq_obj = “../test.fa”)

return amino acid

tt <- find_orf(seq_obj = “../test.fa”,return_AA = T)

返回所有的可能情况:
tt <- find_orf(seq_obj = “../test.fa”,return_AA = T,filter_orf = F)

可视化一下:
orf_plot(data = tt)


可以看到,有很多 具有相同的 start_codon 或者 stop_codon 的序列。所以需要有一个筛选的规则的。filter_orf 的原则是:相同的 start_codon 选最短的序列,相同的 stop_codon 选最长的序列。 具体翻译的情况还得需要实验去验证的。
结尾


路漫漫其修远兮,吾将上下而求索。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。

声明:文中观点不代表本站立场。本文传送门:https://eyangzhen.com/423684.html

联系我们
联系我们
分享本页
返回顶部