-1

私はバイオインフォマティクスの初心者で、ペアエンドの MiSeq データ (現在は 1 つの fastq ファイル) を 2 つのファイルに分割する小さな Bio Perl コードに取り組んでおり、各ファイルにはペアの一方の端が含まれています。ペアエンドリードの異なるエンドは、fastq ヘッダーのスペースの後の1または2で区別できます。このファイルは、コマンド ラインで「head」を使用する例のように、典型的な fastq 形式に従います。

@M00763:6:000000000-A1U80:1:1101:12620:1732 1:N:0:1
TTATACTC
+
@A@AA@A@
@M00763:6:000000000-A1U80:1:1101:12620:1732 2:N:0:1
T
+
E

一致を使用してヘッダーの 1 または 2 をターゲットにしようとするコードを作成しました。私は Bio::SeqIO を使用していますが、perl は fastq 形式を認識していないようで、このエラーが発生し続けます:

MSG: Could not guess format from file/fh
STACK: Error::throw
STACK: Bio::Root::Root::throw /sw/lib/perl5/5.12.3/Bio/Root/Root.pm:472
STACK: Bio::SeqIO::new /sw/lib/perl5/5.12.3/Bio/SeqIO.pm:389
STACK: SplitPairedEndReads.pl:7

誰かが私のエラーを見つけて修正するのを手伝ってくれますか? BioPerl Web サイトから入手できる情報は、Bio::SeqIO が fastq 形式を認識できる必要があることを示しています。

ここに私が書いたコードがあります：

#!/usr/bin/perl 

use Bio::SeqIO;
use Bio::SeqIO::fastq;


$seqout1 = Bio::SeqIO->new(-file => ">peread1.fastq" -format => "fastq",);
$seqout2 = Bio::SeqIO->new(-file => ">peread2.fastq" -format => "fastq",);

$seqio_obj = Bio::SeqIO->new(-file => "AIS351_Strin1edit.fastq", -format => "fastq",
                         -alphabet => "dna" );
$seq_obj = $seqio_obj->next_seq;

while ($seq_obj = $seqio_obj->next_seq) { 
    $name = $seq_obj->desc; if($name=~ / 1:/) {$seqout1->write_seq($seq_obj);
     } else { $seqout2->write_seq($seq_obj); 

    }
}

私の初心者の知識に助けてくれてありがとう。

〜アル

質問の更新:

行のコンマ エラーを修正しましたnewが、コードを実行すると次のエラーが発生します。

------------- EXCEPTION: Bio::Root::Exception -------------
MSG: No description line parsed
STACK: Error::throw
STACK: Bio::Root::Root::throw /sw/lib/perl5/5.12.3/Bio/Root/Root.pm:472
STACK: Bio::SeqIO::fastq::next_dataset /sw/lib/perl5/5.12.3/Bio/SeqIO/fastq.pm:71
STACK: Bio::SeqIO::fastq::next_seq /sw/lib/perl5/5.12.3/Bio/SeqIO/fastq.pm:29
STACK: samplesettrim.pl:10
-----------------------------------------------------------

私が行ったすべての読み取りは、BioPerl 自体の FASTQ パーサーにいくつかの問題があることを示しているようです。私は初心者であり、プログラミングのスキルを向上させようとしているので (私は完全に独学です)、このコードを機能させることを望んでいました。これは遅く、おそらく大きな FASTQ ファイルを操作するための最良の方法ではないというコメントに同意します。

+ 記述子に関しては、私のファイルを他のソフトウェア プログラム (例: CLC) で使用できるようにするために必要ですか、それとも FASTQ でその行を削除することで問題を解決できますか? + には、読み取りに関する品質情報は実際には含まれていませんね。

入力していただきありがとうございます。

4

3 に答える 3

2

への呼び出しでは、すべてのリスト項目の間にコンマを追加する必要がありますnew。変化する：

$seqout1 = Bio::SeqIO->new(-file => ">peread1.fastq" -format => "fastq",);
$seqout2 = Bio::SeqIO->new(-file => ">peread2.fastq" -format => "fastq",);

に：

$seqout1 = Bio::SeqIO->new(-file => ">peread1.fastq", -format => "fastq",);
$seqout2 = Bio::SeqIO->new(-file => ">peread2.fastq", -format => "fastq",);
于 2013-12-04T17:21:31.150 に答える
1

BioPerl は非常に遅いため、Fastq データには使用しないことをお勧めします (以下の私のコメントを参照してください)。このタスクにはPairfqを使用できます。これは、この目的のために設計されたものの 1 つであるためです (完全な開示: 私は作成者です)。仕組みは次のとおりです。

pairfq splitpairs -i AIS351_Strin1edit.fastq -f AIS351_Strin1edit_1.fastq -r AIS351_Strin1edit_2.fastq

私のベンチマークでは、これは BioPerl で同等のタスクを実行するよりも約 300 倍高速です。たとえば、Bio::SeqIO で 100 万件の Fastq レコードを読み取るのに 465 秒かかることを測定しましたが、上記のコードでは約 1.5 秒で完了します。5 億件のレコードがある場合、64 時間と 11 分の差になります。そのため、NGS データに BioPerl を使用することを強くお勧めしません。私は BioPerl を毎日使っているのでバッシングしているわけではありませんが、この問題には注意してください。

あなたのコメントのエラーについて、BioPerl パーサーはあなたの「+」行にあるものが好きではありません。「+」の後に何もないか、シーケンス ヘッダーと一致する必要があります。実際のデータを見ずに具体的に言うのは難しいですが、行末の問題やその他の問題である可能性もあります。

編集:すべてのスクリプトの先頭にuse strict;andを配置する必要があります。use warnings;また、何かをしようとする前に、ファイルが存在するかどうかをテストすることをお勧めします (BioPerl でファイルを読み取ろうとするなど)。最後の質問については、 FASTQ形式についてお読みになることをお勧めします。レコードから行を削除することはできません。そうしないと、有効な FASTQ になりません。マイナーな点は、適切なクラスのロードを処理するuse Bio::SeqIO::fastq;ため、必要がないことです。Bio::SeqIO

あなたが投稿したものは実際のデータのようには見えないため、問題の原因を特定するのは簡単ではありません.

于 2013-12-04T23:23:33.103 に答える
0

このスニペットを使用すると、目的を達成できます。

#!/usr/bin/perl
use warnings;
use strict; 

my @array = ('@M00763:6:000000000-A1U80:1:1101:12620:1732 1:N:0:1
TTATACTC
+
@A@AA@A@',
'@M00763:6:000000000-A1U80:1:1101:12620:1732 2:N:0:1
T
+
E');

foreach (@array){
        if (/\s+1:/) {
            print "1st pair: $_\n"; # You could redirect this to first.OUTFILE
         }
        if (/\s+2:/) {
            print "2nd pair: $_\n"; # You could redirect this to second.OUTFILE
         }

}

どちらが印刷されますか:

1st pair: @M00763:6:000000000-A1U80:1:1101:12620:1732 1:N:0:1
TTATACTC
+
@A@AA@A@
2nd pair: @M00763:6:000000000-A1U80:1:1101:12620:1732 2:N:0:1
T
+
于 2013-12-04T17:17:11.703 に答える