![[PukiWiki]](image/tomarufront.png "[PukiWiki]")
#author("2021-01-22T14:59:08+09:00","","") #author("2021-01-22T15:00:41+09:00","","") * perl BLASTの出力結果を 加工する 22 January 2021 [#w1abec6f] 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,000,000 だったので((DDBJ は最大 100,000 とひと桁小さく作業が面倒なことがわかり、NCBI を使うことにした。DDBJ なら結果をメールで受け取れるのが便利。))、 今回は 16,271 個を 1700 個ずつに分割((FASTA形式は２行で１組)) $ 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 を行う [#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 blastout4.txt ［改行せず］ 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=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 プログラム [#d5e9e5a4] 今回は、query配列がゲノム上でプラス鎖かマイナス鎖のどちらかを 知りたいだけなので、それだけを取り出し、タブ区切りテキストを 標準出力に吐き出すことにした また、念のために、染色体（腕）とゲノムの Sequence ID も同じ行に出力させる なお、print の部分を push などに変えれば、更なる加工ができるはず *** NCBI BLAST から必要なものを抽出 [#d3e74c0f] perl プログラムを使って、NCBI BLAST の結果から必要なものを抽出 $ perl extractfromblast.pl NCBI_BLAST_all.txt > NCBI_BLAST_exstract.txt *** 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: 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 の結果、複数のゲノム染色体がヒットする場合は…… [#u8754d3d] このプログラムは BLAST の結果の出力ファイルの上から処理していくので、 確率が大きい（余り確からしくない）ものには、Query名が付与されない。 そこで、そういったデータを後から削る必要がある。~ BLAST の確率が小さい（より確からしい）ものは正しい出力となり、 それ以外は１列左にずれて出力するので、 Query の列が染色体名（X、2L、2R、3L、3R、4、Y）になる。~ そこで、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:&counter(today);|Yesterday:&counter(yesterday);|Total:&counter(); since 22 January 2021|
