我有多个制表符分隔的 fastq 文件。我想匹配每次读取的第二行,如果匹配则将其旁边的值相加。例如:
file1.fq
>1
ATGCCGTT file1:1
+
HHHHKKKK
file2.fq
>2
ATGCCGTT file2:3
+
JJKHHTTT
>3
ATTCCAAC file2:1
+
=#GJLMNB
我想要的输出是这样的:
output.txt
ATGCCGTT file1:1 file2:3 count:4
ATTCCAAC file2:1 count:1
我写的代码是:
#!/usr/bin/env perl
use strict;
use warnings;
no warnings qw( numeric );
my %seen;
$/ = "";
while () {
chomp;
my ($key, $value) = split ('\t', $_);
my @lines = split /\n/, $key;
my $key1 = $lines[1];
$seen{$key1} //= [ $key ];
push (@{$seen{$key1}}, $value);
}
foreach my $key1 ( sort keys %seen ) {
my $tot = 0;
my $file_count = @ARGV;
for my $val ( @{$seen{$key1}} ) {
$tot += ( split /:/, $val )[0];
}
if ( @{ $seen{$key1} } >= $file_count) {
print join( "\t", @{$seen{$key1}});
print "\tcount:". $tot."\n\n";
}
}
该代码适用于小文件,但当我想比较大文件时,它会占用整个内存,导致脚本运行而没有结果。我想修改脚本,使其不占用内存。我不想使用任何模块。我认为如果我一次只加载一个文件到内存中,它会节省内存,但无法做到这一点。请帮助修改我的脚本。
答案1
你有没有尝试过awk
?不确定它能更好地处理大文件perl
,但可能值得一试:
在你的 awk 脚本中:
BEGIN {
RS=">[0-9]+"
}
FNR==1{next}
NR==FNR {
a[$1]++
next
}
$1 in a {
b[$1]++
next
}
{
c[$1]++
}
END {
for (key in a) {
if (b[key] == "") {
printf key"\tfile1:"a[key]"\t\tcount:"a[key]"\n"
} else {
printf key"\tfile1:"a[key]"\tfile2:"b[key]"\tcount:"a[key]+b[key]"\n"
}
}
for (key in c) {
printf key"\t\tfile2:"c[key]"\tcount:"c[key]"\n"
}
}
运行它:
$ awk -f myscript.awk file1 file2 > output.txt
测试它:
文件1
>1
ATGCCGTT file1:1
+
HHHHKKKK
>2
ATTCCAACg file2:1
+
=#GJLMNB
文件2
>2
ATGCCGTT file2:3
+
JJKHHTTT
>3
ATTCCAAC file2:1
+
=#GJLMNB
终端输出:
ATTCCAACg file1:1 count:1
ATGCCGTT file1:1 file2:1 count:2
ATTCCAAC file2:1 count:1
答案2
将这些神秘咒语添加到您的程序中
use DB_File;
my %seen;
unlink '/tmp/translation.db';
sleep 2;
tie ( %seen, 'DB_File', '/tmp/translation.db' )
or die "Can't open /tmp/translation.db\n";
并且您的哈希将不再驻留在内存中,而是驻留在磁盘上的数据库中。您可以将其余代码原样保留。确实,我使用了 DB_File 模块,但确实没有理由不这样做。它伴随着每一个珀尔开箱即用安装,因此您无需安装它或任何东西。
如果我的哈希值变得非常大,我会一直使用这种方法,并且我发现,在通过一些模糊定义的休斯点后,事情会加快很多。