How to apply function on bedfile based on a reference bedfile

Question

I have two bed-like files:

File1

chr1     1     10    .    1    Sample1
chr1    20    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    1    Sample3

File2 (reference)

chr1     1     5     .    0 
chr1     5    20     .    1
chr1    50   120     .   -1

I'd like to find a way to apply a function (in my case a subtraction) between the first file and the second, using the second as a reference, and producing a results file containing also the indication of the sample name. For example, in my case the results file would be:

chr1     1      5    .    1    Sample1
chr1     5     10    .    0    Sample1
chr1    20     50    .    1    Sample1
chr1    50    120    .    2    Sample1
chr1   120    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    2    Sample3

I was thinking about implementing it in pandas, but it would be long to do and probably not so efficient. I looked for bedmap but it seems to me that this operation can be carried on (with the sum operation) only within a single file during aggregation, not between two different bed files. So maybe it is not the right tool to use in this case.

Also, it would be important in the final file to retain the sample name.

Answer 1

You could munge your first file into a BED5 file:

$ awk -vOFS="\t" '{ print $1, $2, $3, $6, $5 }' /tmp/f1.bed | sort-bed - > /tmp/f1.bed5

Sort the second file, if necessary:

$ sort-bed /tmp/f2.unknownSort.bed5 > /tmp/f2.bed5

Take the union of the two sorted files:

$ bedops -u /tmp/f1.bed5 /tmp/f2.bed5 > /tmp/union.bed

Here's what the union looks like:

chr1    1       5       .       0       
chr1    1       5       Sample2 1
chr1    1       10      Sample1 1
chr1    5       20      .       1
chr1    20      180     Sample1 1
chr1    50      120     .       -1
chr1    70      100     Sample3 1

We have one name scheme for reference element IDs (.) and a second name scheme for sample IDs (Sample*).

We will exploit this difference in name schemes when we apply operations in a certain order. For this example, we will subtract a reference element score from a sample element score.

We can partition the union, mapping the disjoint elements back to the unioned set:

$ bedops --partition /tmp/union.bed | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed
chr1    1   5   .;Sample2;Sample1   0.000000;1.000000;1.000000
chr1    5   10  Sample1;.   1.000000;1.000000
chr1    10  20  .   1.000000
chr1    20  50  Sample1 1.000000
chr1    50  70  Sample1;.   1.000000;-1.000000
chr1    70  100 Sample1;.;Sample3   1.000000;-1.000000;1.000000
chr1    100 120 Sample1;.   1.000000;-1.000000
chr1    120 180 Sample1 1.000000

Each of the fourth and fifth columns are semi-colon-delimited strings. Each element is ordered by the order in which it is discovered in a bedmap operation.

Both columns have the same ordering or same one-to-one correspondence of ID to signal score. This is very useful, because we can turn these into a pair of ordered arrays and apply some logic, depending on the answer we want back.

For example, here is some awk logic to parse those two strings and subtract a reference element score from a sample element score:

$ bedops --partition /tmp/union.bed \
    | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed \
    | awk -vFS="\t" -vOFS="\t" \
        '{ \
             i=split($4,a,";"); \
             j=split($5,b,";"); \
             if ((i==1) && (a[i]==".")) { pass; } \
             else { \
                 for (k=1; k<=i; k++) { \
                     if (a[k]==".") { refIdx=k; } \
                 } \
                 for (k=1; k<=i; k++) { \
                     if (k!=refIdx) { print $0, ".", b[k]-b[refIdx], a[k]; } \
                 } \
             } \
         }' \
    | cut -f1-3,6-8 \
    | sort-bed - \
    > answer.bed

The file answer.bed will look like this:

chr1    1   5   .   1   Sample1
chr1    1   5   .   1   Sample2
chr1    5   10  .   0   Sample1
chr1    20  50  .   1   Sample1
chr1    50  70  .   2   Sample1
chr1    70  100 .   2   Sample1
chr1    70  100 .   2   Sample3
chr1    100 120 .   2   Sample1
chr1    120 180 .   1   Sample1

This partitioning of intervals is different from the one in your question. Perhaps this is due to how you partitioned your set. In any case, perhaps this would give you an idea of how to apply custom operations over a disjoint partitioning.

The main idea is that this part of the awk script finds the score associated with the reference file ("File2"):

for (k=1; k<=i; k++) {
    if (a[k]==".") { refIdx=k; }
}

Any index other than refIdx corresponds to a Sample* element. That means we can do operations between Sample* scores and the specific reference score.

We also leave out mappings where a reference score does not associate with any Sample* element:

if ((i==1) && (a[i]==".")) { pass; }

Thus, all the operations are always between a pairing of a sample and the mapped reference element.

The cut and sort-bed statements at the end just clean up the result so that you get something that looks like what you are expecting at the end, which is also sorted correctly and ready to use for downstream set operations.