从多序列比对中提取和连接外显子

Question

使用我(公平)的基本编码技巧，我已经将一个脚本组合在一起，它将解析一个对齐的多fasta文件(一个多序列对齐)，并提取两个指定列之间的所有数据。

use Bio::SimpleAlign;  
use Bio::AlignIO;
$str = Bio::AlignIO->new(-file => $inputfilename, -format => 'fasta');  
$aln = $str->next_aln();
$mini = $aln->slice($array[0], $array[1]);
$out = Bio::AlignIO->new(-file => $array[3], -format => 'fasta');
$out->write_aln($mini);

我所拥有的问题是，我希望能够从相同的对齐中分割多个区域，然后在写入一个外部文件之前加入这些区域。复杂的是我想要提供一个文件的协调列表，其中每一行包含两个或更多的协调，其中数据应该被提取和连接。

下面是一个协调文件的示例

ORF1, 10, 50, exon1 # The above line should produce a slice between columns 10-50 and write to an outfile  
ORF2, 70, 140, exon1  
ORF2, 190, 270, exon2  
ORF2, 500, 800, exon3 # Data should be extracted between the ranges specified here and in the above two lines and then joined (side by side) to produce the outfile.  
ORF3, 1200, 1210, exon1  
etc etc 

下面是一个对齐fasta文件的(小)示例

\>Sample1
ATGGCGACCGTGCACTACTCCCGCCGACCTGGGACCCCGCCGGTCACCCTCACGTCGTCC
CCCAGCATGGATGACGTTGCGACCCCCATCCCCTACCTACCCACATACGCCGAGGCCGTG
GCAGACGCGCCCCCCCCTTACAGAAGCCGCGAGAGTCTGGTGTTCTCCCCGCCTCTTTTT
CCTCACGTGGAGAATGGCACCACCCAACAGTCTTACGATTGCCTAGACTGCGCTTATGAT
GGAATCCACAGACTTCAGCTGGCTTTTCTAAGAATTCGCAAATGCTGTGTACCGGCTTTT
TTAATTCTTTTTGGTATTCTCACCCTTACTGCTGTCGTGGTCGCCATTGTTGCCGTTTTT
CCCGAGGAACCTCCCAACTCAACTACATGA  
\>Sample2
ATGGCGACCGTGCACTACTCCCGCCGACCTGGGACCCCGCCGGTCACCCTCACGTCGTCC
CCCAGCATGGATGACGTTGCGACCCCCATCCCCTACCTACCCACATACGCCGAGGCCGTG
GCAGACGCGCCCCCCCCTTACAGAAGCCGCGAGAGTCTGGTGTTCTCCCCGCCTCTTTTT
CCTCACGTGGAGAATGGCACCACCCAACAGTCTTACGATTGCCTAGACTGCGCTTATGAT
GGAATCCACAGACTTCAGCTGGCTTTTCTAAGAATTCGCAAATGCTGTGTACCGGCTTTT
TTAATTCTTTTTGGTATTCTCACCCTTACTGCTGTCGTGGTCGCCATTGTTGCCGTTTTT
CCCGAGGAACCTCCCAACTCAACTACATGA  

我认为应该有一个很简单的方法来解决这个问题，有可能使用第一列中的信息，加上外显子号，但我无法用我的一生来想出这是如何做到的。

有人能帮我吗？

Answer 1

1. 您发布的对齐fasta文件--至少在堆栈溢出网页上显示的文件--没有编译。根据format，描述行应该以>开始，而不是以\>开始。
2. 确保使用use strict; use warnings;运行所有Perl程序。这将方便调试。
3. 您没有填充@array。因此，您可能会得到这样的错误： 在perl-5.24.0/lib/site_perl/5.24.0/Bio/SimpleAlign.pm第1086行的模式匹配(m//)中使用未初始化的值$start，在第16行使用未初始化的值$start (.)或者perl-5.24.0/lib/site_perl/5.24.0/Bio/SimpleAlign.pm第1086行的字符串，第16行。不[]堆栈:错误：：抛出堆栈:Bio：：Root：：perl-5.24.0/lib/site_perl/5.24.0/Bio/Root/Root.pm:444堆栈:Bio：：SimpleAlign：：片perl-5.24.0/lib/site_perl/5.24.0/Bio/SimpleAlign.pm:1086堆栈: fasta.pl:26

一旦你分配了合理的价值，例如，

@array = (1,17);

..。你会得到更合理的结果：

$ perl fasta.pl
>Sample1/1-17
ATGGCGACCGTGCACTA
>Sample2/1-17
ATGGCGACCGTGCACTA

哈哈！