gpt4 book ai didi

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

转载 作者:太空宇宙 更新时间:2023-11-03 19:35:02 25 4
gpt4 key购买 nike

我有两个文件:

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 亿行。

我基本上想获得 regions.txt 中每个区域中 coverage.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 脚本中的所有循环,我也会很高兴。

谢谢

最佳答案

这是一个 Perl 解决方案。我使用散列(又名字典)通过染色体访问各种范围,从而减少循环迭代的次数。

这可能很有效,因为我不会在每个输入行上对 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, 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 - 从文件中获取属于其他文件中区域的区域(无循环),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/12872469/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com