6

私は2つのファイルを持っています：

region.txt：最初の列は染色体名、2番目と3番目は開始位置と終了位置です。

1  100  200
1  400  600
2  600  700

Coverage.txt：最初の列は染色体名、2番目と3番目は開始位置と終了位置、最後の列はスコアです。

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スクリプトのすべてのループを取り除く方法を誰かが教えてくれれば嬉しいです。

ありがとう

4

4 に答える 4

4

これが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…）

于 2012-10-13T12:43:15.543 に答える
1

通常、ファイルをRにロードして計算しますが、そのうちの1つが非常に大きいことを考えると、これは問題になります。ここにあなたがそれを解決するのを助けるかもしれないいくつかの考えがあります。

  1. coverage.txt染色体による分割を検討してください。これにより、計算の負担が軽減されます。

  2. ループする代わりにcoverage.txt、最初にregions.txt完全なものをメモリに読み込みます（私はそれがはるかに小さいと思います）。地域ごとに、スコアと数値を保持します。

  3. coverage.txt行ごとに処理します。各行について、この特定のストレッチが属する染色体と領域を決定します。これにはある程度のフットワークが必要になりますが、regions.txt大きすぎない場合は、より効率的である可能性があります。地域のスコアにスコアを追加し、番号を1つ増やします。

別の最も効率的な方法では、両方のファイルを最初に染色体で、次に位置でソートする必要があります。

  1. から行を取りregions.txtます。染色体と位置を記録します。前のループから残っている行がある場合は、3に進みます。; それ以外の場合は2に進みます。

  2. から行を取りcoverage.txtます。

  3. 現在のリージョン内にあるかどうかを確認してください。

    • はい：スコアをリージョンに追加し、番号を増やします。2に移動します。
    • いいえ：スコアを数値で割り、現在の領域を出力に書き込み、1に進みます。

この最後の方法では、微調整が必​​要ですが、最も効率的です。各ファイルを1回だけ処理する必要があり、メモリにほとんど何も保存する必要はありません。

于 2012-10-13T11:47:42.003 に答える
1

とを使用する1つの方法がjoinありawkます。次のように実行します：

join regions.txt coverage.txt | awk -f script.awk - regions.txt

内容script.awk：

FNR==NR && $4>=$2 && $5<=$3 { 
    sum[$1 FS $2 FS $3]+=$6
    cnt[$1 FS $2 FS $3]++
    next
}

{
    if ($1 FS $2 FS $3 in sum) {
        printf "%s  %.1f\n", $0, sum[$1 FS $2 FS $3]/cnt[$1 FS $2 FS $3]
    }
    else if (NF == 3) {
        print $0 "  0"
    }
}

結果：

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

または、これがワンライナーです。

join regions.txt coverage.txt | awk 'FNR==NR && $4>=$2 && $5<=$3 { sum[$1 FS $2 FS $3]+=$6; cnt[$1 FS $2 FS $3]++; next } { if ($1 FS $2 FS $3 in sum) printf "%s  %.1f\n", $0, sum[$1 FS $2 FS $3]/cnt[$1 FS $2 FS $3]; else if (NF == 3) print $0 "  0" }' - regions.txt
于 2012-10-13T14:25:53.593 に答える
0

カバレッジをリージョンにビニングする簡単なMATLABの方法は次のとおりです。

% extract the regions extents
bins = regions(:,2:3)';
bins = bins(:);

% extract the coverage - only the start is needed
covs = coverage(:,2);

% use histc to place the coverage start into proper regions
% this line counts how many coverages there are in a region
% and assigns them proper region ids.
[h, i]= histc(covs(:), bins(:));

% sum the scores into correct regions (second output of histc gives this)
total = accumarray(i, coverage(:,4), [numel(bins),1]);

% average the score in regions (first output of histc is useful)
avg = total./h;

% remove every second entry - our regions are defined by start/end
avg = avg(1:2:end);

これは、リージョンが重複していないことを前提として機能しますが、そうであると思います。また、coverageファイル内のすべてのエントリは、特定の領域に分類される必要があります。

また、ファイル全体の読み取りを避けたい場合は、カバレッジに対してこのアプローチを「ブロック」するのは簡単です。必要なのはbins、おそらく小さいリージョンファイルだけです。カバレッジをブロック単位で処理し、段階的に追加しtotalて、最終的に平均を計算できます。

于 2012-10-13T12:07:17.483 に答える