perl - 複数の配列で重複する二量体をカウントするにはどうすればよいですか?

Question

Perl を使用して、複数のシーケンスで重複するダイマー (AA、AG、AC、AT、GA、GG、GC、GT、CC、CG、CA、CT、TT、TA、TG、TC) の数をカウントする必要があります。次のコードを書きましたが、1 つのシーケンスでしか機能しません。複数のシーケンスに拡張するにはどうすればよいですか?

#!/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

私が最初に考えたのは、 1 つをバックアップするためにpos()を調整した、スカラー コンテキストでのグローバル マッチでした。このように、すべてのダイマーに対して文字列を 1 回スキャンするだけで済みます (これまでの他の回答では、すべてのダイマーに対してシーケンスを 1 回スキャンしています)。私もハッシュを維持しています。1 つは合計用、もう 1 つは順序用です。現在の入力行番号を保持する特別な変数を使用して$.、シーケンスにラベルを付けています。

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からの出力は非常に優れています。YAPC::Brasil で見つけて以来、Data::Dumperよりも優先して使用しています。(ちなみに、その作者は今年、マディソンの 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

これの鍵は、構築しようとしていたものの機能単位としてこれまでに得たものを利用することです。あなたのコードは 1 つのシーケンスで機能しました。それを複数のシーケンスに拡張することは、1 つのシーケンスのみが考慮されているという仮定を取り除くことです。

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)]
);

読者への演習として出力を残しておきましたが、変換の一般的な形は明らかなはずです: かつてプログラム全体だったものは、今ではシーケンスごとに 1 回呼び出され、呼び出しごとの結果が保存され、結果が渡される関数になっています。すべての通話の合計。