Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simplify BedToIntervalList by not reimplementing coordinate conversion #1292

Merged
merged 2 commits into from
Mar 12, 2019
Merged

Simplify BedToIntervalList by not reimplementing coordinate conversion #1292

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d5efb57
Simplify BedToIntervalList by not reimplementing coordiant conversion
lbergelson Mar 7, 2019
eddfc69
add more tests and fix for zero-length edge case at the start of the …
lbergelson Mar 11, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 4 additions & 14 deletions src/main/java/picard/util/BedToIntervalList.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,25 +124,14 @@ protected int doWork() {
header.setSortOrder(SAMFileHeader.SortOrder.coordinate);
final IntervalList intervalList = new IntervalList(header);

/**
* NB: BED is zero-based, but a BEDCodec by default (since it is returns tribble Features) has an offset of one,
* so it returns 1-based starts. Ugh. Set to zero.
*/
final FeatureReader<BEDFeature> bedReader = AbstractFeatureReader.getFeatureReader(INPUT.getAbsolutePath(), new BEDCodec(BEDCodec.StartOffset.ZERO), false);
final FeatureReader<BEDFeature> bedReader = AbstractFeatureReader.getFeatureReader(INPUT.getAbsolutePath(), new BEDCodec(), false);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

will this be able to read a bed that has position 0?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(the tests are not very good...they just test that the dictionary is right and that the program didn't crash...)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good question.. The tests do actually test the output interval see testBedToIntervalListDataProvider, that's how we noticed that there was a bug with merging 0-length intervals.

With some experimentation: 0 intervals that start at 0 work, 0-length intervals that start as 0 have never worked... I added a fix and a new test.

To be complete we might want to allow the weird edge case of the zero-length interval AFTER the end of the contig... My guess is that doesn't work either... GATK will almost certainly explode if it sees one of those as well.

final CloseableTribbleIterator<BEDFeature> iterator = bedReader.iterator();
final ProgressLogger progressLogger = new ProgressLogger(LOG, (int) 1e6);

while (iterator.hasNext()) {
final BEDFeature bedFeature = iterator.next();
final String sequenceName = bedFeature.getContig();
/**
* NB: BED is zero-based, so we need to add one here to make it one-based. Please observe we set the start
* offset to zero when creating the BEDCodec.
*/
final int start = bedFeature.getStart() + 1;
/**
* NB: BED is 0-based OPEN (which, for the end is equivalent to 1-based closed).
*/
final int start = bedFeature.getStart();
final int end = bedFeature.getEnd();
// NB: do not use an empty name within an interval
String name = bedFeature.getName();
Expand All @@ -157,7 +146,8 @@ protected int doWork() {
throw new PicardException(String.format("Start on sequence '%s' was less than one: %d", sequenceName, start));
} else if (sequenceRecord.getSequenceLength() < start) {
throw new PicardException(String.format("Start on sequence '%s' was past the end: %d < %d", sequenceName, sequenceRecord.getSequenceLength(), start));
} else if (end < 1) {
} else if ((end == 0 && start != 1 ) //special case for 0-length interval at the start of a contig
|| end < 0 ) {
throw new PicardException(String.format("End on sequence '%s' was less than one: %d", sequenceName, end));
} else if (sequenceRecord.getSequenceLength() < end) {
throw new PicardException(String.format("End on sequence '%s' was past the end: %d < %d", sequenceName, sequenceRecord.getSequenceLength(), end));
Expand Down
6 changes: 5 additions & 1 deletion src/test/java/picard/util/BedToIntervalListTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,11 @@ public Object[][] testBedToIntervalListDataProvider() {
{"overlapping.bed"},
{"extended.bed"},
{"one_base_interval.bed"},
{"zero_base_interval.bed"}
{"zero_base_interval.bed"},
{"first_base_in_contig.bed"},
{"zero_length_interval_at_first_position_in_contig.bed"},
{"last_base_in_contig.bed"},
{"multi_contig.bed"}
};
}

Expand Down
1 change: 1 addition & 0 deletions testdata/picard/util/BedToIntervalListTest/first_base_in_contig.bed
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
chr1 0 100
10 changes: 10 additions & 0 deletions testdata/picard/util/BedToIntervalListTest/first_base_in_contig.bed.interval_list
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:1000000
@SQ SN:chr2 LN:1000000
@SQ SN:chr3 LN:1000000
@SQ SN:chr4 LN:1000000
@SQ SN:chr5 LN:1000000
@SQ SN:chr6 LN:1000000
@SQ SN:chr7 LN:1000000
@SQ SN:chr8 LN:1000000
chr1 1 100 + .
1 change: 1 addition & 0 deletions testdata/picard/util/BedToIntervalListTest/last_base_in_contig.bed
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
chr1 0 1000000
10 changes: 10 additions & 0 deletions testdata/picard/util/BedToIntervalListTest/last_base_in_contig.bed.interval_list
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:1000000
@SQ SN:chr2 LN:1000000
@SQ SN:chr3 LN:1000000
@SQ SN:chr4 LN:1000000
@SQ SN:chr5 LN:1000000
@SQ SN:chr6 LN:1000000
@SQ SN:chr7 LN:1000000
@SQ SN:chr8 LN:1000000
chr1 1 1000000 + .
4 changes: 4 additions & 0 deletions testdata/picard/util/BedToIntervalListTest/multi_contig.bed
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
chr1 100 2000 chr1_100_2000+ 11 +
chr1 3000 4000 chr1_3000_4000- 12 -
chr2 100 2000 chr2_100_2000+ 11 +
chr2 3000 4000 chr2_3000_4000- 12 -
13 changes: 13 additions & 0 deletions testdata/picard/util/BedToIntervalListTest/multi_contig.bed.interval_list
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:1000000
@SQ SN:chr2 LN:1000000
@SQ SN:chr3 LN:1000000
@SQ SN:chr4 LN:1000000
@SQ SN:chr5 LN:1000000
@SQ SN:chr6 LN:1000000
@SQ SN:chr7 LN:1000000
@SQ SN:chr8 LN:1000000
chr1 101 2000 + chr1_100_2000+
chr1 3001 4000 - chr1_3000_4000-
chr2 101 2000 + chr2_100_2000+
chr2 3001 4000 - chr2_3000_4000-
1 change: 1 addition & 0 deletions ...ta/picard/util/BedToIntervalListTest/zero_length_interval_at_first_position_in_contig.bed
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
chr1 0 0
10 changes: 10 additions & 0 deletions .../BedToIntervalListTest/zero_length_interval_at_first_position_in_contig.bed.interval_list
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:1000000
@SQ SN:chr2 LN:1000000
@SQ SN:chr3 LN:1000000
@SQ SN:chr4 LN:1000000
@SQ SN:chr5 LN:1000000
@SQ SN:chr6 LN:1000000
@SQ SN:chr7 LN:1000000
@SQ SN:chr8 LN:1000000
chr1 1 0 + .