ログイン 概要 よくある質問

初めて質問します。 solexaのreadをde novo assembleして出来たcontigをbwaもしくはbwa-swを用いてmappingしたいと考えています。 bwaにinputするためには、contigの.fastaファイルを.fastqに変換する必要があると思いますが、そのような変換ツールを知っている方がおられたら教えていただけないでしょうか。調べたのですが調べ方が悪いのか見つかりませんでした。

具体例

【変換前(.fasta)】

>SEQ_ID

GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT

【変換後(.fastq)】

@SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

qualityはphred_score(Phred+33)=40(I)で固定で構いません。

ご存知の方がいたら、教えていただけると助かります。 よろしくお願いします。

質問日 Jun 04 '11 at 02:35

yasugoyasu's gravatar image

yasugoyasu
2623


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などを使った方が良いかと思います。

回答日 Jun 04 '11 at 03:00

gaou's gravatar image

gaou ♦♦
22125

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さんの指摘は、もっともな気がしますので、ご留意ください。

回答日 Jun 05 '11 at 19:03

nob_fj's gravatar image

nob_fj ♦
50761328

edited Jun 05 '11 at 19:19

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のどちらがいいか、もしご存知の方がおられたら教えていただけますでしょうか。

よろしくお願いします。

回答日 Jun 05 '11 at 22:06

yasugoyasu's gravatar image

yasugoyasu
2623

あなたの回答
プレビューをトグルする

この質問をフォローする

By Email:

Once you sign in you will be able to subscribe for any updates here

By RSS:

回答

回答とコメント

タグ:

×47
×4
×4

質問日: Jun 04 '11 at 02:35

閲覧数: 5,081 回

最終更新日: Jun 05 '11 at 22:06

powered by OSQA