###前提・実現したいこと
ダウンロードした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
回答1件
あなたの回答
tips
プレビュー