バイオインフォマティクス関連のプログラム

Question

###前提・実現したいこと
ダウンロードしたmultifastaを一行ずつのものに分割し、NGSによるreadのようなものを作成したい。
御教授頂けたら幸いです。

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

multifastaの例（今回の配列部分は、A,G,C,Tの4文字のみから成る）

gi|1035248882|ref|XM_017038806.1|

GTCGCAACCTTTTGAATCTCCTCGATAAACTGTTACTGCAGAGACGTTTTCCCCAACCAGCCCAGACAGAATGCGCCATGATACCTCAGCTCACCGTCGTAAGGGAAAAGCCCCTCAGAGACGCTCACCGGTGCCAGGTCCGTCCACGTCGACACAACGTTCACCGAGACGCAGCGGGCCTTCAGGACTCCATCATCCAGCATCTACTAAGAAAAGGAAGTTTAGGCCAGGAACCAAGGCCTTGATGGAAATCCGAAAGTACCAGAAAAGCACAGACCTTTTGCTAAGGAAGGCACCATTCTCGCGTCTGGTTCGAGAAGTTTGTCAAAGTTTTTCCCGAGGAGCCTTAAAATGGCAGGTGTATGCTCTCCTGGCTGTGCAGGAGGCTGCAGAGGCTTTTCTTGTCAGACTGTTTGCGGACTCAAACCTGTGTGCCATCCACGCCAAGAGAGTGACCTTGTTCCCTCGTGACATCCAGTTGGCAAGGAGGATCCGTGGAGCAGAAAACCTTTAGAGACAGAGAAAAACACACTGTCTGACATGACGTGCAGCAATGCGAGAACTGTTTCTTTATTTGGTTGGTTACCTTTTGACCTGTTTTCATAGCTGTTGTGTTCATGATGCCTCTCAGTTGTATTATGCAAATAAAGATTTATTTTT

gi|1035248881|ref|XM_008324367.2|

TTCCCGGCACTCAGCCCGGACTTCCTGATTACTTCTTCTTTCTTCCAGCCGGGCGGCAGAGCGGCACATTCCGTCCTCTTCCTCTCCATCACCATACAAGAAGCGATTGGCGTAACTGTTACTGCAGAGACGTTTTCCCCAACCAGCCCAGACAGAATGCGCCATGATACCTCAGCTCACCGTCGTAAGGGAAAAGCCCCTCAGAGACGCTCACCGGTGCCAGGTCCGTCCACGTCGACACAACGTTCACCGAGACGCAGCGGGCCTTCAGGACTCCATCATCCAGCATCTACTAAGAAAAGGAAGTTTAGGCCAGGAACCAAGGCCTTGATGGAAATCCGAAAGTACCAGAAAAGCACAGACCTTTTGCTAAGGAAGGCACCATTCTCGCGTCTGGTTCGAGAAGTTTGTCAAAGTTTTTCCCGAGGAGCCTTAAAATGGCAGGTGTATGCTCTCCTGGCTGTGCAGGAGGCTGCAGAGGCTTTTCTTGTCAGACTGTTTGCGGACTCAAACCTGTGTGCCATCCACGCCAAGAGAGTGACCTTGTTCCCTCGTGACATCCAGTTGGCAAGGAGGATCCGTGGAGCAGAAAACCTTTAGAGACAGAGAAAAACACACTGTCTGACATGACGTGCAGCAATGCGAGAACTGTTTCTTTATTTGGTTGGTTACCTTTTGACCTGTTTTCATAGCTGTTGTGTTCATGATGCCTCTCAGTTGTATTATGCAAATAAAGATTTATTTTT

gi|1035248878|ref|XM_008324364.2|

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

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

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

GTCGCAACCTTTTGAATCTCCTCGATAAACTGTTACTGCAGAGACGTTTTCCCCAACCAGCCCAGACAGAATGCGCCATGATACCTC（80文字ずつ）

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

AGCTCACCGTCGTAAGGGAAAAGCCCCTCAGAGACGCTCACCGGTGCCAGGTCCGTCCACGTCGACACAACGTTCACCGAGACGCAG（80文字ずつ）

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

CGGGCCTTCAGGACTCCATCATCCAGCATCTACTAAGAAAAGGAAGTTTAGGCCAGGAACCAAGGCCTTGATGGAAATCCGAAAGTA（80文字ずつ）
・
・
・

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

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

コマンドは、"$0 ファイル名 80 20"と入力

###発生している問題・エラーメッセージ

gi|664708653|ref|XM_008513042.1|-47

TGTTCAGCATCAGCTCCTATTTACGCAACTAAATAAACGATCAAATAAAGTTTA

gi|742122539|ref|XM_010840600.1|-8

TTCAATCCCTGGGTCAGGAAGATCCTCTGGAGAAAGGAATGGCAATCCACTCCAGTATTCTTGCT

gi|1016637223|ref|XM_016210042.1|-14

TATAACTGCTTTGTGGAGAGATAATAAAAATAATAATTCTTGAATGAGAAAA

gi|742122543|ref|XM_010840601.1|-9

GGATGGGGGGACACATGTATACCTGTGACTGATTCATGTTGATGTATGAAAAAACA
・
・
・
1つのfastaから1行しか出力されておらず、その行がどのような基準で出力されているのか理解できない。
更にgrepコマンドで探しても対応するfastaにその配列が存在しているかが疑わしい。

エラーメッセージ
特になし

###該当のソースコード
プログラムを御教授頂いた先生に感謝致します。頂いたプログラムを一部改変させて頂いたものが
本プログラムです。

#! /usr/bin/perl

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

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

open($fh, $file) or die "cannot open\n";
while (($name, $seq) = fasta_get($fh)) {
        @con = split(/[AGCT]{$bp_len}/, $seq);
        $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);
}

###試したこと
sub fasta_getを本分野の参考書にあるような一般的なファイル全体を読み込むような
ものに書き換えて実行（上記のものは1行ずつ処理するものであるはずので）したが、メモリの問題（？）で動かなかった。


###補足情報(言語/FW/ツール等のバージョンなど)
perl v5.18.2

Comments

[退会済みユーザー]

俺らにもわかるような質問にしていただけませんでしょうか？
multifastaもNGSもreadのようなものも　何のことかわからないのです。
プログラムが何をしてるのかはわかりますが、どんなデータをどうしたいのかさっぱりです。

[退会済みユーザー]

質問内容に詳細を追加したいと思います。有難う御座いました。

Answer 1

@con = split(/[AGCT]{$bp_len}/, $seq); 

を

@con = $seq =~ /[AGCT]{$bp_len}/g; 

でどうでしょうか。

Comments

[退会済みユーザー]

解決致しました。
なるほど、改変前のものがsplitを用いていたので、それにずっと拘っていたことがダメでした。
このプログラムではsplitを用いなくても分割されますね。
御多忙のところを有難う御座いました。