perlスクリプトの不具合

Question

いつもお世話になっております。自身で書いたperlスクリプトが思うように動かずに困っています。どなたか御教授頂けたら幸いです。

入力ファイルについて
・ file.fa

1-5853

AAACATCCTTGACTGAAAGCT

2-3440

ATGACCTATGATTTGACAGCT

3-2117

TGAGAACTGAATTCCATAGGCTG

4-2013

TGTAAACATCCTTGACTGAAAGCT

5-1966

AACCCGTAGATCCGAACTTGTG
・
・
・
数十万行続く
">"以下の"数字-数字"がIDに相当し、1行下の文字列（15〜30字）と対応

・ blast_result.tsv
1-5853 ssa-miR-30d-5p 100.000 21 0 0 1 21 4 24 8.09e-07 42.1
3-2117 chi-miR-146b-5p 100.000 23 0 0 1 23 1 23 6.22e-08 46.1
3-2117 ggo-miR-146b 100.000 23 0 0 1 23 1 23 6.22e-08 46.1
4-2013 ssa-miR-30d-5p 100.000 24 0 0 1 24 1 24 1.70e-08 48.1
4-2013 oha-miR-30e-5p 95.833 24 1 0 1 24 1 24 4.15e-06 40.1
4-2013 chi-miR-30e-5p 95.833 24 1 0 1 24 1 24 4.15e-06 40.1
4-2013 ssa-miR-30c-5p 95.833 24 1 0 1 24 1 24 4.15e-06 40.1
4-2013 sha-miR-30e 95.833 24 1 0 1 24 1 24 4.15e-06 40.1
5-1966 aca-miR-100 100.000 22 0 0 1 22 1 22 2.25e-07 44.1
・
・
・
数十万行続く
1列目がfile.fa中のいずれかのIDとなっている。

実現したいこと
blast_result.tsvの1列目を、IDに対応する文字列に置換して出力したい。そのために以下のようなスクリプトを書いてみた。

(@ARGV == 3) or die "usage: $0 file.fa blast_result.tsv DNAorRNA\n";

($ARGV[2] == "DNA || dna || RNA || rna") or die "please select DNA(dna) or RNA(rna)\n";

$file = $ARGV[0];

open($fh, $file) or die "cannot open the fasta!!\n";
while (($name, $seq) = fasta_get($fh)) {
        @con = $seq;
        $name =~ /^(\S+)/;
        $name = $1;
      for (0..$#con) {
         if ($ARGV[2] == "RNA || rna") {
             $con =~ s/T/U/;
         }
    }
}

open($in, $ARGV[1]);
while (chomp($line = <$in>)) {
       $miseq = get_seq($con[$_], 80);
       $line =~ s/$name/$miseq/;
       print "$line\n";
}

sub fasta_get
{
        my $in = shift;
        my($l, $name, $seq);

        while ($l = <$in>) {
                if (substr($l, 0, 1) eq '>') {
                        chomp $l;
                        last;
                }
        }
        if ($l eq '') {
                return ();
        }

        $name = substr($l, 1);

        while ($l = <$in>) {
                if (substr($l, 0, 1) eq '>') {
                        seek($in, -length($l), 1);
                        return ($name, $seq);
                }
                chomp $l;
                $seq .= $l 
        }
        return ($name, $seq);
}

sub get_seq
{
    my $seq = shift;
    my $line_len = shift;
    my $seq_len = length $seq;

    unless ($line_len) {
            $line_len = 80;
    }

    for (my $i = 0; $i < $seq_len; $i += $line_len) {
            substr($seq, $i, $line_len);
    }
}

また、3つ目の引数を"RNA"か"rna"と入力した際は、同時に文字列中の"T"を"U"に置換したい。

生じてしまう不具合、考えられること
blast_result.tsvが何も変わらずにそのまま出力されてしまう。
恐らく配列や変数の宣言において、間違いや不足しているものがあるために、思うように動かないのではないか。
また、sub get_seq中の"80"の部分は、以前自身で書いたスクリプトの一部を転用しているので、必要ないかもしれないが、
今回あっても特に問題ないものであると私は思っている。

以上です。プログラムを完全に独学でやっている私には、プロの皆様には常識的なものが欠けている可能性がありますが、
御教授頂けたら幸いです。

Answer 1

10万行位は配列に入れられると思って

(@ARGV == 3) or die "usage: $0 file.fa blast_result.tsv DNAorRNA\n";

($ARGV[2] eq "DNA" || $ARGV[2] eq "dna" || $ARGV[2] eq "RNA" || $ARGV[2] eq "rna") or die "please select DNA(dna) or RNA(rna)\n";

$file = $ARGV[0];

open($fh, $file) or die "cannot open the fasta!!\n";
while (($name, $seq) = fasta_get($fh, $ARGV[2])) {
        $con[$name] = $seq;
}

open($in, $ARGV[1]);
while ($line = <$in>) {
       chomp($line);
       ($name) = split (/\s+/, $line);
       $line =~ s/$name/$con[$name]/;
       print "$line\n";
}

sub fasta_get
{
        my $in = shift;
        my $pm = shift;
        my $seq;
        while (my $l = <$in>) {
            chomp($l);
            if (substr($l, 0, 1) eq '>') {
                $name = substr($l,1);
            }else{
                  $seq = $l;
                  if ($pm eq "RNA" || $pm eq "rna"){
                      $seq =~ s/T/U/g;
                  }
                  return ($name, $seq);
            }
        }
        return ();
}

Comments

[退会済みユーザー]

御多忙のところを態々スクリプトを書き直して頂き、誠に有難う御座いました。
お陰様で実現できました。
スクリプトについては復習したいと思います。

Answer 2

これでは意図したとおりの挙動にはならないだろうという箇所が、非常に沢山見受けられます。また、アルゴリズムも根本的に見直したほうが良さそうです。

まず、

$ARGV[2] == "DNA || dna || RNA || rna";

これは全く意味をなしません。==は数値比較演算子です。"DNA || dna || RNA || rna"を数値化して（つまりゼロとみなして）$ARGV[2]と比較しようとしています。これはおそらく貴方の意図とは全く違うでしょう。たぶん、正規表現を用いて、

$ARGV[2]=~/^(DNA|RNA$)/i;

あたりが妥当な書き方となるでしょう。

 @con = $seq;

$seq には単一の文字列が入っているはずです。にもかかわらずリストとして処理しようとしています。この後のループは、先述の通り演算子の使い方がおかしいことも有り、文字列の置き換えを全く行わないでしょう。

$con =~ s/T/U/;

これでは最初に出てきたTしか置換しません。正しくは、

$con =~ s/T/U/g;

このように、書法上いろいろ不備があることに加え、より根本的な問題として、ID-文字列の対応関係がどこにも保存されること無く読み捨てられているため（$nameを使って置き換えようとしていますが、最初のループに於いて、次々に新しいデータによって書き換えられてしまいますから、役割を果たしていません）、csvファイルの文字列置き換えが機能しようもないという状況になっています。
ここは、連想配列を一つ用意して$id2seq{"4-2013"}のような形で"TGTAAACATCCTTGACTGAAAGCT"を参照できるよう算段せねばならないでしょう。

まずは、デバッガの使い方を覚えて下さい。ステップバイステップで、処理が自分の思っていたとおりに進んでいるか、変数の値をチェックしていって下さい。
次に、リストコンテキストとスカラーコンテキストの違いを理解してください。一部の関数は引数にリストを取ったときとスカラーをとったときとで全く異なる挙動を取ります。
スカラーについては、それを数値として扱うのか文字列として扱うのかをきちんと区別して考えて下さい。文字列比較演算子は==ではなくeqです。

Comments

[退会済みユーザー]

早々の御回答、誠に有難う御座います。
デバッガなるものを知りませんでした。是非活用したいと思います。
$con =~ s/T/U/;の部分は、ループ文中に入れたのでこれでいいつもりでいましたが、勘違いのようでしたね。