从作为其他文件中的区域一部分的文件中获取区域(无循环)

Question

我有两个文件:

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中有一条具有相同染色体的行,则start-coverage为> = start-region,end-coverage为<= end-region,然后将其分数保存到新阵列.在所有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或者等等做出替代方案的建议如果有人可以教我如何摆脱我的matlab脚本中的所有循环,我会很高兴.

谢谢

Answer 1

这是一个 Perl 解决方案。我使用哈希（又名字典）通过染色体访问各个范围，从而减少循环迭代的次数。

\n\n

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

\n\n

#!/usr/bin/perl\n\nmy ($rangefile) = @ARGV;\nopen my $rFH, \'<\', $rangefile    or die "Can\'t open $rangefile";\n\n# construct the ranges. The chromosome is used as range key.\nmy %ranges;\nwhile (<$rFH>) {\n    chomp;\n    my @field = split /\\s+/;\n    push @{$ranges{$field[0]}}, [@field[1,2], 0, 0];\n}\nclose $rFH;\n\n# iterate over all the input\nwhile (my $line = <STDIN>) {\n    chomp $line;\n    my ($chrom, $lower, $upper, $value) = split /\\s+/, $line;\n    # only loop over ranges with matching chromosome\n    foreach my $range (@{$ranges{$chrom}}) {\n        if ($$range[0] <= $lower and $upper <= $$range[1]) {\n            $$range[2]++;\n            $$range[3] += $value;\n            last; # break out of foreach early because ranges don\'t overlap\n        }\n    }\n}\n\n# create the report\nforeach my $chrom (sort {$a <=> $b} keys %ranges) {\n    foreach my $range (@{$ranges{$chrom}}) {\n        my $value = $$range[2] ? $$range[3]/$$range[2] : 0;\n        printf "%d %d %d %.1f\\n", $chrom, @$range[0,1], $value;\n    }\n}\n

\n\n

调用示例：

\n\n

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

\n\n

示例输入的输出：

\n\n

1 100 200 6.7\n1 400 600 0.0\n2 600 700 12.5\n

\n\n

（因为（5+7+8）/3 = 6.66\xe2\x80\xa6）

\n