-1

我是生物信息学的初学者，我一直在编写一些 Bio Perl 代码来将配对的末端 MiSeq 数据（当前位于 1 个 fastq 文件中）拆分为 2 个文件，每个文件包含该对的一个末端。配对末端 reads 的不同末端可以通过fastq 标头中空格后的1或2来区分。该文件遵循典型的 fastq 格式，例如在命令行中使用“head”：

@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 网站提供的信息表明 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用于此任务，因为这是它的设计目的之一（完全披露：我是作者）。以下是它的工作原理：

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 分钟的差异。这就是为什么强烈反对将 BioPerl 用于 NGS 数据的原因。我不是因为我每天都在使用 BioPerl 而抨击它，但请注意这个问题。

关于您评论中的错误，BioPerl 解析器不喜欢您的“+”行中的内容。'+' 之后必须没有任何内容，否则它必须与序列头匹配。没有看到真实数据很难说具体，也可能是行尾问题或其他问题。

编辑：您需要将use strict;anduse 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 回答