perl バイオインフォマティクス

Question

###実現させたいこと
ダウンロードしたmultifastaを一行ずつのものに分割し、NGSによるreadのようなものを作成したい。
先日同じような質問をさせて頂きましたが、今回もそれに似たようなものです。
御教授頂けたら幸いです。

fastaとは生物の核酸あるいはタンパク質の配列を文字列で示したもので、
multifastaは1つのファイルに複数の配列情報を含むもの。
'>'で始まる行はその配列の名称やIDで、その次の行からそれに対応する実際の配列である。
次の'>'で始まる行からは別の配列情報となる。

multifastaの例

gi|1035248882|ref|XM_017038806.1|

GTCGCAACCTTTTGAATCTCCTCGATAAACTGTTACTGCAGAGACGTTTTCCCCAACCAGCCCAGACAGAATGCGCCATGATACCTCAGCTCACCGTCGTAAGGGAAAAGCCCCTCAGAGACGCTCACCGGTGCCAGGTCCGTCCACGTCGACACAACGTTCACCGAGACGCAGCGGGCCTTCAGGACTCCATCATCCAGCATCTACTAAGAAAAGGAAGTTTAGGCCAGGAACCAAGGCCTTGATGGAAATCCGAAAGTACCAGAAAAGCACAGACCTTTTGCTAAGGAAGGCACCATTCTCGCGTCTGGTTCGAGAAGTTTGTCAAAGTTTTTCCCGAGGAGCCTTAAAATGGCAGGTGTATGCTCTCCTGGCTGTGCAGGAGGCTGCAGAGGCTTTTCTTGTCAGACTGTTTGCGGACTCAAACCTGTGTGCCATCCACGCCAAGAGAGTGACCTTGTTCCCTCGTGACATCCAGTTGGCAAGGAGGATCCGTGGAGCAGAAAACCTTTAGAGACAGAGAAAAACACACTGTCTGACATGACGTGCAGCAATGCGAGAACTGTTTCTTTATTTGGTTGGTTACCTTTTGACCTGTTTTCATAGCTGTTGTGTTCATGATGCCTCTCAGTTGTATTATGCAAATAAAGATTTATTTTT

gi|1035248881|ref|XM_008324367.2|

TTCCCGGCACTCAGCCCGGACTTCCTGATTACTTCTTCTTTCTTCCAGCCGGGCGGCAGAGCGGCACATTCCGTCCTCTTCCTCTCCATCACCATACAAGAAGCGATTGGCGTAACTGTTACTGCAGAGACGTTTTCCCCAACCAGCCCAGACAGAATGCGCCATGATACCTCAGCTCACCGTCGTAAGGGAAAAGCCCCTCAGAGACGCTCACCGGTGCCAGGTCCGTCCACGTCGACACAACGTTCACCGAGACGCAGCGGGCCTTCAGGACTCCATCATCCAGCATCTACTAAGAAAAGGAAGTTTAGGCCAGGAACCAAGGCCTTGATGGAAATCCGAAAGTACCAGAAAAGCACAGACCTTTTGCTAAGGAAGGCACCATTCTCGCGTCTGGTTCGAGAAGTTTGTCAAAGTTTTTCCCGAGGAGCCTTAAAATGGCAGGTGTATGCTCTCCTGGCTGTGCAGGAGGCTGCAGAGGCTTTTCTTGTCAGACTGTTTGCGGACTCAAACCTGTGTGCCATCCACGCCAAGAGAGTGACCTTGTTCCCTCGTGACATCCAGTTGGCAAGGAGGATCCGTGGAGCAGAAAACCTTTAGAGACAGAGAAAAACACACTGTCTGACATGACGTGCAGCAATGCGAGAACTGTTTCTTTATTTGGTTGGTTACCTTTTGACCTGTTTTCATAGCTGTTGTGTTCATGATGCCTCTCAGTTGTATTATGCAAATAAAGATTTATTTTT

gi|1035248878|ref|XM_008324364.2|

GTGAGTCGCAACCTTTTGAATCTCCTCGATAAACTGGTACTATTCCGATTCCATTTTTTAAAATACTCGATGCGCAGTTTATTTGTTTTCGAGTTTTATTAATTTTTAAAAAGGTCGGTTATGTTTTAAGGACACAACAGTTACTGCAGAGACGTTTTCCCCAACCAGCCCAGACAGAATGCGCCATGATACCTCAGCTCACCGTCGTAAGGGAAAAGCCCCTCAGAGACGCTCACCGGTGCCAGGTCCGTCCACGTCGACACAACGTTCACCGAGACGCAGCGGGCCTTCAGGACTCCATCATCCAGCATCTACTAAGAAAAGGAAGTTTAGGCCAGGAACCAAGGCCTTGATGGAAATCCGAAAGTACCAGAAAAGCACAGACCTTTTGCTAAGGAAGGCACCATTCTCGCGTCTGGTTCGAGAAGTTTGTCAAAGTTTTTCCCGAGGAGCCTTAAAATGGCAGGTGTATGCTCTCCTGGCTGTGCAGGAGGCTGCAGAGGCTTTTCTTGTCAGACTGTTTGCGGACTCAAACCTGTGTGCCATCCACGCCAAGAGAGTGACCTTGTTCCCTCGTGACATCCAGTTGGCAAGGAGGATCCGTGGAGCAGAAAACCTTTAGAGACAGAGAAAAACACACTGTCTGACATGACGTGCAGCAATGCGAGAACTGTTTCTTTATTTGGTTGGTTACCTTTTGACCTGTTTTCATAGCTGTTGTGTTCATGATGCCTCTCAGTTGTATTATGCAAATAAAGATTTATTTTT
・
・
・ （今回は1つのファイルにfastaが約2000存在）

これを次のように1文字ずつずらして分割したい。ID行部分は連番をふって出力したい

gi|1035248882|ref|XM_017038806.1|-1

GTCGC・・・（>gi|1035248882|ref|XM_017038806.1|の1～80の文字）

gi|1035248882|ref|XM_017038806.1|-2

TCGCA・・・（>gi|1035248882|ref|XM_017038806.1|の2～81の文字）

gi|1035248882|ref|XM_017038806.1|-3

CGCAA・・・（>gi|1035248882|ref|XM_017038806.1|の3～82の文字）
・
・
・

gi|1035248878|ref|XM_008324364.2|-n

・・・ATTTTT （80文字ない行はそのまま出力）

コマンドは "$0 「ファイル名」 20"と入力し、20文字未満のものは出力されないようにした。

###ソースコードを以下に示す。

(@ARGV == 2) or die "usage: $0 file.fa min_contig_length\n";

$file = $ARGV[0];
$min_len = $ARGV[1];

open($fh, $file) or die "cannot open\n";
while (($name, $seq) = fasta_get($fh)) {
        foreach my $con (@seq) { 
         substr($seq, 0, 80);
         seek($seq, -79, 1);
         }
        $name =~ /^(\S+)/;
        $name = $1;
        for (0..$#con) {
                if (length($con[$_]) >= $min_len) {
                        print ">$name-$_\n";
                        $con_name = $name;
                        print_seq($con[$_], 80);
            }
      }
}

sub print_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) {
                print(substr($seq, $i, $line_len) . "\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);
}

###生じた問題点・間違っていると思われる部分
実際に用いたところ、何も出力されなかった。"perl -cw $0"を入力したところ、以下のエラーを確認

Useless use of substr in void context at NGSreads_maker.pl line 11.
Name "main::con_name" used only once: possible typo at NGSreads_maker.pl line 19.
NGSreads_maker.pl syntax OK

Name "・・・の箇所は以前このプログラムを改変して用いていた時から出ていたが、特に気にするものではないはずである。

恐らく9～12行目のforeach文が悪いと思われますが、今の自分では対処方が分かりません。

Comments

[moonphase]

ソースが読みにくいので ``` で囲ってください

[退会済みユーザー]

コード表示の機能を知りませんでした。有り難う御座います

Answer 1

foreach ループの中だけで有効な $con がループの外で使われています。
###追記
理解できました。
foreach 構文を間違って覚えておられます。
次のように使います。

my @arr = ("あいうえお", "かきくけこ", "さしすせそ");

foreach my $item (@arr)
{
	print "($item)¥n";
}

結果は次のようになります。

(あいうえお)
(かきくけこ)
(さしすせそ)

$item は、@arr の要素をすべて含んだ配列になるのではなく、ループの中で @arr の現在の要素が順に代入されます。

Comments

[Zuishin]

あと、これは直接問題には関係ない質問のテクニックなのですが、「バイオインフォマティクス」なる多くの人にとって耳慣れない言葉をタイトルに入れたり、いくつもの専門用語を本文に入れたりすると、それだけで回答者は離れていきます。

問題はバイオインフォマティクスの中になく、Perl プログラミングにあります。
この分野に全く無知の素人でも Perl さえ使えれば回答できます。

ですので、タイトルは「Perl を使ってのデータ解析が失敗します」で十分です。fasta や multifasta の説明もいりません。「次のようなデータ」で意味は通ります。
専門用語の説明をするなら、冒頭ではなく、最後の最後に補足として入れたらいいでしょう。

Answer 2

@con が作成されていないのだから、何もプリントされないのは当たり前です。

該当部分に最低限の変更を加えたコードを記載します。

while (($name, $seq) = fasta_get($fh)) {
	my @con ;
	for ( 0 .. length $seq - 1 ){
		$s = substr $seq, $_, 80 ;
		push @con, $s
	}
	$name =~ /^(\S+)/;
	$name = $1;
	for (0..$#con) {
		if (length($con[$_]) >= $min_len) {
			my $i = $_ + 1 ;
			print ">$name-$i\n";
			$con_name = $name;
			print_seq($con[$_], 80);
		}
	}
}

• 本来では、 push @con, $s の前に長さ判定などを加えたい所ですが、元のコードで後程判定が入る為に割愛しました。
• 追記、連番が 0 からスタートになっていたので、後半にも、ちょっと手を加えました。
• さらに一応追記、index で回しているのに最終値が長さ。orz...

Comments

[退会済みユーザー]

回答有り難う御座います。投稿頂いた方法で解決致しました。
なるほど、@conを作成していませんでしたし、1文字ずつずらすプロセスも（失礼ですが）単純に
pushを用いれば良かったですね。また、連番を1からスタートさせるプロセスも追加頂き有難う御座いました。こういった一手間が分かりやすいプログラム作成に繋がるのだと思います。
御多忙のところを誠に有難う御座いました。

Answer 3

正直データの流れが少し分かりづらいですね。
私の力不足を公言しているようなものですが……。

シンプルに書き直してみました。
オプションは"入力ファイル名 出力ファイル名 単位あたりの文字数 最低文字数"です。
例：test.pl input.txt output.txt 80 20

私の環境ではデータ数が2000個の場合、処理時間は5秒程でした。

use strict;
use warnings;
use utf8;
use open IN    => ":encoding(utf8)";
use open OUT   => ":encoding(utf8)";
binmode STDIN  => ":encoding(utf8)";
binmode STDOUT => ":encoding(utf8)";
binmode STDERR => ":encoding(utf8)";
use Encode qw/encode decode/;

(@ARGV == 4) or die "usage: input_file_name output_file_name chunk min_len\n";
my $input_file_name = $ARGV[0];
my $output_file_name = $ARGV[1];
my $chunk = $ARGV[2];
my $min_len = $ARGV[3];
my $name = '';

open my $INPUT, "<", $input_file_name or die "Can't open $input_file_name : $!";
open my $OUTPUT, ">", $output_file_name or die "Can't open $output_file_name : $!";
while(<$INPUT>) {
	chomp $_;
	$_ =~ s/\s//ig;
	next if (!$_);

	if ($_ =~ /^\>/) {
		$name = $_;
	} else {
		&getFasta($_);
	}
}
close $INPUT;
close $OUTPUT;

sub getFasta()
{
	my ($fasta) = @_;
	my $num = length($fasta)-$chunk;

	return if (length($fasta) < $min_len);
	&printFasta($name,$fasta,0) if ($num < 0);

	for(my $i=0; $i<$num; $i++) {
		&printFasta($name,substr($fasta,$i,$chunk),$i);
	}
}

sub printFasta()
{
	my ($name,$fasta,$i) = @_;
	print $OUTPUT $name."-".($i+1)."\n".$fasta."\n";
}

Comments

[退会済みユーザー]

態々初めからプログラムを作成頂き、誠に有難う御座いました。
このプログラムのほうが置換演算子などの初心者には分かりやすい形になっておりますね。
実は私自身、初心者中の初心者であり元のプログラムも、プロの先生に作成頂いたものを
今回の目的のために一部改変したものであります（申し訳御座いません）。早く自分で全て
書けるようになるためにも、作成頂いたプログラムを見て勉強させて頂きます。
御多忙のところを本当に有難う御座いました。