如何切取特定数目的fasta序列并生成ORFs

Question

我有如下所示的文件：

 CDS             join(36..56,37..67)
 CDS             36..183
 CDS             457..565
 CDS             join(505..519,521..596)
 CDS             join(577..591,725..770)
 CDS             join(516..591,725..899)
 CDS             508..556
 CDS             571..841
 CDS             complement(619..788)
 CDS             843..863

我想打印特定数量的核苷酸范围，如文件(序列是从另一个文件“sequence.fasta”读取)。例如，对于sequence.fasta文件，如：

>gi1234 HIVgenome|NC_909999.1
AACTGCGTGTGTGTCCACACAACACTGGGGGACACACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACAGGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTGTACAACACTGGGGGACAGTGACGACGACAACACTGGGGGACACGAGCGTTGTGAGCAGGTGACAACACTGGGGGACAGTGTTTTTACAACACTGGGGGACATTTTTGAGCAGCGACGCAGCGTTGTGGGGTGTGTCGGAAGGTGTGTCGTGTGTCGTGTGTC

输出应该是

36  -  56   ACAACAACAACACTGGGGGAC 

37  -  67   CAACAACAACACTGGGGGACAACACTGGGAC

等等.

直到

843 - 863   GTGT....

通过shell脚本最简单的方法是什么？

Answer 1

这个问题需要比这个论坛可能提供的更大的编程努力(我是以这种编程为生)。

DDBJ/ENA/GenBank文件格式 (问题中的第一个文件)是复杂的，它允许CDS(基因组序列的编码部分)不仅是简单的或连接的，而且可以补充和组合它们。此外，位置坐标可能有修饰符，对于一般的解决方案，它需要处理。

最好是询问当地的生物信息专家(或程序员)或生物信息学论坛，如StackExchange 生物信息学站点。他们将为您指出现有的工具来完成这类事情，或者，了解生物信息人员，给您一些古怪的BioPerl/BioPython脚本，这些脚本可能会更经常地工作;-)

一种可能的方法是使用GenBank特征抽取器，但是在线使用它很可能是除了小数据集以外的最佳选择。