perl のワンライナーに関しての、適当な解説 #1/#2/#3/#4
FASTQ と FASTA について
201501 現在
FASTQ に関しては、リードの配列と qual はそれぞれ一行に表示されています。従って複数行になっている場合は考慮しません。
逆に FASTA に関しては複数行の可能性があるものに関しては、その旨明記の上、コードを書きます。
FASTA => FASTQ(dummy qual)
qual が無いものをダミー qual 40 で。
ワンライナー
改行を考慮せず。
$ perl -Mvars='$n' -lne '$n = ( />(.*?)$/ )[0] and next if $. % 2 == 1 ; print join "\n", "\@$n", $_, "+", q{I} x length' seq.fa > fastq.fq
改行を考慮。
$ perl -F'>' -0777 -lane 'shift @F ; for ( @F ){ my @d = split /\n/ ; printf "@%s\n%s\n+\n%s\n", shift @d, $_ = +(join "", @d), q{I} x length } ' seq.fa > fastq.fq
前者は至って素直に fastq を作ってるだけです。q{I} x length の x を知らなくて面倒な事をやってたのは過去の思い出。
B::Deparse
後者に関して。 B::Deparse 出力をちょっと整形
BEGIN { $/ = undef; $\ = ""; }
LINE:
while (defined($_ = <ARGV>)) {
chomp $_;
our(@F) = split(/>/, $_, 0);
shift @F;
foreach $_ (@F) {
my(@d) = split(/\n/, $_, 0);
printf "\@%s\n%s\n+\n%s\n",
shift @d,
$_ = join('', @d),
'I' x length($_)
;
}
}
-e syntax OK
-
>で split すると配列@Fの最初の要素は何も入ってないのでshift @Fして捨てる - @F の各要素、すなわちリード毎に改行で
splitする - 最初の要素はヘッダ行なので
shift @dしてprintfに渡す - 他の要素は塩基配列なので
joinして$_を上書きすると同時に、printfに渡す -
q{I} x lengthで更新された$_の長さのq{I}を出力
FASTA => FASTQ
こっちは普通に、改行無しの FASTA を想定
下準備
丁寧な人は、 seq.fa と seq.qual のヘッダ数と位置が合ってる事くらいは確認。
$ perl -Mvars='@d' -lne 'next if $. % 2 == 0 ; my $n = ( />(.*?)$/ )[0] ; push @d, $n and next if $ARGV eq q{seq.fa} ; die $_ if $n ne shift @d }{ die if @d != 0' seq.fa seq.qual
現在開いてるファイルのファイル名は、 $ARGV で把握出来る
それさえ分かれば、seq.fa のヘッダをストック。seq.qual のヘッダの時に確認のコードを実行。ってのは難しい話ではない。
丁寧じゃない人でも最初と最後のリードヘッダーくらいは確認しておきたい。
$ diff <( head -1 seq.fa ) <( head -1 seq.qual)
$ diff <( tail -2 seq.fa | head -1) <( tail -2 seq.qual | head -1)
bash, zsh のプロセス置換
ワンライナー
$ perl -Mvars='@d,$n' -lne 's/>// if $. % 2 == 1 ; push @d, $_ and next if $ARGV eq q{seq.fa} ; printf "@%s\n%s\n+\n", splice @d, 0, 2 and next if $. % 2 == 1 ; print' seq.fa seq.qual > fastq.fq
- ファスタヘッダの '>' は邪魔なだけなので
s/>//する。 - seq.fa の時は
@dに追加するだけ。 - seq.qual のヘッダ行の時は、seq.fa の当該リードデータの出力を先にする。当然、ヘッダ行の頭に '@' を追加
- qual のヘッダは
+のみ - qual 行は普通にプリント