6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Possible duplicate of Merging bed records based on name $\endgroup$ – Alex Reynolds Mar 20 '18 at 14:59
  • 1
    $\begingroup$ This post is already linked in my entry, and they are looking for different outputs,so not a duplicate $\endgroup$ – bapors Mar 20 '18 at 15:17
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you for your answer. Do you mean that the input file should be sorted by 1st column or 4th? $\endgroup$ – bapors Mar 20 '18 at 9:56
  • 1
    $\begingroup$ You might as well sort it by the first column. You could alternatively rewrite the program to handle files sorted by the fourth column, but the current script isn't written to handle that. $\endgroup$ – Devon Ryan Mar 20 '18 at 9:57
  • $\begingroup$ Sure. The code is clean and working well. I am glad that it is adaptable $\endgroup$ – bapors Mar 20 '18 at 11:01
6
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ Thank you for the answer, it works well. However, would bedtools complain if the $4th column has special characters within? $\endgroup$ – bapors Mar 20 '18 at 11:01
  • $\begingroup$ @bapors why would the 4th column have special characters? $\endgroup$ – terdon Mar 20 '18 at 11:14
  • $\begingroup$ @bapors if you mean that bedtools would complain if the modified 1st field would have special characters, then please edit your question to give a representative example of your input data. That said I just tried using the characters &^%*()@#^ and bedtools had no problem. $\endgroup$ – terdon Mar 20 '18 at 12:24
  • $\begingroup$ No, my input does not have special characters, my question was solely theoretical. Perfect! $\endgroup$ – bapors Mar 20 '18 at 12:25
4
$\begingroup$

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
$\endgroup$
3
$\begingroup$

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')
$\endgroup$
  • $\begingroup$ Thank you for your answer. It assumes that the input will be sorted by which column? $\endgroup$ – bapors Mar 20 '18 at 15:19
  • 1
    $\begingroup$ @bapors effectively by chromosome, then feature then start value, you could easily add code to sort based on these criteria if this assumption is violated $\endgroup$ – Chris_Rands Mar 20 '18 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.