Our documentation websites are currently offline due to a data center fire. We do not yet have an ETA for restoring service; we’ll update this message when we know more.

Odd behavior in GenomeLocSortedSet class

rkanwarrkanwar Member
edited May 2013 in Ask the GATK team

Hello,

I am trying to use the GenomeLocSortedSet class for a program that I am writing. I have a set of genomic intervals [some of them overlapping]. I am interested in finding a set of intervals which correspond to the union of these intervals. I tried to use GenomeLocSortedSet for this purpose, but am getting odd behavior. Here is a toy example :

public class UnionInterval {

// main
public static void main(String[] args) throws Exception {
    // setup the parser
    ReferenceSequenceFile ref_fa = new CachingIndexedFastaSequenceFile(
            new File("ref.fa"));
    GenomeLocParser parser = new GenomeLocParser(ref_fa);

    // testing out the locset
    GenomeLocSortedSet loc_set = new GenomeLocSortedSet(parser);
    loc_set.addRegion(parser.createGenomeLoc("chr1", 10, 100));
    loc_set.addRegion(parser.createGenomeLoc("chr1", 150, 200));
    System.out.println(loc_set);

    loc_set.addRegion(parser.createGenomeLoc("chr1", 90, 105));
    System.out.println(loc_set);

    loc_set.addRegion(parser.createGenomeLoc("chr1", 110, 115));
    System.out.println(loc_set);

    loc_set.addRegion(parser.createGenomeLoc("chr1", 108, 210));
    System.out.println(loc_set);
}
}

Here is the output I get:

[ chr1:10-100 chr1:150-200]
[ chr1:10-105 chr1:150-200]
[ chr1:10-105 chr1:110-115 chr1:150-200]
[ chr1:10-105 chr1:108-210 chr1:108-210]

The last line has the interval chr1:108-210 repeated 2 times, whereas it should have been repeated only once.
Similarly if you insert a region which overlaps two regions then only one of the regions gets update and not the other. Is this behavior wrong or am I missing something here. Thanks !

regards,
Rahul

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi Rahul,

    Sorry to respond to you so late; unfortunately at this time we simply don't have the resources to help developers with coding problems. We're hoping to improve the dev documentation soon but until then I'm afraid you're on your own. I hope you find the solution to your issue. Good luck!

Sign In or Register to comment.