当 DNA 序列中存在某种模式时，检索编码氨基酸

Question

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

输入文件：

>sequence1
ATGGCGCATAGTAATGC
>sequence2
ATGATAGTAATGCGCGC

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

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

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

---

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

请忘记原帖和sequence1和sequence2。

大家好，很抱歉造成混乱。输入的 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。我不希望程序检索四个或更多重复的编码氨基酸。

再次感谢，伯纳多

Answer 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

Answer 2

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

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

input_seq.txt：

GC
ATA
ATAGTA
ATGATA

dna.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")}

test.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"}

文件必须全部在同一个目录中。

祝你好运！