Skip to content

Fixing bug in Chromosome Mapping Tools #701

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

Merged
josemduarte merged 4 commits into from
Aug 22, 2017
Merged

Fixing bug in Chromosome Mapping Tools #701

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
033d92e
CDS length calculation bug fix
valasatava Aug 15, 2017
e5813e8
unit test fix
valasatava Aug 15, 2017
b69ddfb
Setter for chromosome coordinate systems: 0-based, 1-based
valasatava Aug 15, 2017
3e2be54
improvements
valasatava Aug 15, 2017
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
121 changes: 28 additions & 93 deletions biojava-genome/src/main/java/org/biojava/nbio/genome/util/ChromosomeMappingTools.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,11 @@ public class ChromosomeMappingTools {

public static final String CHROMOSOME = "CHROMOSOME";
public static final String CDS = "CDS";


private static int base = 1;
public static void setCoordinateSystem(int baseInt) {
base = baseInt;
}

/** Pretty print the details of a GeneChromosomePosition to a String
*
Expand Down Expand Up @@ -202,25 +206,12 @@ private static String formatExonStructureReverse(GeneChromosomePosition chromPos
*/
public static int getCDSLength(GeneChromosomePosition chromPos) {

logger.debug(chromPos.toString());

logger.debug("chromosomal information: ");

logger.debug("Gene:" + chromPos.getGeneName());
logger.debug(" Transcription (including UTRs): " + chromPos.getTranscriptionStart() + " - " + chromPos.getTranscriptionEnd() + " (length:" + (chromPos.getTranscriptionEnd() - chromPos.getTranscriptionStart()) + ")");
logger.debug(" Orientation: " + chromPos.getOrientation());
logger.debug(" CDS: " + (chromPos.getCdsStart()) + " - " + chromPos.getCdsEnd() + " (length: " + (chromPos.getCdsEnd() - chromPos.getCdsStart()) + ")");


List<Integer> exonStarts = chromPos.getExonStarts();
List<Integer> exonEnds = chromPos.getExonEnds();

logger.debug("Exons:" + exonStarts.size());

int cdsStart = chromPos.getCdsStart();
int cdsEnd = chromPos.getCdsEnd();


int codingLength;
if (chromPos.getOrientation().equals('+'))
codingLength = ChromosomeMappingTools.getCDSLengthForward(exonStarts, exonEnds, cdsStart, cdsEnd);
Expand Down Expand Up @@ -504,16 +495,14 @@ public static ChromPos getChromPosForward(int cdsPos, List<Integer> exonStarts,
*/
public static int getCDSLengthReverse(List<Integer> exonStarts, List<Integer> exonEnds, int cdsStart, int cdsEnd) {

boolean inCoding = false;
int codingLength = 0;

if (cdsEnd < cdsStart) {
int tmp = cdsEnd;
cdsEnd = cdsStart;
cdsStart = tmp;
}

int lengthExons = 0;
cdsStart = cdsStart + base;

// map reverse
for (int i = exonStarts.size() - 1; i >= 0; i--) {
Expand All @@ -526,52 +515,19 @@ public static int getCDSLengthReverse(List<Integer> exonStarts, List<Integer> ex
end = start;
start = tmp;
}
lengthExons += end - start;

if (start <= cdsEnd && end >= cdsEnd) {
inCoding = true;

int tmpstart = start;
if (start < cdsStart) {
tmpstart = cdsStart;
}
codingLength += (cdsEnd - tmpstart);
start = start + base;

boolean debug = logger.isDebugEnabled();

if ( debug) {

StringBuffer b = new StringBuffer();

b.append(" UTR :" + (cdsEnd + 1) + " - " + (end) + newline);
if (tmpstart == start)
b.append(" -> ");
else
b.append(" <-> ");
b.append("Exon :" + tmpstart + " - " + cdsEnd + " | " + (cdsEnd - tmpstart) + " | " + codingLength + " | " + (codingLength % 3) + newline);
// single exon with UTR on both ends
if (tmpstart != start)
b.append(" UTR :" + (cdsStart - 1) + " - " + start + newline);
logger.debug(b.toString());
if ((start < cdsStart && end < cdsStart) || (start > cdsEnd && end > cdsEnd))
continue;

}
} else if (start <= cdsStart && end >= cdsStart) {
inCoding = false;
codingLength += (end - cdsStart);
if (start < cdsStart)
start = cdsStart;

logger.debug(" <- Exon : " + (cdsStart+1) + " - " + end + " | " + (end - cdsStart) + " | " + (codingLength-3) + " | " + (codingLength % 3));
logger.debug(" UTR : " + start + " - " + (cdsStart ));
if (end > cdsEnd)
end = cdsEnd;

} else if (inCoding) {
// full exon is coding
codingLength += (end - start);
logger.debug(" Exon : " + start + " - " + end + " | " + (end - start) + " | " + codingLength + " | " + (codingLength % 3));
} else {
// e.g. see UBQLN3
logger.debug(" no translation!");
}
codingLength += (end - start + 1);
}
logger.debug("length exons: " + lengthExons + " codin length: " + (codingLength - 3));
return codingLength - 3;
}

Expand All @@ -585,47 +541,26 @@ public static int getCDSLengthReverse(List<Integer> exonStarts, List<Integer> ex
* @return
*/
public static int getCDSLengthForward(List<Integer> exonStarts, List<Integer> exonEnds, int cdsStart, int cdsEnd) {
boolean inCoding = false;

int codingLength = 0;

int lengthExons = 0;
// map forward
for (int i = 0; i < exonStarts.size(); i++) {

int start = exonStarts.get(i);
int start = exonStarts.get(i)+base;
int end = exonEnds.get(i);
lengthExons += end - start;

logger.debug("forward exon: " + (start+1) + " - " + end + " | " + (end - start));
if ( (start < cdsStart+base && end < cdsStart) || (start > cdsEnd && end > cdsEnd) )
continue;

if (start+1 <= cdsStart +1 && end >= cdsStart+1) {
if (start < cdsStart+base)
start = cdsStart+base;

inCoding = true;
codingLength += (end - cdsStart);

logger.debug(" UTR : " + start + " - " + (cdsStart ));
logger.debug(" -> Exon : " + (cdsStart+1) + " - " + end + " | " + (end - cdsStart+1) + " | " + codingLength + " | " + (codingLength % 3));

} else if (start+1 <= cdsEnd && end >= cdsEnd) {

inCoding = false;
codingLength += (cdsEnd - start);

logger.debug(" <- Exon : " + (start +1)+ " - " + cdsEnd + " | " + (cdsEnd - start+1) + " | " + codingLength + " | " + (codingLength % 3));
logger.debug(" UTR : " + cdsEnd + 1 + " - " + end);
if (end > cdsEnd)
end = cdsEnd;

} else if (inCoding) {
// full exon is coding
codingLength += (end - start);

logger.debug(" Exon :" + (start+1) + " - " + end + " | " + (end - start+1) + " | " + codingLength + " | " + (codingLength % 3));
}
codingLength += (end - start + 1);
}

logger.debug("length exons: " + Integer.toString(lengthExons));
logger.debug("CDS length:" + Integer.toString((codingLength-3)));

return codingLength-3 ;
return codingLength - 3;
}

/** Extracts the exon boundaries in CDS coordinates. (needs to be divided by 3 to get AA positions)
Expand Down Expand Up @@ -868,7 +803,7 @@ public static int getCDSPosForward(int chromPos, List<Integer> exonStarts, List<
int cdsStart, int cdsEnd) {

// the genetic coordinate is not in a coding region
if ( (chromPos < (cdsStart+1) ) || ( chromPos > (cdsEnd+1) ) ) {
if ( (chromPos < (cdsStart+base) ) || ( chromPos > (cdsEnd+base) ) ) {
logger.debug("The "+format(chromPos)+" position is not in a coding region");
return -1;
}
Expand All @@ -887,7 +822,7 @@ public static int getCDSPosForward(int chromPos, List<Integer> exonStarts, List<

lengthExon = end - start;

if (start+1 <= chromPos && end >= chromPos ) {
if (start+base <= chromPos && end >= chromPos ) {
return codingLength + (chromPos-start);
}
else {
Expand All @@ -914,7 +849,7 @@ public static int getCDSPosReverse(int chromPos, List<Integer> exonStarts, List<
int cdsStart, int cdsEnd) {

// the genetic coordinate is not in a coding region
if ( (chromPos < (cdsStart+1)) || ( chromPos > (cdsEnd+1) ) ) {
if ( (chromPos < (cdsStart+base)) || ( chromPos > (cdsEnd+base) ) ) {
logger.debug("The "+format(chromPos)+" position is not in a coding region");
return -1;
}
Expand All @@ -933,7 +868,7 @@ public static int getCDSPosReverse(int chromPos, List<Integer> exonStarts, List<

lengthExon = end - start;
// +1 offset to be a base 1
if (start+1 <= chromPos && end >= chromPos ) {
if (start+base <= chromPos && end >= chromPos ) {
return codingLength + (end-chromPos+1);
}
else {
Expand Down
61 changes: 61 additions & 0 deletions biojava-genome/src/test/java/org/biojava/nbio/genome/TestChromosomeMappingTools.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
package org.biojava.nbio.genome;

import org.biojava.nbio.genome.util.ChromosomeMappingTools;
import org.junit.Test;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import static org.junit.Assert.assertEquals;

/**
* Created by Yana Valasatava on 8/14/17.
*/
public class TestChromosomeMappingTools {

@Test
public void testGetCDSLengthForward() {

List<Integer> exonStarts = new ArrayList<>(Arrays.asList(10, 30, 50, 70));
List<Integer> exonEnds = new ArrayList<>(Arrays.asList(20, 40, 60, 80));
int cdsStart = 35;
int cdsEnd = 75;

int cdsDesired = 23 - 3;
ChromosomeMappingTools.setCoordinateSystem(0);
int cdsTest = ChromosomeMappingTools.getCDSLengthForward(exonStarts, exonEnds, cdsStart, cdsEnd);

assertEquals(cdsDesired, cdsTest);
}

@Test
public void testGetCDSLengthReverseAsc() {

List<Integer> exonStarts = new ArrayList<>(Arrays.asList(10, 50, 70));
List<Integer> exonEnds = new ArrayList<>(Arrays.asList(20, 60, 80));
int cdsStart = 55;
int cdsEnd = 75;

int cdsDesired = 12 - 3;
ChromosomeMappingTools.setCoordinateSystem(0);
int cdsTest = ChromosomeMappingTools.getCDSLengthReverse(exonStarts, exonEnds, cdsStart, cdsEnd);

assertEquals(cdsDesired, cdsTest);
}

@Test
public void testGetCDSLengthReverseDesc() {

List<Integer> exonStarts = new ArrayList<>(Arrays.asList(70, 50, 10));
List<Integer> exonEnds = new ArrayList<>(Arrays.asList(80, 60, 20));
int cdsStart = 75;
int cdsEnd = 50;

int cdsDesired = 17 - 3;
ChromosomeMappingTools.setCoordinateSystem(0);
int cdsTest = ChromosomeMappingTools.getCDSLengthReverse(exonStarts, exonEnds, cdsStart, cdsEnd);

assertEquals(cdsDesired, cdsTest);
}
}