いつもお世話になっております。自身で書いた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"の部分は、以前自身で書いたスクリプトの一部を転用しているので、必要ないかもしれないが、
今回あっても特に問題ないものであると私は思っている。
以上です。プログラムを完全に独学でやっている私には、プロの皆様には常識的なものが欠けている可能性がありますが、
御教授頂けたら幸いです。
回答2件
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
退会済みユーザー
2017/06/27 07:46