Go to Qiita Advent Calendar Top
LoginSignup
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@yoho

コロナウィルス の配列を比較するために

Posted at

はじめに

コロナウィルスが猛威を振るっています。ウィルスゲノムにどんどん変異が蓄積しているとのことを聞きますが、どれくらい変異が入るものなんでしょうか。

コロナウィルスの配列データ(アクセッション番号 NC_045512)のDNA配列の一部を、NCBIのウェブサイトでblastnにかけ、ヒットした配列を100件程度含むGenBankファイルを入手しました。

この記事では、複数のエントリを含むGenBankファイルから、情報を抜き出すperlスクリプトを記述します。

perlでGenBankファイルをエントリごとに処理する

マルチGenBankファイルでは、LOCUSから始まって//までが１つのエントリです。

コロナウィルス についてはこの3ヶ月ほどでたくさんのデータが登録されており、いつ、どこで単離されたものであるかが重要になりそうです。いつ、どこで、に当たるデータはsourceのところに書いてあるようです。

まずはマルチGenBankファイルを読み込んで、1つずつ処理できるようにしたいと思います。

001.pl
open (inputdata,$ARGV[0]);#ファイルのオープン
$flag1 = 0;#初期化(この記述は必要ではない)
$flag2 = 0;#初期化(この記述は必要ではない)
$genbankText = "";#初期化(この記述は必要ではない)
while (<inputdata>){
    $genbankText .= $_;#データを取っておく
    if($_ =~ /     source          /){#source行にたどり着いたらフラグをあげる
        $flag1 = 1;
        $header .= $_;
        next;
    }
    if( $flag1 == 1 && $_ !~ /                     /){#source行にたどり着いてから、先頭に21個のスペースがない行にたどり着いたら
        $flag2 = 1;
    }
    if($flag2 != 1){
        $header .= $_;#sourceデータまでを記録
    }  
    if($_ =~ /\/\/\n/){
        #//までたどり着いた時の処理
    } 
}

冒頭で変数の初期化をしていますが、perlでは必要ないはずですが、しないと気持ちが悪いので初期化します。

これで、1つのGenBankエントリ全体を、変数$genbankTextに、GenBankエントリのヘッダー部分+sourceの部分を変数$headerに入れることができました。

いくつかサブルーチンを書きます。

GenBankエントリからDNA配列を書き出すサブルーチン

002.pl
# GenBankエントリを渡すと、含まれるDNA配列を大文字にして返す(改行などは含まない)。
sub dnaSequence{
    my $data = @_[0];
    $data =~ m/ORIGIN (.+)\/\//s;
    $seq = $1;
    $seq =~ s/\d//g;#数字を取り除く
    $seq =~ s/\s//g;#空白文字を取り除く
    return uc($seq);#大文字にする。
}

ヘッダー部分からlocusデータ(アクセッション番号)を返すサブルーチン。

003.pl
sub locus_data{
    my $data = @_[0];
    $data =~ /LOCUS       (.+?)\s/;
    return $1;
}

ヘッダー部分からDEFINITIONの記述を取り出すサブルーチン。
DEFINITION部分の記述は複数行に渡る場合があり、ACCESSIONで始まる行の直前までをDEFINITIONの記述であるとみなす。複数のスペースがあれば、スペース1つに変換する。

004.pl
sub definition{
    my $data = @_[0];
    $data =~ m/DEFINITION  (.+)\nACCESSION   /s;
    my $def = $1;
    while($def =~ s/\n/ /g){}
    while($def =~ s/  / /g){}
    return $def;
}

ヘッダー情報から、各種情報を標準出力へprintするようにする。

005.pl
sub infoFromHeader{
    my $data = @_[0];
    #Locusの出力
    $data =~ /LOCUS       (.+?)\s/;
    print $1."\t";
    #DEFINITIONの出力
    print definition($data)."\t";
    #authorsの出力
    $data =~ m/  AUTHORS   (.+?)\n  TITLE     /s;
    my $authors = $1;
    while($authors =~ s/\n/ /g){}
    while($authors =~ s/  / /g){}
    print $authors."\t";
    #  JOURNALの出力(所属機関？？)
    $data =~ m/  JOURNAL   Submitted \(.+\)(.+)\nCOMMENT     /s;
    my $affiliation = $1;
    while($affiliation =~ s/\n/ /g){}
    while($affiliation =~ s/  / /g){}
    print $affiliation."\t";
    #isolateの出力
    if($data =~ /                     \/isolate=\"(.+)\"/){
        print $1."\t";
    }else{
        print "\t";
    }
    #countryの出力
    if($data =~ /                     \/country=\"(.+)\"/){
        print $1."\t";
    }else{
        print "\t";
    }
    #collection_dateの出力
    if($data =~ /                     \/collection_date=\"(.+)\"/){
        print $1."\t";
    }else{
        print "\t";
    }
    
}

配列データにNが含まれているか調べるサブルーチン

006.pl
sub containAnyN{
    my $data = @_[0];
    $data =~ m/ORIGIN (.+)\/\//s;
    if($1 !~ /n/){
        return 0;#Nなし
    }else{
        return 1;#Nあり
    }
}

DEFINITION部分にpartialと書いてあるかどうか調べるサブルーチン。完全には決まっていない配列は使いたくないので。

007.pl
sub isPartial{
    my $data = @_[0];
    my $def = definition($data);
    if($def =~ /partial/){
        return 1;
    }else{
        return 0;
    } 
}

//までたどり着いた時の処理。配列にNが含まれておらず、かつpartialな配列でない場合に、標準出力に各種データを、ファイルにGenBankファイルを書き出す。

008.pl

    if($_ =~ /\/\/\n/){
        $containAnyN = containAnyN($genbankText);
        $isPartial = isPartial($header);
        $dnaSequence = dnaSequence($genbankText);
        if($containAnyN == 0 && $isPartial == 0){#partialでなくNもない
            #標準出力へ書き出し
            infoFromHeader($header);
            print "$dnaSequence";
            print "\n";

            #ファイルに書き出し
            $locus = locus_data($header);
            open (resultFile,">$locus".".gb");
            print resultFile $genbankText;
            close(resultFile);#ファイルを閉じる。
        }
        
        #初期化
        $dnaSequence = "";
        $genbankText = "";
        $header = "";
        $flag1 = 0;
        $flag2= 0;
    }

コロナウィルス のゲノム配列は30 kb弱。ということで1つの配列をエクセルのセル1つに入れることができる。

以上をまとめる

009.pl
open (inputdata,$ARGV[0]);
$flag1 = 0;
$flag2= 0;
$genbankText = "";
while (<inputdata>){
    $genbankText .= $_;
    
    if($_ =~ /     source          /){
        $flag1 = 1;
        $header .= $_;
        next;
    }
    
    if( $flag1 == 1 && $_ !~ /                     /){
        $flag2 = 1;
    }

    if($flag2 != 1){
        $header .= $_;
    }
    
    if($_ =~ /\/\/\n/){
        $containAnyN = containAnyN($genbankText);
        $isPartial = isPartial($header);
        $dnaSequence = dnaSequence($genbankText);
        if($containAnyN == 0 && $isPartial == 0){#partialでなくNもない
            infoFromHeader($header);
            
            print "$dnaSequence";
            print "\n";
            
            $locus = locus_data($header);
            open (resultFile,">$locus".".gb");
            print resultFile $genbankText;
            close(resultFile);#ファイルを閉じる。
        }
        
        #初期化
        $dnaSequence = "";
        $genbankText = "";
        $header = "";
        $flag1 = 0;
        $flag2= 0;
    }
    
}
sub isPartial{
    my $data = @_[0];
    my $def = definition($data);
    if($def =~ /partial/){
        return 1;
    }else{
        return 0;
    }
    
}

sub containAnyN{
    my $data = @_[0];
    $data =~ m/ORIGIN (.+)\/\//s;
    if($1 !~ /n/){
        return 0;#Nなし
    }else{
        return 1;#Nあり
    }
}
sub dnaSequence{
    my $data = @_[0];
    $data =~ m/ORIGIN (.+)\/\//s;
    $seq = $1;
    $seq =~ s/\d//g;#数字を取り除く
    $seq =~ s/\s//g;#空白文字を取り除く
    return uc($seq);
}

sub locus_data{
    my $data = @_[0];
    $data =~ /LOCUS       (.+?)\s/;
    return $1;
}

sub definition{
    my $data = @_[0];
    $data =~ m/DEFINITION  (.+)\nACCESSION   /s;
    my $def = $1;
    while($def =~ s/\n/ /g){}
    while($def =~ s/  / /g){}
    return $def;
}


sub infoFromHeader{
    my $data = @_[0];
    #Locusの出力
    $data =~ /LOCUS       (.+?)\s/;
    print $1."\t";
    #DEFINITIONの出力
    print definition($data)."\t";
    #authorsの出力
    $data =~ m/  AUTHORS   (.+?)\n  TITLE     /s;
    my $authors = $1;
    while($authors =~ s/\n/ /g){}
    while($authors =~ s/  / /g){}
    print $authors."\t";
    #  JOURNALの出力(所属機関？？)
    $data =~ m/  JOURNAL   Submitted \(.+\)(.+)\nCOMMENT     /s;
    my $affiliation = $1;
    while($affiliation =~ s/\n/ /g){}
    while($affiliation =~ s/  / /g){}
    print $affiliation."\t";
    #isolateの出力
    if($data =~ /                     \/isolate=\"(.+)\"/){
        print $1."\t";
    }else{
        print "\t";
    }
    #countryの出力
    if($data =~ /                     \/country=\"(.+)\"/){
        print $1."\t";
    }else{
        print "\t";
    }
    #collection_dateの出力
    if($data =~ /                     \/collection_date=\"(.+)\"/){
        print $1."\t";
    }else{
        print "\t";
    }
    
}

できました。

あとはperlの第一引数にGenBankファイルを指定して実行するだけです。ファイルに書き出すなら、>でリダイレクトですね。

出力したデータはエクセルに貼り付けて使う予定です。

これで複数のエントリを含むGenBankファイルから情報を抜き出すことが簡単にできるようになりました。

1
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?