1

当 DNA 序列中存在某种模式时，我想检索编码氨基酸。例如，模式可以是：ATAGTA。所以，当有：

输入文件：

>sequence1
ATGGCGCATAGTAATGC
>sequence2
ATGATAGTAATGCGCGC

理想的输出将是一个表格，其中每个氨基酸的次数由模式编码。在序列 1 中，模式只编码一个氨基酸，但在序列 2 中，它编码两个。我想让这个工具可以扩展到数千个序列。我一直在考虑如何完成这项工作，但我只想：替换所有与模式不同的核苷酸，翻译剩下的内容并获得编码氨基酸的摘要。

请让我知道是否可以通过现有工具执行此任务。

谢谢你的帮助。一切顺利，贝尔纳多

编辑（由于我的帖子产生的混乱）：

请忘记原始帖子以及序列1和序列2。

大家好，很抱歉造成混乱。输入的 fasta 文件是使用“FeatureExtract”工具（ http://www.cbs.dtu.dk/services/FeatureExtract/download.php ）从 GenBank 文件派生的 *.ffn 文件，因此可以想象它们已经在帧（+1），并且不需要在与+1不同的帧中编码氨基酸。

我想知道以下序列编码的是哪种氨基酸：

AGAGAG
GAGAGA
CTCTCT
TCTCTC

我想获得编码氨基酸的唯一字符串是三个 AG、GA、CT 或 TC 的重复，分别是 (AG)3、(GA)3、(CT)3 和 (TC)3。我不希望程序检索四个或更多重复的编码氨基酸。

再次感谢，贝尔纳多

4

2 回答 2

1

这是一些至少可以帮助您入门的代码。例如，您可以像这样运行：

./retrieve_coding_aa.pl file.fa ATAGTA

内容retrieve_coding_aa.pl：

#!/usr/bin/perl 

use strict;
use warnings;

use File::Basename;
use Bio::SeqIO;
use Bio::Tools::CodonTable;
use Data::Dumper;

my $pattern = $ARGV[1];

my $fasta = Bio::SeqIO->new ( -file => $ARGV[0], -format => 'fasta');

while (my $seq = $fasta->next_seq ) {

    my $pos = 0;

    my %counts;

    for (split /($pattern)/ => $seq->seq) {

        if ($_ eq $pattern) {

            my $dist = $pos % 3;

            unless ($dist == 0) {

                my $num = 3 - $dist;

                s/.{$num}//;

                chop until length () % 3 == 0;
            }

            my $table = Bio::Tools::CodonTable->new();

            $counts{$_}++ for split (//, $table->translate($_));
        }

        $pos += length;
    }

    print $seq->display_id() . ":\n";

    map {

        print "$_ => $counts{$_}\n"
    }
    sort {

        $counts{$a} <=> $counts{$b}
    }
    keys %counts;

    print "\n";
}

以下是使用示例输入的结果：

sequence1:
S => 1

sequence2:
V => 1
I => 1

该类Bio::Tools::CodonTable还支持非标准密码子使用表。id您可以使用指针更改表格。例如：

$table = Bio::Tools::CodonTable->new( -id => 5 );

或者：

$table->id(5);

有关更多信息，包括如何检查这些表格，请参阅此处的文档：http: //metacpan.org/pod/Bio ::Tools::CodonTable

于 2013-11-13T13:24:32.803 回答
0

我会坚持你想要的第一个版本，因为附录只会让我更加困惑。（帧？）我只在 sequence2 中找到了一次 ATAGTA，但我假设您也想要镜像/反向序列，在这种情况下是 ATAGATA。好吧，我的脚本没有这样做，所以你必须在 input_sequences 文件中写两次，但我认为这应该没问题。

我使用像你这样的文件，我称之为“dna.txt”和一个名为“input_seq.txt”的输入序列文件。结果文件是 dna.txt 文件中模式及其出现的列表（包括重叠结果，但它可以设置为非重叠，如 awk 中所述）。

input_seq.txt：

GC
ATA
ATAGTA
ATGATA

脱氧核糖核酸.txt：

>sequence1
ATGGCGCATAGTAATGC
>sequence2
ATGATAGTAATGCGCGC

结果.txt：

GC,6
ATA,2
ATAGTA,2
ATGATA,1

代码是 awk 调用另一个 awk（但其中一个很简单）。您必须运行“./match_patterns.awk input_seq.txt”才能生成结果文件。：

*match_patterns.awk:*

#! /bin/awk -f
{return_value= system("awk -vsubval="$1" -f test.awk dna.txt")}

测试.awk：

#! /bin/awk -f
{string=$0
do
{
where = match(string, subval)
# code is for overlapping matches (i.e ATA matches twice in ATATAC)
# for non-overlapping replace +1 by +RLENGTH in following line
if (RSTART!=0){count++; string=substr(string,RSTART+1)}
}
while (RSTART != 0)
}
END{print subval","count >> "results.txt"}

文件必须全部位于同一目录中。

祝你好运！

于 2013-11-13T14:48:45.570 回答