###実現させたいこと
ダウンロードした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文が悪いと思われますが、今の自分では対処方が分かりません。
回答3件
あなたの回答
tips
プレビュー