MJさむらいの実験室

ＦＡＳＴＡシークエンスのいっぱい入ったファイルの中から、ある特定のＦＡＳＴＡシークエンスを抽出したい時につかえるプログラムを書いてみました。

ここまで書くには結構勉強した。

ＦＡＳＴＡフォーマット
＞１
ＡＴＧＡＴＡＴＧＡＧＧＡＴＧＣＧＴＡＧＴＡ
＞２
ＡＡＡＡＡＴＴＴＴＧＧＧＧＣＣＣＣＣ
＞３
ＴＴＴＴＴＣＣＣＧＴＧＴＡＧＴＧＡＴＧＴＧＴＣＧＴＧＣＴＧＡＴＣＧＴＡＣＧＴＣＧ

という遺伝子が入ったファイルがある。

その中から３のみを抽出したいときは、
＞３
と書かれたファイルを用意して下のプログラムを実行すると。

ＦＡＳＴＡフォーマットのファイルの中から
＞３
ＴＴＴＴＴＣＣＣＧＴＧＴＡＧＴＧＡＴＧＴＧＴＣＧＴＧＣＴＧＡＴＣＧＴＡＣＧＴＣＧ

だけを抽出してくれる。

今５万個とか扱っているのでこれを勉強して書く手間と５万個コピペする手間を考えたら
これを書くほうが楽チン。

プログラムの内容は何でもいいから目的が果たせたらいいという代物なので
褒められたものではないだろう。

ここから下がプログラム

#!/usr/bin/perl -w
use strict;
use warnings;
use POSIX;

# enter the fastafile to hash.
my ($fastafile) = @ARGV;
open FASTA, "<$fastafile";
my %hash=(); # initializes a hash

while (<FASTA>)
{
 if ($_ =~ /^>/)
 {

        my $header = $_;
        $header =~ s/\s//g;
        my $read_id = $_;
      
        $hash{$header}{name}=$read_id;
        my $line = <FASTA>; 
        $hash{$header}{sequence}= $line;
 }
}
close FASTA;



my $counter = 0;

my $name = <STDIN>;

open (FASTANAME, $name);


while (<FASTANAME>)
{

 my $filename =  $_;

  $filename =~ s/\s//g;
my $header = $filename;
if (exists $hash{"$header"})
  {
    $counter = $counter + 1;
            print "$hash{$header}{name}";
        print "$hash{$header}{sequence}";
  }
}
close (FASTANAME);

print "$counter fasta sequences are here\n";