Как подсчитать перекрывающиеся димеры для нескольких последовательностей?

Question

Я должен посчитать количество перекрывающихся димеров (AA, AG, AC, AT, GA, GG, GC, GT, CC, CG, CA, CT, TT, TA, TG, TC) в нескольких последовательностях, используя Perl , Я написал следующий код, но он работает только для одной последовательности. Как я могу расширить его на несколько последовательностей?

#!/usr/bin/perl -w
open FH, "sample.txt";
$genome=<FH>;
%count=();
$count{substr($genome, $_, 2)}++ for 0..length($genome)-2;
print "AA: $count{AA}\n";
print "AG: $count{AG}\n";
print "AC: $count{AC}\n";
print "AT: $count{AT}\n";

print "TA: $count{TA}\n";
print "TG: $count{TG}\n";
print "TC: $count{TC}\n";
print "TT: $count{TT}\n";

print "GA: $count{GA}\n";
print "GG: $count{GG}\n";
print "GC: $count{GC}\n";
print "GT: $count{GT}\n";

print "CA: $count{CA}\n";
print "CG: $count{CG}\n";
print "CC: $count{CC}\n";
print "CT: $count{CT}\n";

Мне нужно:

1. рассчитывает для каждой последовательности и
2. всего насчитывается

пример ввода: sample.txt

ATGGGCTCCTCCGCCATCACCGTGAGCTTCCTCCTCTTTCTGGCATTTCAGCTCCCAGGGCAAACAGGAGCAAATCCCGTGTATGGCTCTGTGTCCAATGCAGACCTGATGGATTTCAAGTAAAAG
ATGGTGAGAAAATGGGCCCTGCTCCTGCCCATGCTGCTCTGCGGCCTGACTGGTCCCGCACACCTCTTCCAGCCAAGCCTGGTGCTGGAGATGGCCCAGGTCCTCTTGGATAACTACTGCTTCCCAGAGAACCTGATGGGGATGCAGGGAGCCATCGAGCAGGCCATCAAAAGTCAGGAGATTCTGTCTATCTCAGACCCTCAGACTCTGGCCCATGTGTTGACAGCTGGGGTGCAGAGCTCCTTGAATGACCCTCGCCTGGTCATCTCCTATGAGCCCAGCACCCTCGAGGCCCCTCCGCGAGCTCCAGCAGTCACGAACCTCACACTAGAGGAAATCATCGCAGGGCTGCAGGATGGCCTTCGCCATGAGATTCTGGAAGGCAATGTGGGCTACCTGCGGGTGGACGACATCCCGGGCCAGGAGGTGATGAGCAAGCTGAGGAGCTTCCTGGTGGCCAACGTCTGGAGGAAGCTCGTGAACACGTCCGCCTTGGTGCTGGACCTCCGGCACTGCACTGGGGGACACGTGTCTGGCATCCCCTATGTCATCTCCTACCTGCACCCAGGGAGCACAGTCTCGCACGTGGACACCGTCTACGACCGCCCCTCCAACACAACCACTGAGATCTGGACCCTGCCTGAAGCCCTGGGAGAGAAGTACAGTGCAGACAAGGATGTGGTGGTCCTCACCAGCAGCCGCACGGGGGGCGTGGCTGAGGACATCGCTTACATCCTCAAACAGATGCGCAGGGCCATCGTGGTGGGCGAGCGGACTGTTGGGGGGGCTCTGAACCTCCAGAAGCTGAGGGTAGGCCAGTCCGACTTCTTTCTCACTGTGCCTGTGTCCAGATCCCTGGGGCCCCTGGGTGAGGGCAGCCAGACGTGGGAGGGCAGTGGGGTGCTGCCCTGTGTGGGGACACCGGCCGAGCAGGCCCTGGAGAAAGCCCTGGCCGTTCTCATGCTGCGCAGGGCCCTGCCAGGAGTCATTCAGCGCCTTCAGGAGGCGCTGCGCGAGTACTACACGCTGGTGGACCGTGTGCCCGCCCTGCTGAGCCACCTGGCCGCCATGGACCTGTCCTCGGTGGTCTCCGAGGACGATCTGGTCACTAAGCTCAATGCTGGCCTGCAGGCTGTGTCTGAGGACCCCAGGCTCCAGGTGCAGGTGGTCAGACCCAAAGAAGCCTCTTCTGGGCCTGAGGAAGAAGCTGAAGAACCTCCAGAGGCGGTCCCGGAAGTGCCCGAGGACGAGGCTGTTCGGCGGGCTCTGGTGGACTCCGTGTTCCAGGTTTCTGTGCTGCCGGGCAACGTGGGCTACCTGCGCTTCGACAGTTTCGCTGATGCCTCTGTCCTGGAGGTGCTGGGCCCCTACATCCTGCACCAGGTGTGGGAGCCCCTGCAGGACACGGAGCACCTCATCATGGACCTGCGGCAGAACCCCGGGGGGCCGTCCTCCGCGGTGCCCCTGCTGCTCTCCTACTTCCAGAGCCCTGACGCCAGCCCCGTGCGCCTCTTCTCCACCTACGACCGGCGCACCAACATCACACGCGAGCACTTCAGCCAGACGGAGCTGCTGGGCCGGCCCTACGGCACCCAGCGTGGCGTGTACCTGCTCACTAGCCACCGCACCGCCACCGCGGCCGAGGAGCTGGCCTTCCTCATGCAGTCACTGGGCTGGGCCACGCTGGTGGGCGAGATCACCGCGGGCAGCCTGCTGCACACACACACAGTATCCCTGCTGGAGACGCCCGAGGGCGGCCTGGCGCTCACGGTGCCTGTGCTCACCTTCATCGACAACCATGGCGAGTGCTGGCTGGGGGGCGGTGTGGTCCCCGATGCCATTGTGCTGGCCGAGGAAGCCCTAGACAGAGCCCAGGAGGTGCTGGAGTTCCACCGAAGCTTGGGGGAGTTGGTGGAAGGCACGGGGCGCCTGCTGGAGGCCCACTACGCTCGGCCAGAGGTCGTGGGGCAGATGGGTGCCCTGCTGCGAGCCAAGCTGGCCCAGGGGGCCTACCGCACCGCGGTGGACCTGGAGTCGCTGGCTTCCCAGCTTACGGCCGACCTGCAGGAGATGTCTGGGGACCACCGTCTGCTGGTGTTCCACAGCCCCGGCGAAATGGTGGCTGAGGAGGCGCCCCCACCGCCTCCCGTCGTCCCCTCCCCGGAGGAGCTGTCCTATCTCATCGAGGCCCTGTTCAAGACTGAGGTGCTGCCCGGCCAGCTGGGCTACCTGCGTTTCGACGCCATGGCTGAGCTGGAGACGGTGAAGGCCGTCGGGCCACAGCTGGTGCAGCTGGTGTGGCAGAAGCTGGTGGACACGGCCGCGCTGGTGGTCGACCTGCGCTACAACCCCGGCAGCTACTCCACAGCCGTGCCTCTACTCTGCTCCTACTTCTTCGAGGCAGAGCCCCGCCGGCACCTCTACTCTGTCTTTGACAGGGCCACGTCAAGGGTCACAGAGGTATGGACCCTGCCCCACGTTACAGGCCAGCGCTATGGCTCCCACAAGGACCTCTACGTTCTGGTGAGCCACACCAGCGGTTCAGCAGCTGAGGCTTTTGCTCACACCATGCAGGATCTGCAGCGAGCCACCATCATCGGGGAGCCCACGGCCGGAGGGGCACTCTCCGTGGGAATCTACCAGGTGGGCAGCAGCGCCTTATACGCCTCCATGCCCACGCAGATGGCCATGAGTGCCAGCACCGGCGAGGCCTGGGATCTGGCTGGGGTGGAGCCGGACATCACTGTGCCCATGAGCGTGGCCCTCTCCACAGCCCGGGACATAGTGACCCTGCGTGCCAAGGTGCCCACTGTGCTGCAGACAGCTGGGAAGCTCGTAGCGGATAACTACGCCTCCCCTGAGCTGGGAGTCAAGATGGCAGCCGAACTGAGTGGTCTGCAGAGCCGCTATGCCAGGGTGACCTCAGAAGCCGCCCTCGCCGAGCTGCTGCAAGCCGACCTGCAGGTGCTGTCCGGGGACCCACACCTGAAGACAGCTCATATACCTGAGGATGCCAAAGACCGCATTCCTGGCATTGTACCCATGCAGTAACAG
ATGGACATGATGGACGGCTGCCAGTTCTCGCCCTCTGAGTACTTCTACGACGGCTCCTGCATCCCATCCCCCGACGGTGAGTTCGGGGACGAGTTTGAGCCGCGAGTGGCTGCTTTCGGGGCTCACAAGGCAGACCTGCAAGGCTCAGACGAGGACGAGCACGTGCGAGCACCCACGGGCCACCACCAGGCCGGCCACTGCCTCATGTGGGCCTGCAAAGCATGCAAAAGGAAGTCCACCACCATGGATCGGCGGAAGGCGGCCACCATGCGCGAGCGGAGACGCCTGAAGAAGGTCAACCAGGCTTTCGACACGCTCAAGCGGTGCACCACGACCAACCCTAACCAGAGGCTGCCCAAGGTGGAGATCCTCAGGAATGCCATCCGCTACATTGAGAGTCTGCAGGAGCTGCTTAGGGAACAGGTGGAAAACTACTATAGCCTGCCGGGGCAGAGCTGCTCTGAGCCCACCAGCCCCACCTCAAGTTGCTCTGATGGCATGTAAATG

Answer 1

#!/usr/bin/perl

use strict;

my @dimers = qw(AA AG AC AT GA GG GC GT CC CG CA CT TT TA TG TC);
my @dimers_totals;

my $i = 1;
for(<>)
{
    my $sequence = $_;
    print("Sequence $i:\n");
    my $j = 0;
    for(@dimers)
    {
        my $number =()= $sequence =~ /$_/gi;
        $dimers_totals[$j++] += $number;
        print "\t$_: $number\n"
    }
    print("\n");
    $i++;
}

print("Totals:\n");
$i = 0;
for(@dimers)
{
    print("\t$_: $dimers_totals[$i++]\n");
}

Используйте как

./dimers_count.pl < sample.txt

Answer 2

Моей первой мыслью было глобальное сопоставление в скалярном контексте с настройкой pos () для резервного копирования.Таким образом, мне нужно сканировать строку только один раз для всех димеров (тогда как другие ответы пока сканируют последовательность один раз для каждого димера).Я тоже поддерживаю хэши;один для всего и один для последовательности.Я использую специальную переменную $., которая содержит текущий номер строки ввода для обозначения последовательности:

use Data::Printer;

while( <DATA> ) {
    while( /\G([ATGC]{2})/g ) {
        $total{$1}++;
        $by_sequence{$.}{$1}++;
        pos() = pos() - 1
        }
    }

p( %total );
p( %by_sequence );

То же самое с substr потребовало немного больше работы, потому чтоЯ не мог ограничить подстроку только символами, которые имели значение.Я должен был обратить внимание на новые строки и нечетные числа в конце:

use Data::Printer;

LINE: while( <DATA> ) {
    chomp;
    my $pos = 0;
    SUB: while( my $sub = substr( $_, $pos, 2 ) ) {
        last SUB unless 2 == length $sub;
        $total{$sub}++;
        $by_sequence{$.}{$sub}++;
        $pos++;
        }
    }

p( %total );
p( %by_sequence );

Вывод из Data :: Printer очень хорош.Я использую его в предпочтении Data :: Dumper , так как я узнал об этом на YAPC :: Brasil.(И, кстати, его автор будет в YAPC :: NA в Мэдисоне в этом году , чтобы поговорить об этом ):

{
    AA   101,
    AC   215,
    AG   268,
    AT   106,
    CA   286,
    CC   388,
    CG   201,
    CT   310,
    GA   239,
    GC   376,
    GG   369,
    GT   168,
    TA   61,
    TC   206,
    TG   317,
    TT   73
}
{
    1   {
        AA   9,
        AC   3,
        AG   8,
        AT   8,
        CA   11,
        CC   12,
        CG   3,
        CT   11,
        GA   5,
        GC   9,
        GG   8,
        GT   6,
        TA   2,
        TC   13,
        TG   10,
        TT   7
    },
    2   {
        AA   68,
        AC   177,
        AG   219,
        AT   80,
        CA   230,
        CC   329,
        CG   168,
        CT   264,
        GA   195,
        GC   316,
        GG   317,
        GT   146,
        TA   50,
        TC   169,
        TG   271,
        TT   54
    },
    3   {
        AA   24,
        AC   35,
        AG   41,
        AT   18,
        CA   45,
        CC   47,
        CG   30,
        CT   35,
        GA   39,
        GC   51,
        GG   44,
        GT   16,
        TA   9,
        TC   24,
        TG   36,
        TT   12
    }
}

Answer 3

Ключом к этому является то, что вы имели в качестве функциональной единицы того, что вы пытались построить.Ваш код работал для одной последовательности;расширение его на несколько последовательностей - это просто вопрос избавления от предположения, что рассматривается только одна последовательность:

use strict;
use warnings;
use 5.010;

use Data::Dumper ();

sub count_dimers {
    my ($sequence) = @_;

    my %counts;
    $counts{substr($sequence, $_, 2)}++ for 0..length($sequence) - 2;

    my @dimers = qw(AA AG AC AT GA GG GC GT CC CG CA CT TT TA TG TC);
    %counts = map { $_ => $counts{$_} } @dimers;

    return %counts;
}

open(my $fh, '<', 'sample.txt');

my @counts_by_sequence;
while (my $sequence = <$fh>) {
    my %sequence_counts = count_dimers($sequence);
    push @counts_by_sequence, \%sequence_counts;
}

my %total_counts;
for my $sequence_counts (@counts_by_sequence) {
    for my $dimer (keys %$sequence_counts) {
        $total_counts{$dimer} += ${ $sequence_counts}{$dimer};
    }
}

say Data::Dumper->Dump(
    [\%total_counts, \@counts_by_sequence],
    [qw(total_count counts_by_sequence)]
);

Я оставил вывод в качестве упражнения для читателя, но общую форму преобразованиядолжно быть очевидно: то, что раньше было целой программой, теперь является функцией, которая вызывается один раз для каждой последовательности с сохранением результатов на каждый вызов и итоговых результатов по всем вызовам.