SINGLE FASTA

perl のワンライナーに関しての、適当な解説 #1/#2/#3/#4

Genbank format から sequence を取り出す

Genbank format ってなんぞ？って人は、こちら

普通に FASTA 選択しなおしてダウンロードすれば済むが、他人様がお持ちの場合は、バージョン違いとか、余計な地雷を踏む事にもなりかねないので、

$ perl -lne 'next if !( /ORIGIN/ .. m{//} ); next if !/\d/ ; s/[\d\s]//g ; printf qq{>%s\n%s\n}, $ARGV, $_ ' GENBANK.txt

• unless は地雷
• $ARGV が気に入らないなら ACCESSION から文字列取ってくるとか、お好きに

アレな作業をしなきゃならん時には、

$ bl2seq -i <( 上のコード ) -j SOMETHING.fa -p blastn

consensus の N の位置と長さを知る

ゲノムプロジェクト毎に公表はしているものですが、、、

$ perl -Mvars='$s' -lne 'next if />/; $s .= $_ }{ $, = "\t";  print $-[0] + 1, length $1 while $s =~ /(?<!N)(N{1,})(?!N)/ig' consensus.fa

大きな N の連続だけを知りたい場合 (N{1,}) を (N{1000,}) などと変更する。

本当に正しく算出出来てるのかどうか不安になったなら、

$ perl -Mvars='$s' -lne ' next if />/ ; $s .= $_ }{ print substr $s, INDEX, LENGTH' consensus.fa
# INDEX は、posision - 1 
# LENGTH は、長さ

などとして、確認すれば良い。

コードゴルフ的には、 $, = "\t" だが、 print join "\t" の方が普通かな、、、

異物混入

もちろん、これは承知の上で。

while $s =~ /(?<!N)(N{1,})(?!N)/ig を　while $s =~ /(?<![^ACGTN])([^ACGTN]{1,})(?![^ACGTN])/ig
とすれば、異物が入ってた場合の、場所と長さを見る事も出来る。

制限酵素断片

例えば、HindIII （ 'A/AGCTT' で切れる。）での断片がどの様に表われるか知りたい場合。

配列そのまま

$ perl -Mvars='$s' -lne 'next if />/; $s .= $_ }{ my @d = split /(?<=A)(?=AGCTT)/, $s ; print for @d' consensus.fa

長さ

$ perl -Mvars='$s' -lne 'next if />/; $s .= $_ }{ my @d = split /(?<=A)(?=AGCTT)/, $s ; print length for @d' consensus.fa

うん？ BssHII（ 'G/CGCGC' で切れる。） とかって dry では、正確にはどう切れるんだっけか？

GC の長い繰替えし配列なんて結構あると思うが、、、まあ、 wet だと誤差範囲の一言で済む話だけど。

次回へ。