perl – 从文件中获取属于其他文件中的区域的区域(无循环)

regions.txt：第一列是染色体名称,第二列是开始和结束位置.

1  100  200
1  400  600
2  600  700

coverage.txt：第一列是染色体的名字,再第二和第三的开始和结束位置,而最后一列是得分.

1 100 101  5
1 101 102  7 
1 103 105  8
2 600 601  10
2 601 602  15

这个文件非常庞大,大约15GB,大约有3亿行.

我基本上想得到coverage.txt中所有得分的均值,这些得分位于regions.txt中的每个区域.

换句话说,开始在regions.txt第一行中,如果在一个coverage.txt线,其具有在相同的染色体上,启动覆盖是&GT =启动区,和最终的覆盖为< =最终区域,然后将其得分保存到新阵列.经过所有coverages.txt完成搜索打印区域染色体,开始,结束,并已发现所有得分的平均值. 预期产量：

1  100 200 14.6   which is (5+7+8)/3
1  400 600 0      no match at coverages.txt
2  600 700 12.5   which is (10+15)/2

我内置了以下MATLAB脚本,因为我不得不遍历coverage.txt很多时间,这需要很长的时间.我不知道如何制作一个快速awk类似的脚本.

我的matlab脚本

fc = fopen('coverage.txt','r');
ft = fopen('regions.txt','r');
fw = fopen('out.txt','w');

while feof(ft) == 0

linet = fgetl(ft);
scant = textscan(linet,'%d%d%d');
tchr = scant{1};
tx = scant{2};
ty = scant{3};
coverages = [];

    frewind(fc);
    while feof(fc) == 0

    linec = fgetl(fc);
    scanc = textscan(linec,'%d%d%d%d');
    cchr = scanc{1};
    cx = scanc{2};
    cy = scanc{3};
    cov = scanc{4};


        if (cchr == tchr) && (cx >= tx) && (cy <= ty)

            coverages = cat(2,coverages,cov);

        end

    end

    covmed = median(coverages);
    fprintf(fw,'%d\t%d\t%d\t%d\n',tchr,tx,ty,covmed);

end

任何建议使用AWK,Perl中,使替代或…等我将ASLO很高兴,如果有人可以教我如何摆脱所有循环的在我的MATLAB程序.

谢谢

解决方法

这是一个Perl解决方案.我使用哈希(aka dictionaries)通过染色体访问各种范围,从而减少循环迭代次数.

这可能是有效的,因为我没有对每个输入行上的regions.txt进行完整循环.使用多线程时,效率可能会进一步提高.

#!/usr/bin/perl

my ($rangefile) = @ARGV;
open my $rFH,'<',$rangefile    or die "Can't open $rangefile";

# construct the ranges. The chromosome is used as range key.
my %ranges;
while (<$rFH>) {
    chomp;
    my @field = split /\s+/;
    push @{$ranges{$field[0]}},[@field[1,2],0];
}
close $rFH;

# iterate over all the input
while (my $line = <STDIN>) {
    chomp $line;
    my ($chrom,$lower,$upper,$value) = split /\s+/,$line;
    # only loop over ranges with matching chromosome
    foreach my $range (@{$ranges{$chrom}}) {
        if ($$range[0] <= $lower and $upper <= $$range[1]) {
            $$range[2]++;
            $$range[3] += $value;
            last; # break out of foreach early because ranges don't overlap
        }
    }
}

# create the report
foreach my $chrom (sort {$a <=> $b} keys %ranges) {
    foreach my $range (@{$ranges{$chrom}}) {
        my $value = $$range[2] ? $$range[3]/$$range[2] : 0;
        printf "%d %d %d %.1f\n",$chrom,@$range[0,1],$value;
    }
}

示例调用：

$perl script.pl regions.txt <coverage.txt >output.txt

示例输入的输出：

1 100 200 6.7
1 400 600 0.0
2 600 700 12.5

(因为(5 7 8)/ 3 = 6.66 …)

perl – 从文件中获取属于其他文件中的区域的区域(无循环)

解决方法

相关文章