Answers to: fasta2fastq

Answer by yasugoyasu

yasugoyasu — Sun, 05 Jun 2011 22:06:34 +0900

gaouさま、nob_fjさま

返信ありがとうござました。 maqに入っているかもしれないと思い、maqも見たのですが、fq_all2std.plは見落としていました。教えていただきありがとうございました。

説明不足でしたが、exomeキャプチャーしてきたsequenceをde novo assemblyしてから、referenceにmappingしたものと、単純にsolexaのreadをmappingしたものとでの効率（% on-targetや% read coverage）の比較をしたいと考えていました。 referenceがない種からどの程度のコーディング情報が取れるかの予備実験です。

assemblyしたcontigは24bp~4000bp程度（exome_assembly）で、(i) bwa+bwa-swと(ii) blatの両方を使ったmappingをやろうと思っています。 blatは最初のmappingまではやり、pslCDnaFilterによるpslのfilteringを行なうと思っていますが、pslCDnaFilterのインストールに手間取っています。 bwaのmappingは200bpまでとmanualに書かれていますが、bwa-swは最大100kbくらいまではmappingできると書いてあります。 bwaは通常のmappingで使っていますが、bwa-swはいまだに使った事がありません。

bwa+bwa-swとblatによるmappingのどちらがいいか、もしご存知の方がおられたら教えていただけますでしょうか。

よろしくお願いします。

Answer by nob_fj

nob_fj — Sun, 05 Jun 2011 19:03:55 +0900

NGS surfer's Wiki(コンバータ・パーサ) にも記載していますが、純粋にFASTAファイルをFASTQに変換するには？という質問に対する質問への回答としては、

MAQというツールをインストールすると付属しているperlスクリプトを使用する。

$ head -n 4 input.fa
>seq1
ATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATG
>seq2
ATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATG
$ /toolPath/maq-0.7.1/scripts/fq_all2std.pl fa2std -q 40 input.fa > output.fastq
$ head -n 8 output.fastq
@seq1
ATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATG
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@seq2
ATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATG
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

また、

自前でスクリプト書く。

多くのUNIX/LINUX環境では、sed、awkは入っていると思います。cygwinでもたぶん動きます。充分テストをしていないので、お勧めはしません。ツール(G-language,BioPerl,maq)などのインストールがどうしても面倒だという場合のみ使用を検討下さい。また、配列行が複数行になることは想定していません。

$ sed -n "/^>/N;s/^>//;s/\n/\t/p;" input.fa | awk '{qual=$2;gsub(/./,"I",qual);print "@"$1"\n"$2"\n+"$1"\n"qual}' > output2.fastq
$ head -n 8 output2.fastq
@seq1
ATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATG
+seq1
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@seq2
ATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATG
+seq2
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

BWAが今回のマッピングに適しているかというgaouさんの指摘は、もっともな気がしますので、ご留意ください。

Answer by gaou

gaou — Sat, 04 Jun 2011 03:00:19 +0900

G-language GAEなら

use G;
my %fasta = readFile('file.fasta', -format=>"fasta");
for my $key (keys %fasta){
    print to_fastq($fasta{$key}, -name=>$key);
}

BioPerlなら http://www.bioperl.org/wiki/Merging_separate_sequence_and_quality_files_to_FASTQ のようにライブラリ毎にいろいろ書き方はありますが、そもそもBWAのアルゴリズムは長い配列には向かないのでアセンブルする前にshort readのままmappingした方が良いですし、アセンブルしたものがTranscriptならExonerate、ゲノムの一部ならMUMERなどを使った方が良いかと思います。