BLAST_out_perl
をテンプレートにして作成
開始行:
* perl BLASTの出力結果を 加工する 22 January 2021 [#w1abe...
FASTA形式の複数(16,000以上)の塩基配列が
ショウジョウバエゲノム上でどちらの向きか
(プラスかマイナスか)を判定したい
[[NCBI BLAST>https://blast.ncbi.nlm.nih.gov/Blast.cgi]]
とちょっとした perl プログラムを使って、
ざっと判定をしたので、そのメモ
この作業内容では、bioinformatics とは言えませんね (笑)
- [[FlyBase BLAST>http://flybase.org/blast/]]
では、複数の塩基配列を対象とした BLAST が行えない~
- [[DDBJ>http://blast.ddbj.nig.ac.jp/blastn?lang=en]]
では、受け入れ可能ファイルサイズが小さい~
との使いづらさがあったので、NCBI BLAST を使った
** fasta ファイルの加工 [#xfa2bd03]
[[NCBI BLAST>https://blast.ncbi.nlm.nih.gov/Blast.cgi]] の
Nucleotide BLAST では、受け入れるファイルサイズが最大 1,0...
だったので((DDBJ は最大 100,000 とひと桁小さく作業が面倒...
今回は 16,271 個を 1700 個ずつに分割((FASTA形式は2行で1...
$ head -n 3400 sequence.fasta > sequence1.fasta
$ head -n 6800 sequence.fasta | tail -n 3400 > sequence...
$ head -n 10200 sequence.fasta | tail -n 3400 > sequence...
$ head -n 13600 sequence.fasta | tail -n 3400 > sequence...
$ head -n 17000 sequence.fasta | tail -n 3400 > sequence...
$ head -n 20400 sequence.fasta | tail -n 3400 > sequence...
$ head -n 23800 sequence.fasta | tail -n 3400 > sequence...
$ head -n 27200 sequence.fasta | tail -n 3400 > sequence...
$ head -n 30600 sequence.fasta | tail -n 3400 > sequence...
$ tail -n 1942 sequence.fasta > sequence10.fasta
区切りが正しいか元データと手動で照合した
** NCBI BLAST を行う [#f7063dc3]
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 blastout...
blastout5.txt blastout6.txt blastout7.txt blastout8.txt ...
blastout9.txt blastout10.txt > NCBI_BLAST_all.txt
** ゲノムの Sequence ID [#tf19c64d]
NCBI BLAST は検索の対象がゲノム配列だけに絞れないので
(この点が Flybase と大きく違って手間になる)、
ショウジョウバエ(D. melanogaster)ゲノムの Sequence ID を
手掛りにして、結果を絞り込む必要がある
https://www.ncbi.nlm.nih.gov/genome/47?genome_assembly_id...
の 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 から必要なものを抽出 [#d3e74c0f]
perl プログラムを使って、NCBI BLAST の結果から必要なもの...
$ perl extractfromblast.pl NCBI_BLAST_all.txt > NCBI_BLA...
*** perl プログラム:extractfromblast.pl [#o08fa35e]
#!/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: ...
#
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/M...
#
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 の確率が小さい(より確からしい)ものは正しい出力と...
それ以外は1列左にずれて出力するので、
Query の列が染色体名(X、2L、2R、3L、3R、4、Y)になる。~
そこで、Query の列(=行頭)が染色体名(X、2L、2R、3L、3R...
ふたつめにヒットしたものとして削除する。
$ grep -v ^X NCBI_BLAST_exstract.txt | grep -v ^2 | grep...
最後に、余分な改行をエディタを使って、手動で削除する。
なお、BLAST の結果が「ない」ものがあるので、それは個別に...
|Today:&counter(today);|Yesterday:&counter(yesterday);|To...
終了行:
* perl BLASTの出力結果を 加工する 22 January 2021 [#w1abe...
FASTA形式の複数(16,000以上)の塩基配列が
ショウジョウバエゲノム上でどちらの向きか
(プラスかマイナスか)を判定したい
[[NCBI BLAST>https://blast.ncbi.nlm.nih.gov/Blast.cgi]]
とちょっとした perl プログラムを使って、
ざっと判定をしたので、そのメモ
この作業内容では、bioinformatics とは言えませんね (笑)
- [[FlyBase BLAST>http://flybase.org/blast/]]
では、複数の塩基配列を対象とした BLAST が行えない~
- [[DDBJ>http://blast.ddbj.nig.ac.jp/blastn?lang=en]]
では、受け入れ可能ファイルサイズが小さい~
との使いづらさがあったので、NCBI BLAST を使った
** fasta ファイルの加工 [#xfa2bd03]
[[NCBI BLAST>https://blast.ncbi.nlm.nih.gov/Blast.cgi]] の
Nucleotide BLAST では、受け入れるファイルサイズが最大 1,0...
だったので((DDBJ は最大 100,000 とひと桁小さく作業が面倒...
今回は 16,271 個を 1700 個ずつに分割((FASTA形式は2行で1...
$ head -n 3400 sequence.fasta > sequence1.fasta
$ head -n 6800 sequence.fasta | tail -n 3400 > sequence...
$ head -n 10200 sequence.fasta | tail -n 3400 > sequence...
$ head -n 13600 sequence.fasta | tail -n 3400 > sequence...
$ head -n 17000 sequence.fasta | tail -n 3400 > sequence...
$ head -n 20400 sequence.fasta | tail -n 3400 > sequence...
$ head -n 23800 sequence.fasta | tail -n 3400 > sequence...
$ head -n 27200 sequence.fasta | tail -n 3400 > sequence...
$ head -n 30600 sequence.fasta | tail -n 3400 > sequence...
$ tail -n 1942 sequence.fasta > sequence10.fasta
区切りが正しいか元データと手動で照合した
** NCBI BLAST を行う [#f7063dc3]
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 blastout...
blastout5.txt blastout6.txt blastout7.txt blastout8.txt ...
blastout9.txt blastout10.txt > NCBI_BLAST_all.txt
** ゲノムの Sequence ID [#tf19c64d]
NCBI BLAST は検索の対象がゲノム配列だけに絞れないので
(この点が Flybase と大きく違って手間になる)、
ショウジョウバエ(D. melanogaster)ゲノムの Sequence ID を
手掛りにして、結果を絞り込む必要がある
https://www.ncbi.nlm.nih.gov/genome/47?genome_assembly_id...
の 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 から必要なものを抽出 [#d3e74c0f]
perl プログラムを使って、NCBI BLAST の結果から必要なもの...
$ perl extractfromblast.pl NCBI_BLAST_all.txt > NCBI_BLA...
*** perl プログラム:extractfromblast.pl [#o08fa35e]
#!/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: ...
#
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/M...
#
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 の確率が小さい(より確からしい)ものは正しい出力と...
それ以外は1列左にずれて出力するので、
Query の列が染色体名(X、2L、2R、3L、3R、4、Y)になる。~
そこで、Query の列(=行頭)が染色体名(X、2L、2R、3L、3R...
ふたつめにヒットしたものとして削除する。
$ grep -v ^X NCBI_BLAST_exstract.txt | grep -v ^2 | grep...
最後に、余分な改行をエディタを使って、手動で削除する。
なお、BLAST の結果が「ない」ものがあるので、それは個別に...
|Today:&counter(today);|Yesterday:&counter(yesterday);|To...
ページ名: