BLAST_out_perl

perl BLASTの出力結果を加工する 22 January 2021 †

FASTA形式の複数（16,000以上）の塩基配列がショウジョウバエゲノム上でどちらの向きか（プラスかマイナスか）を判定したい

NCBI BLAST とちょっとした perl プログラムを使って、ざっと判定をしたので、そのメモ

この作業内容では、bioinformatics とは言えませんね (笑)

FlyBase BLAST では、複数の塩基配列を対象とした BLAST が行えない
DDBJ では、受け入れ可能ファイルサイズが小さい
との使いづらさがあったので、NCBI BLAST を使った

fasta ファイルの加工 †

NCBI BLAST の Nucleotide BLAST では、受け入れるファイルサイズが最大 1,000,000 だったので*1、今回は 16,271 個を 1700 個ずつに分割*2

$ head -n  3400 sequence.fasta > sequence1.fasta
$ head -n  6800 sequence.fasta | tail -n 3400 > sequence2.fasta
$ head -n 10200 sequence.fasta | tail -n 3400 > sequence3.fasta
$ head -n 13600 sequence.fasta | tail -n 3400 > sequence4.fasta
$ head -n 17000 sequence.fasta | tail -n 3400 > sequence5.fasta
$ head -n 20400 sequence.fasta | tail -n 3400 > sequence6.fasta
$ head -n 23800 sequence.fasta | tail -n 3400 > sequence7.fasta
$ head -n 27200 sequence.fasta | tail -n 3400 > sequence8.fasta
$ head -n 30600 sequence.fasta | tail -n 3400 > sequence9.fasta
$ tail -n  1942 sequence.fasta > sequence10.fasta

区切りが正しいか元データと手動で照合した

↑

NCBI BLAST を行う †

NCBI BLAST に sequence1.fasta, ..., sequence10.fasta を投げた

Organism を Drosophila melanogaster (taxid:7227) に指定

「Show results in a new window」のチェックを入れれば、次から次へと10回 BLAST検索ができる。順序よく待てばよい。といっても、大して待たされない。

Download All より、text をダウンロード。
ファイル名（blastout1.txt, ..., blastout10.txt）を付けて保存。

出力結果ファイルを結合

$ cat blastout1.txt blastout2.txt blastout3.txt blastout4.txt ［改行せず］
blastout5.txt blastout6.txt blastout7.txt blastout8.txt  ［改行せず］
blastout9.txt blastout10.txt > NCBI_BLAST_all.txt

↑

ゲノムの Sequence ID †

NCBI BLAST は検索の対象がゲノム配列だけに絞れないので（この点が Flybase と大きく違って手間になる）、ショウジョウバエ（D. melanogaster）ゲノムの Sequence ID を手掛りにして、結果を絞り込む必要がある

https://www.ncbi.nlm.nih.gov/genome/47?genome_assembly_id=204923 の Replicon Info にゲノム配列に使われている Sequence ID の情報がある

ショウジョウバエ（D. melanogaster）ゲノムの Sequence ID は以下の通り

X:	AE014298.5
2L:	AE014134.6
2R:	AE013599.5
3L:	AE014296.5
3R:	AE014297.3
4:	AE014135.4
Y:	CP007106.1

今回は、核ゲノムのみ扱うのでこれだけでよい

↑

NCBI BLAST から必要なものを抽出する perl プログラム †

今回は、query配列がゲノム上でプラス鎖かマイナス鎖のどちらかを知りたいだけなので、それだけを取り出し、タブ区切りテキストを標準出力に吐き出すことにした

また、念のために、染色体（腕）とゲノムの Sequence ID も同じ行に出力させる

なお、print の部分を push などに変えれば、更なる加工ができるはず

↑

NCBI BLAST から必要なものを抽出 †

perl プログラムを使って、NCBI BLAST の結果から必要なものを抽出

$ perl extractfromblast.pl NCBI_BLAST_all.txt > NCBI_BLAST_exstract.txt

↑

perl プログラム：extractfromblast.pl †

#!/usr/bin/perl  
use strict;  
use warnings;  
  
my $myinput = $ARGV[0];  
open(BLASTOUT,$myinput);  
  
my @mydata = <BLASTOUT>;  
close BLASTOUT;  
  
# Print Header	  
print "Query\tChromosome\tGenomeID\tStrand\n";	  
#	  
  
my @extractedlines = ();  
for (my $myi = 0; $myi < scalar(@mydata); $myi++) {  
    my $line = $mydata[$myi];  
    if ($line =~ m/Query #/g ) {  
#	  
#	  
# Example  
# Query #1: AB*****.1 Query ID: ***|Query_***** Length: 60  
#	  
	my @splitline = split(/ +/, $line);  
	my $queryID = $splitline[2];  
	print "\n$queryID\t";  
    } elsif ($line =~ m/>/g ) {  
#	  
# Example  
# >Drosophila melanogaster chromosome X  
# Sequence ID: AE014298.5 Length: 23542271   
# Range 1: 21318987 to 21319046  
#   
# Score:111 bits(60), Expect:6e-24,   
# Identities:60/60(100%),  Gaps:0/60(0%), Strand: Plus/Minus  
#	  
	my $line2 = $mydata[$myi + 1];  
	my @splitline = split(/ +/, $line2);  
	my $sequenceID = $splitline[2];  
  
	my $line3 = $mydata[$myi + 5];  
	@splitline = split(/ +/, $line3);  
	my $strand = $splitline[3];  
#  
# Sequence IDs of Drosophila melanogaster genome:   
#  
# X  
# AE014298.5  
#   
# 2L  
# AE014134.6  
#   
# 2R  
# AE013599.5  
#   
# 3L  
# AE014296.5  
#   
# 4  
# AE014135.4  
#   
# 3R  
# AE014297.3  
#   
# Y  
# CP007106.1  
#   
	if ($sequenceID =~ m/AE014298.5/ ) { # X chromosome  
	    print "X\t$sequenceID\t";  
	    print "$strand\n";  
	} elsif ($sequenceID =~ m/AE014134.6/ ) { # 2L  
	    print "2L\t$sequenceID\t";  
	    print "$strand\n";  
	} elsif ($sequenceID =~ m/AE013599.5/ ) { # 2R  
	    print "2R\t$sequenceID\t";  
	    print "$strand\n";  
	} elsif ($sequenceID =~ m/AE014296.5/ ) { # 3L  
	    print "3L\t$sequenceID\t";  
	    print "$strand\n";  
	} elsif ($sequenceID =~ m/AE014297.3/ ) { # 3R  
	    print "3R\t$sequenceID\t";  
	    print "$strand\n";  
	} elsif ($sequenceID =~ m/AE014135.4/ ) { # 4  
	    print "4\t$sequenceID\t";  
	    print "$strand\n";  
	} elsif ($sequenceID =~ m/CP007106.1/ ) { # Y  
	    print "Y\t$sequenceID\t";  
	    print "$strand\n";  
	}  
    }  
}

print の部分を push などに変えるのは、例えば：

print "X\t$sequenceID\t";  
print "$strand\n";

を

push(@extractedlines, "X\t$sequenceID\t" );
push(@extractedlines, "$strand");

などとする。そのまま出力するだけならば

for (my $myj=0; $myj<scalar(@extractedlines); $myj++) {
    print $extractedlines[$myj];
}

などとすればよい。更に、加工／解析するのならば $extractedlines[$i] を上手く利用するようにプログラムを作る、とか。

↑

BLAST の結果、複数のゲノム染色体がヒットする場合は…… †

このプログラムは BLAST の結果の出力ファイルの上から処理していくので、確率が大きい（余り確からしくない）ものには、Query名が付与されない。そこで、そういったデータを後から削る必要がある。
BLAST の確率が小さい（より確からしい）ものは正しい出力となり、それ以外は１列左にずれて出力するので、 Query の列が染色体名（X、2L、2R、3L、3R、4、Y）になる。
そこで、Query の列（＝行頭）が染色体名（X、2L、2R、3L、3R、4、Y）のものを、ふたつめにヒットしたものとして削除する。

$ grep -v ^X NCBI_BLAST_exstract.txt | grep -v ^2 | grep -v ^3 | grep -v ^4 | grep -v ^Y > NCBI_BLAST_exstract_fin.txt

最後に、余分な改行をエディタを使って、手動で削除する。

なお、BLAST の結果が「ない」ものがあるので、それは個別に見る必要あり。

Today:1

Yesterday:3

Total:3609 since 22 January 2021

Menu

最新の10件

perl BLASTの出力結果を加工する 22 January 2021 †

fasta ファイルの加工 †

NCBI BLAST を行う †

ゲノムの Sequence ID †

NCBI BLAST から必要なものを抽出する perl プログラム †

NCBI BLAST から必要なものを抽出 †

perl プログラム：extractfromblast.pl †

BLAST の結果、複数のゲノム染色体がヒットする場合は…… †

BLAST_out_perl

Menu

最新の10件

perl BLASTの出力結果を 加工する 22 January 2021 †

fasta ファイルの加工 †

NCBI BLAST を行う †

ゲノムの Sequence ID †

NCBI BLAST から必要なものを抽出する perl プログラム †

NCBI BLAST から必要なものを抽出 †

perl プログラム：extractfromblast.pl †

BLAST の結果、複数のゲノム染色体がヒットする場合は…… †

perl BLASTの出力結果を加工する 22 January 2021 †