Merging regions according to their identifier

Question

I have a huge file (20 GB) which has a range of genomic locations, and for each location there is an identifier(4th column), which is sometimes the same.

file1.txt

chr1  10  20  ABC
chr1  13  20  ABC
chr1  14  21  ABC
chr1  22  27  ABC
chr1  29  37  ABC
chr2  15  21  JJJ
chr2  21  31  JJJ
chr2  23  27  JJJ
chr2  35  56  JJJ
chr2  25  26  MMM
chr3  32  42  MMM
chr3  45  76  MMM
chr3  88  101 MMM
chr3  101 105 MMM

I am trying to merge the ranges of column 2 and 3, if they are in the same chromosome and having the same identifier ( having the 1st and 4th columns the same).

I have tried to use bedtools, as shown in this post, which looks very similar but not the same with mine, as they are looking for different outputs ( range vs merge )

So the answers in that post was with using groupby in bedtools2,source code from github.

I have applied the same as the following:

./groupBy -i ../file1.txt -g 1,4 -c 2,3 -o min,max | awk -v OFS='\t' '{print $1, $3, $4, $2}'

chr1    10      37      ABC
chr2    15      56      JJJ
chr2    25      26      MMM
chr3    32      105     MMM

But it only connects the ranges of identifiers.. It groups by, and does not merge.

And, if I use bedtools v2.26 merge function :

>sort -k1,1 -k2,2n file1.txt > file2.txt
>cat file2.txt
chr1    10      20      ABC
chr1    13      20      ABC
chr1    14      21      ABC
chr1    22      27      ABC
chr1    29      37      ABC
chr2    15      21      JJJ
chr2    21      31      JJJ
chr2    23      27      JJJ
chr2    25      26      MMM
chr2    35      56      JJJ
chr3    32      42      MMM
chr3    45      76      MMM
chr3    88      101     MMM
chr3    101     105     MMM


bedtools merge -i file2.txt 
chr1    10      21
chr1    22      27
chr1    29      37
chr2    15      31
chr2    35      56
chr3    32      42
chr3    45      76
chr3    88      105

Which, merges without seeing the identifiers as it also gives me this output.

What I want is to merge the ranges if they overlap, if having the same 1st and 4th column as:

expected_output.txt


chr1  10  21  ABC
chr1  22  27  ABC
chr1  29  37  ABC
chr2  15  31  JJJ
chr2  35  56  JJJ
chr2  25  26  MMM
chr3  32  42  MMM
chr3  45  76  MMM
chr3  88  105 MMM

I only have access to 15GB of RAM, so I can't store the entire file in memory.

Answer 1

The following python script will do what you want and should be relatively memory efficient. It processes a single chromosome at a time, so either sort the input or ensure that at least entries in a group are next to each other.

#!/usr/bin/env python
import sys

# Process a sorted list of tuples (e.g., [(0, 10), (10, 20), ...])
# Overlapping items will be merged and a list of the merged regions is
# returned (e.g., [(0, 20), ...])
def merge(values):
    start, end = values[0]
    o = []
    for v in values[1:]:
        if v[0] <= end:
            if v[1] > end:
                end = v[1]
        else:
            o.append((start, end))
            start, end = v
    o.append((start, end))
    return o


# Process a single chromosome, one group at a time
def processChrom(d, chrom):
    for k, v in d.items():
        v.sort()
        tuples = merge(v)
        for t in tuples:
            print("{}\t{}\t{}\t{}".format(chrom, t[0], t[1], k))


def main():
    currentChrom = None
    d = dict()
    for line in sys.stdin:
        cols = line.split()
        if cols[0] == currentChrom:
            if cols[3] not in d:
                d[cols[3]] = []
            d[cols[3]].append((int(cols[1]), int(cols[2])))
        else:
            if currentChrom:
                processChrom(d, currentChrom)
            d = dict()
            d[cols[3]] = []
            d[cols[3]].append((int(cols[1]), int(cols[2])))
            currentChrom = cols[0]

    if currentChrom:
        processChrom(d, currentChrom)

if __name__ == "__main__":
    main()

Usage is python foo.py < foo.txt and the output is:

chr1    10  21  ABC
chr1    22  27  ABC
chr1    29  37  ABC
chr2    25  26  MMM
chr2    15  31  JJJ
chr2    35  56  JJJ
chr3    32  42  MMM
chr3    45  76  MMM
chr3    88  105 MMM

You'll notice that this isn't necessarily sorted, so you may want to pass the output through bedtools sort or sort.

Answer 2

If I understand correctly, you could create dummy chromosomes made by merging chromosome and identifier, merge with bedtools, split back chromosomes and identifiers. E.g.

awk -v OFS='\t' '{print $1"_"$4, $2, $3}' file1.txt \
| mergeBed \
| sed 's/_/\t/' \
| awk -v OFS='\t' '{print $1, $3, $4, $2}'

(Assuming "_" is not part of chromosome or identifier name, adjust accordingly if so).

Answer 3

Just to add this can be done with bedtools by processing each identifier separately which may or may not be practical given your large file.

for ID in $(perl -ane 'print $F[3] . "\n";' file1.txt | sort | uniq);do 
    grep $ID file1.txt | \
    perl -ane 'chomp;print join("\t",@F[0..$#F]) . "\n";' | \
    sort -V | \
    bedtools merge -i stdin -c 4 -o distinct
done | sort -V

Answer 4

Here is a fairly simple (and hopefully readable) native Python solution. It assumes that the bed file is sorted prior to parsing:

from itertools import groupby

class BedRecord(object):
    def __init__(self, chr, start, end, feature):
        self.chr = chr
        self.start = start
        self.end = end
        self.feature = feature

def parse_records(in_file):
    with open(in_file) as f:
        for line in f:
            yield BedRecord(*line.strip().split())

def merge_groups(group):
    output = [next(group)]
    for record in group:
        if record.start <= output[-1].end:
            output[-1].end = max(output[-1].end, record.end)
        else:
            output.append(record)
    return output

def main(in_file):
    for _, group in groupby(parse_records(in_file), lambda x: (x.chr, x.feature)):
        for r in merge_groups(group):
            print('\t'.join([r.chr, r.start, r.end, r.feature]))

if __name__ == '__main__':
    main('test.bed')