質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

92.01%

  • Perl

    275questions

    Perlは多目的に使用される実用性が高い動的プログラミング言語のひとつです。

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

解決済

回答 1

投稿 2016/07/16 12:27 ・編集 2016/07/16 14:02

  • 評価
  • クリップ 0
  • VIEW 201

y-fk

score 3

前提・実現したいこと

ダウンロードした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

  • 気になる質問をクリップする

    クリップした質問に回答があった場合に通知・メールを受け取ることができます。

    クリップした質問はマイページの「クリップ」タブからいつでも見ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正の依頼

  • noricyan2

    2016/07/16 12:43

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

    キャンセル

  • y-fk

    2016/07/16 13:17

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

    キャンセル

回答 1

checkベストアンサー

0

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


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


でどうでしょうか。

投稿 2016/07/16 15:30

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    以下のような回答は評価を下げられます

    • 間違っている回答
    • 質問の回答になっていない投稿
    • 不快な投稿

    評価を下げる際はその理由をコメントに書き込んでください。

  • 2016/07/16 16:16

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

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

ただいまの回答率

92.01%

関連した質問

同じタグがついた質問を見る

  • Perl

    275questions

    Perlは多目的に使用される実用性が高い動的プログラミング言語のひとつです。

閲覧数の多いPerlの質問