ここから本文です

fastaファイル(fasta.fa)からkeyword.txtで入っている配列名の配列を抽出したい...

アバター

ID非公開さん

2016/7/412:01:45

fastaファイル(fasta.fa)からkeyword.txtで入っている配列名の配列を抽出したいです。

keyward.txt
>TD25TR23901|c0_g1_i5_len=2106_path=[9987:0-346_10211:347-377_10194:378-424_10222:425-439_10208:440-463_10192:464-473_10172:474-497_10207:498-538_10191:539-562_10171:563-722_10120:723-746_10221:747-762_10170:763-786_10206:787-797_10169:798-1001_@10228@!:1002-1287_10112:1288-1310_10205:1311-1390_10074:1391-1391_10075:1392-1394_10076:1395-1412_10165:1413-1537_10115:1538-1538_10116:1539-1601_10117:1602-1602_10118:1603-2105]_[-1,_9987,_10211,_10194,_10222,_10208,_10192,_10172,_10207,_10191,_10171,_10120,_10221,_10170,_10206,_10169,_10228,_10112,_10205,_10074,_10075,_10076,_10165,_10115,_10116,_10117,_10118,_-2]_evgclass=main,okay,match:TD25TR23901|_c0_g1_i4,pct:100/97/.;_aalen=475,67%,complete;
>497-[117_R2.fq]_(paired)_contig_5707_evgclass=noclass,okay;_aalen=44,61%,complete;
以下続く
fasta.fa
>TD25TR23901|c0_g1_i5_len=2106_path=[9987:0-346_10211:347-377_10194:378-424_10222:425-439_10208:440-463_10192:464-473_10172:474-497_10207:498-538_10191:539-562_10171:563-722_10120:723-746_10221:747-762_10170:763-786_10206:787-797_10169:798-1001_@10228@!:1002-1287_10112:1288-1310_10205:1311-1390_10074:1391-1391_10075:1392-1394_10076:1395-1412_10165:1413-1537_10115:1538-1538_10116:1539-1601_10117:1602-1602_10118:1603-2105]_[-1,_9987,_10211,_10194,_10222,_10208,_10192,_10172,_10207,_10191,_10171,_10120,_10221,_10170,_10206,_10169,_10228,_10112,_10205,_10074,_10075,_10076,_10165,_10115,_10116,_10117,_10118,_-2]_evgclass=main,okay,match:TD25TR23901|_c0_g1_i4,pct:100/97/.;_aalen=475,67%,complete;
AAAAAAAAAAAAAAAAA
>497-[117_R2.fq]_(paired)_contig_5707_evgclass=noclass,okay;_aalen=44,61%,complete;
GGGGGGGGGGGGGG
以下続く

今使っているスクリプトでは
配列名で[ ]など様々な記号が入っているせいかうまく認識してくれません。
どこを直せばいいのでしょうか?
もしくはrubyなどでこれを解決できますでしょうか?

また
>A length=200
などAの後にスペースが入っているため、
>Aというkeywordでは引っかからないです。
これを解決するようなもお願いします

perlスクリプト

##################### Please make sure the name of your files.
my $nameList = "keyword.txt"; ## Keyword list of focal records.
my $db = "fasta.fa"; ## Data base file (Fasta format).
my $out = "test.fa"; ## Output file
#####################




## name list
open(INFILE,"$nameList") or die "Cannot find the keyword list, $nameList\n";
my @inFileCont = <INFILE>;
close INFILE;

## db
open(DB, "$db") or die "Cannot find the database, $db\n";
my @dbFileCont = <DB>;
close DB;


### Main program ####

my ($ref_fasHash, $ref_nameArray, $seqLength) = &fasHashReader(\@inFileCont);
my ($ref_fasHashDB, $ref_nameArrayDB, $seqLengthDB) = &fasHashReader(\@dbFileCont);


# Printout
open(OUT, ">$out");
foreach my $name (@$ref_nameArray){
if(my @tempNameLineDB = grep(/$name/, @$ref_nameArrayDB)){
print "FOUND $name\n";
my $nameLineDB = join('',@tempNameLineDB);
print OUT "$nameLineDB\n";

#my $seqTOed = &turnover($$ref_fasHashDB{$nameLineDB});
#print OUT "$seqTOed\n";
print OUT "$$ref_fasHashDB{$nameLineDB}\n";

} else {
print "NOT FOUND $name\n";
}
}
close OUT;

exit;




#### サブルーチン ####

sub turnover {
my ($seqTO) = @_;
$seqTO =~ s/(.{80})/$1\n/g;
#$seqTO = "|".$seqTO."|\n\n\n";
return($seqTO);
}


# Fasta 形式をハッシュに読み込む Subroutine
sub fasHashReader
{
my $name = "";
my $seq = "";
my $seqLength = "";
my $stone = 0;
my @nameArraySR = ( );
my %fasHashSR = ( );

my ($inFileContSR) = @_;

foreach my $inFileLine (@$inFileContSR) {
chomp($inFileLine);
if($inFileLine =~ /^>/){
if($stone == 0){
$name = $inFileLine;
$stone = 1;
} else {
$fasHashSR{$name} = $seq;
push(@nameArraySR,$name);
$name = "$inFileLine";
$seq = "";
}
} else {
$seq .= $inFileLine;
}
}
$fasHashSR{$name} = $seq;
$seqLength = length($seq);
push(@nameArraySR,$name);

# ハッシュと配列をメインプログラムに返す.
return (\%fasHashSR,\@nameArraySR,\$seqLength);

}

閲覧数:
162
回答数:
1

違反報告

ベストアンサーに選ばれた回答

pis********さん

2016/7/420:25:42

> 配列名で[ ]など様々な記号が入っているせいかうまく認識してくれません。

このエラーメッセージ自体は
Invalid [] range "7-3" in regex; marked by

> if(my @tempNameLineDB = grep(/$name/, @$ref_nameArrayDB)){

ここで $name に正規表現が含まれてるのが原因なので、quotemeta すればとりあえずは回避できます。

if(my @tempNameLineDB = grep(quotemeta($name), @$ref_nameArrayDB)){

これで欲しい結果が出るかどうかまでは説明されてないので判断できません。ですが、fasta.fa ってアミノ酸配列の FASTA ですよね?
それなら CPAN にも色々公開されているので、そういうモジュールを上手く利用すると便利だと思います。

  • アバター

    質問者

    ID非公開さん

    2016/7/510:06:06

    ありがとうございます。教えていただいた方法でできました。
    このスクリプトでは以下のスペース入っているものは解決しませんでした。こちらも解決する方法があれば教えていただきたいです。

    >TCONS_00002679 gene=XLOC_001745
    >TCONS_00002680 gene=XLOC_001746
    >TCONS_00002681 gene=XLOC_001747
    ...
    以下続く

    アウトプットは
    >TCONS_00002679 gene=XLOC_001745>TCONS_00002680 gene=XLOC_001746>TCONS_00002681 gene=XLOC_001747>TCONS_00002682 gene=XLOC_001748>TCONS_00002683 gene=XLOC_001749>TCONS_00002684 gene=XLOC_001750>TCONS_00002685 gene=XLOC_001751

  • その他の返信を表示

返信を取り消しますが
よろしいですか?

  • 取り消す
  • キャンセル

みんなで作る知恵袋 悩みや疑問、なんでも気軽にきいちゃおう!

Q&Aをキーワードで検索:

Yahoo! JAPANは、回答に記載された内容の信ぴょう性、正確性を保証しておりません。
お客様自身の責任と判断で、ご利用ください。
本文はここまでです このページの先頭へ

「追加する」ボタンを押してください。

閉じる

※知恵コレクションに追加された質問は選択されたID/ニックネームのMy知恵袋で確認できます。

不適切な投稿でないことを報告しました。

閉じる