Source code

001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.genome.util;
022
023import com.google.common.collect.Range;
024import org.biojava.nbio.core.sequence.DNASequence;
025import org.biojava.nbio.core.sequence.compound.NucleotideCompound;
026import org.biojava.nbio.core.sequence.template.SequenceView;
027import org.biojava.nbio.genome.parsers.genename.ChromPos;
028import org.biojava.nbio.genome.parsers.genename.GeneChromosomePosition;
029import org.biojava.nbio.genome.parsers.twobit.TwoBitFacade;
030import org.slf4j.Logger;
031import org.slf4j.LoggerFactory;
032
033import java.io.StringWriter;
034import java.util.ArrayList;
035import java.util.List;
036
037/**
038 *  A class that can map chromosomal positions to mRNA (coding sequence) positions.
039 *
040 *  @author Andreas Prlic
041 */
042
043public class ChromosomeMappingTools {
044
045    private static final Logger logger = LoggerFactory.getLogger(ChromosomeMappingTools.class);
046
047    private static final String newline = System.getProperty("line.separator");
048
049    public static final String CHROMOSOME = "CHROMOSOME";
050    public static final String CDS = "CDS";
051
052    private static int base = 1;
053    public static void setCoordinateSystem(int baseInt) {
054        base = baseInt;
055    }
056
057    /** 
058     * Pretty print the details of a GeneChromosomePosition to a String
059     *
060     * @param chromosomePosition
061     * @return
062     */
063    public static String formatExonStructure(GeneChromosomePosition chromosomePosition ){
064        if ( chromosomePosition.getOrientation() == '+')
065            return formatExonStructureForward(chromosomePosition);
066        return formatExonStructureReverse(chromosomePosition);
067    }
068
069    private static String formatExonStructureForward(GeneChromosomePosition chromPos) {
070
071        StringWriter s = new StringWriter();
072
073        List<Integer> exonStarts = chromPos.getExonStarts();
074        List<Integer> exonEnds   = chromPos.getExonEnds();
075
076        int cdsStart = chromPos.getCdsStart();
077        int cdsEnd   = chromPos.getCdsEnd();
078
079        boolean inCoding = false;
080        int codingLength = 0;
081
082        for (int i = 0; i < exonStarts.size(); i++) {
083
084            int start = exonStarts.get(i);
085            int end = exonEnds.get(i);
086
087            if (start <= cdsStart +1 && end >= cdsStart+1) {
088
089                inCoding = true;
090                codingLength += (end - cdsStart);
091                s.append("     UTR         : ").append(format(start)).append(" - ").append(format(cdsStart));
092                s.append(newline);
093                s.append(" ->  Exon        : ").append(format(cdsStart + 1)).append(" - ").append(format(end)).append(" | ").append(Integer.toString(end - cdsStart)).append(" | ").append(Integer.toString(codingLength)).append(" | ").append(Integer.toString(codingLength % 3));
094                s.append(newline);
095
096            } else if (start <= cdsEnd && end >= cdsEnd) {
097                //logger.debug(" <-- CDS end at: " + cdsEnd );
098                inCoding = false;
099                codingLength += (cdsEnd - start);
100
101                s.append(" <-  Exon        : ").append(format(start + 1)).append(" - ").append(format(cdsEnd)).append(" | ").append(Integer.toString(cdsEnd - start)).append(" | ").append(Integer.toString(codingLength)).append(" | ").append(Integer.toString(codingLength % 3));
102                s.append(newline);
103                s.append("     UTR         : " + (cdsEnd +1) + " - " + format(end));
104                s.append(newline);
105
106            } else if (inCoding) {
107                // full exon is coding
108                codingLength += (end - start);
109
110                s.append("     Exon        : ").append(format(start + 1)).append(" - ").append(format(end)).append(" | ").append(Integer.toString(end - start)).append(" | ").append(Integer.toString(codingLength)).append(" | ").append(Integer.toString(codingLength % 3));
111                s.append(newline);
112            }
113        }
114        s.append("Coding Length: ");
115        s.append((codingLength-3)+"");
116        s.append(newline);
117        return s.toString();
118    }
119
120    private static String formatExonStructureReverse(GeneChromosomePosition chromPos) {
121        StringWriter s = new StringWriter();
122
123        List<Integer> exonStarts = chromPos.getExonStarts();
124        List<Integer> exonEnds   = chromPos.getExonEnds();
125
126
127        int cdsStart = chromPos.getCdsStart();
128        int cdsEnd   = chromPos.getCdsEnd();
129
130        // logger.debug("CDS START:" +format(cdsStart) + " - " + format(cdsEnd));
131
132        boolean inCoding = false;
133        int codingLength = 0;
134
135        if (cdsEnd < cdsStart) {
136            int tmp = cdsEnd;
137            cdsEnd = cdsStart;
138            cdsStart = tmp;
139        }
140
141        // map reverse
142        for (int i = exonStarts.size() - 1; i >= 0; i--) {
143
144            int end = exonStarts.get(i);
145            int start = exonEnds.get(i);
146
147            if (end < start) {
148                int tmp = end;
149                end = start;
150                start = tmp;
151            }
152
153            if (start <= cdsEnd && end >= cdsEnd) {
154                inCoding = true;
155
156                int tmpstart = start;
157                if (start < cdsStart) {
158                    tmpstart = cdsStart;
159                }
160                codingLength += (cdsEnd - tmpstart);
161
162                s.append("     UTR         :" + format(cdsEnd + 1) + " | " + format(end));
163                s.append(newline);
164                if (tmpstart == start)
165                    s.append(" ->  ");
166                else
167                    s.append(" <-> ");
168                s.append("Exon        :").append(format(tmpstart + 1)).append(" - ").append(format(cdsEnd)).append(" | ").append(Integer.toString(cdsEnd - tmpstart)).append(" | ").append(Integer.toString(codingLength)).append(" | ").append(Integer.toString(codingLength % 3));
169                s.append(newline);
170                // single exon with UTR on both ends
171                if (tmpstart != start)
172                    s.append("     UTR         :" + format(cdsStart ) + " - " + format(start + 1));
173                s.append(newline);
174
175            } else if (start <= cdsStart && end >= cdsStart) {
176                inCoding = false;
177                codingLength += (end - cdsStart);
178
179                s.append(" <-  Exon        : " + format(cdsStart+1) + " - " + format(end) + " | " + (end - cdsStart) + " | " + codingLength + " | " + (codingLength % 3));
180                s.append(newline);
181                s.append("     UTR         : " + format(start+1) + " - " + format(cdsStart ));
182                s.append(newline);
183
184
185            } else if (inCoding) {
186                // full exon is coding
187                codingLength += (end - start);
188
189                s.append("     Exon        : " + format(start+1) + " - " + format(end) + " | " + (end - start) + " | " + codingLength + " | " + (codingLength % 3));
190                s.append(newline);
191            } else {
192                // e.g. see UBQLN3
193                s.append(" no translation! UTR: ").append(format(start)).append(" - ").append(format(end));
194                s.append(newline);
195            }
196        }
197
198        s.append("CDS length: ").append(Integer.toString(codingLength - 3));
199        s.append(newline);
200
201        return s.toString();
202    }
203
204    /**
205     * Get the length of the CDS in nucleotides.
206     *
207     * @param chromPos
208     * @return length of the CDS in nucleotides.
209     */
210    public static int getCDSLength(GeneChromosomePosition chromPos) {
211
212        List<Integer> exonStarts = chromPos.getExonStarts();
213        List<Integer> exonEnds = chromPos.getExonEnds();
214
215        int cdsStart = chromPos.getCdsStart();
216        int cdsEnd = chromPos.getCdsEnd();
217
218        int codingLength;
219        if (chromPos.getOrientation().equals('+'))
220            codingLength = ChromosomeMappingTools.getCDSLengthForward(exonStarts, exonEnds, cdsStart, cdsEnd);
221        else
222            codingLength = ChromosomeMappingTools.getCDSLengthReverse(exonStarts, exonEnds, cdsStart, cdsEnd);
223        return codingLength;
224    }
225
226    /**
227     * Maps the position of a CDS nucleotide back to the genome
228     *
229     * @param cdsNucleotidePosition
230     * @return a ChromPos object
231     */
232    public static ChromPos getChromosomePosForCDScoordinate(int cdsNucleotidePosition, GeneChromosomePosition chromPos) {
233
234        logger.debug(" ? Checking chromosome position for CDS position " + cdsNucleotidePosition);
235
236        List<Integer> exonStarts = chromPos.getExonStarts();
237        List<Integer> exonEnds = chromPos.getExonEnds();
238
239        logger.debug(" Exons:" + exonStarts.size());
240
241        int cdsStart = chromPos.getCdsStart();
242        int cdsEnd = chromPos.getCdsEnd();
243
244
245        ChromPos chromosomePos = null;
246
247        if (chromPos.getOrientation().equals('+'))
248
249            chromosomePos = ChromosomeMappingTools.getChromPosForward(cdsNucleotidePosition, exonStarts, exonEnds, cdsStart, cdsEnd);
250        else
251            chromosomePos = ChromosomeMappingTools.getChromPosReverse(cdsNucleotidePosition, exonStarts, exonEnds, cdsStart, cdsEnd);
252
253        logger.debug("=> CDS pos " + cdsNucleotidePosition + " for " + chromPos.getGeneName() + " is on chromosome at  " + chromosomePos);
254        return chromosomePos;
255
256    }
257
258    /** 
259     * Returns a nicely formatted representation of the position
260     *
261     * @param chromosomePosition
262     * @return
263     */
264    private static String format(int chromosomePosition){
265        return String.format("%,d", chromosomePosition);
266    }
267
268    /**
269     * Get the CDS position mapped on the chromosome position
270     *
271     * @param exonStarts
272     * @param exonEnds
273     * @param cdsStart
274     * @param cdsEnd
275     * @return
276     */
277    public static ChromPos getChromPosReverse(int cdsPos, List<Integer> exonStarts, List<Integer> exonEnds, int cdsStart, int cdsEnd) {
278
279        boolean inCoding = false;
280        int codingLength = 0;
281
282        if (cdsEnd < cdsStart) {
283            int tmp = cdsEnd;
284            cdsEnd = cdsStart;
285            cdsStart = tmp;
286        }
287
288        int lengthExons = 0;
289
290        // map reverse
291        for (int i = exonStarts.size() - 1; i >= 0; i--) {
292
293            logger.debug("Exon #" + (i+1) + "/" + exonStarts.size());
294            int end = exonStarts.get(i);
295            int start = exonEnds.get(i);
296
297            if (end < start) {
298                int tmp = end;
299                end = start;
300                start = tmp;
301            }
302            lengthExons += end - start;
303
304            logger.debug("     is " + cdsPos + " part of Reverse exon? " + format(start+1) + " - " + format(end) + " | " + (end - start+1));
305            logger.debug("     CDS start: " + format(cdsStart+1) + "-" + format(cdsEnd) + " coding length counter:" + codingLength);
306
307            if (start+1 <= cdsEnd && end >= cdsEnd) {
308
309                // FIRST EXON
310                inCoding = true;
311
312                int tmpstart = start;
313                if (start < cdsStart) {
314                    tmpstart = cdsStart;
315                }
316
317                // here one of the few places where we don't say start+1
318                int check = codingLength + cdsEnd - tmpstart ;
319
320                logger.debug("First Exon    | " + (check) + " | " + format(start+1) + " " + format(end) + " | " + (cdsEnd - tmpstart) + " | " + cdsPos );
321
322
323                if ( ( check > cdsPos)  )  {
324                    int tmp = cdsPos - codingLength ;
325                    logger.debug(" -> found position in UTR exon:  " + format(cdsPos) + " " + format(tmpstart+1) + " tmp:" + format(tmp) + " cs:" + format(cdsStart+1) + " ce:" + format(cdsEnd) + " cl:" + codingLength);
326                    return new ChromPos((cdsEnd - tmp), -1) ;
327                }
328
329                // don't add 1 here
330                codingLength += (cdsEnd - tmpstart );
331
332                boolean debug = logger.isDebugEnabled();
333
334                if ( debug ) {
335
336                    StringBuffer b = new StringBuffer();
337
338                    b.append("     UTR         :" + format(cdsEnd + 1) + " - " + format(end) + newline);
339                    if (tmpstart == start)
340                        b.append(" ->  ");
341                    else
342                        b.append(" <-> ");
343                    b.append("Exon        :" + format(tmpstart + 1) + " - " + (cdsEnd) + " | " + format(cdsEnd - tmpstart + 1) + " - " + codingLength + " | " + (codingLength % 3) + newline);
344
345                    // single exon with UTR on both ends
346                    if (tmpstart != start)
347                        b.append("     UTR         :" + format(cdsStart) + " - " + format(start + 1) + newline);
348
349                    logger.debug(b.toString());
350                }
351            } else if (start <= cdsStart && end >= cdsStart) {
352
353                // LAST EXON
354                inCoding = false;
355
356                if (codingLength + end - cdsStart >= cdsPos) {
357
358                    // how many remaining coding nucleotides?
359                    int tmp = codingLength + end - cdsStart - cdsPos ;
360
361                    logger.debug("cdl: " +codingLength + " tmp:" + tmp + " cdsStart: " + format(cdsStart));
362
363                    logger.debug(" -> XXX found position noncoding exon:  cdsPos:" + cdsPos + " s:" + format(start + 1) + " tmp:" + format(tmp) + " cdsStart:" + (cdsStart + 1) + " codingLength:" + codingLength + " cdsEnd:" + format(cdsEnd));
364
365                    return new ChromPos((cdsStart + tmp),-1);
366                }
367
368                codingLength += (end - cdsStart);
369
370                logger.debug(" <-  Exon        : " + format(cdsStart+1) + " - " + format(end) + " | " + format(end - cdsStart+1) + " | " + codingLength + " | " + (codingLength % 3));
371                logger.debug("     UTR         : " + format(start+1) + " - " + format(cdsStart ));
372
373            } else if (inCoding) {
374
375                if (codingLength + end - start -1  >= cdsPos) {
376
377                    int tmp = cdsPos - codingLength ;
378
379                    if ( tmp > (end - start ) ) {
380                        tmp = (end - start );
381                        logger.debug("changing tmp to " + tmp);
382                    }
383                    logger.debug("     " + cdsPos + " " + codingLength + " | " + (cdsPos - codingLength) + " | " + (end -start) + " | " + tmp);
384                    logger.debug("     Exon        : " + format(start+1) + " - " + format(end) + " | " + format(end - start) + " | " + codingLength + " | " + (codingLength % 3));
385                    logger.debug(" ->  RRR found position coding exon:  " + cdsPos + " " + format(start+1) + " " + format(end) + " " + tmp + " " + format(cdsStart+1) + " " + codingLength);
386
387                    return new ChromPos((end - tmp),cdsPos %3);
388                }
389                // full exon is coding
390                codingLength += (end - start) ;
391
392                logger.debug("     Exon        : " + format(start+1) + " - " + format(end) + " | " + format(end - start+1) + " | " + codingLength + " | " + (codingLength % 3));
393            } else {
394                // e.g. see UBQLN3
395                logger.debug(" no translation!");
396            }
397            logger.debug("     coding length: " + codingLength + "(phase:" + (codingLength % 3) + ") CDS POS trying to map:" + cdsPos);
398        }
399
400        logger.debug("length exons: " + lengthExons);
401        // could not map, or map over the full length??
402        return new ChromPos(-1,-1);
403
404    }
405
406    /**
407     * Get the CDS position mapped onto the chromosome position
408     *
409     * @param exonStarts
410     * @param exonEnds
411     * @param cdsStart
412     * @param cdsEnd
413     * @return
414     */
415    public static ChromPos getChromPosForward(int cdsPos, List<Integer> exonStarts, List<Integer> exonEnds, int cdsStart, int cdsEnd) {
416        boolean inCoding = false;
417        int codingLength = 0;
418
419        @SuppressWarnings("unused")
420                int lengthExons = 0;
421        // map forward
422        for (int i = 0; i < exonStarts.size(); i++) {
423
424            // start can include UTR
425            int start = exonStarts.get(i);
426            int end = exonEnds.get(i);
427
428            lengthExons += end - start;
429
430            if (start <= cdsStart +1 && end >= cdsStart+1) {
431                // first exon with UTR
432                if (codingLength + (end - cdsStart-1) >= cdsPos) {
433                    // we are reaching our target position
434                    int tmp = cdsPos - codingLength;
435
436
437                    logger.debug(cdsStart + " | " + codingLength + " | " + tmp);
438                    logger.debug(" -> found position in UTR exon:  #"+(i+1)+ " cdsPos:" + cdsPos +
439                            " return:"+(cdsStart +1 + tmp) +" start:" + format(start + 1) + " " + format(tmp) + " " + cdsStart + " " + codingLength);
440
441                    // we start 1 after cdsStart...
442                    return new ChromPos((cdsStart +1 + tmp),-1);
443                }
444                inCoding = true;
445                codingLength += (end - cdsStart);
446
447                logger.debug("     UTR         : " + format(start+1) + " - " + (cdsStart ));
448                logger.debug(" ->  Exon        : " + format(cdsStart+1) + " - " + format(end) + " | " + format(end - cdsStart) + " | " + codingLength + " | " + (codingLength % 3));
449
450            } else if (start+1 <= cdsEnd && end >= cdsEnd) {
451                // LAST EXON with UTR
452                //logger.debug(" <-- CDS end at: " + cdsEnd );
453                inCoding = false;
454                if (codingLength + (cdsEnd - start-1) >= cdsPos) {
455                    int tmp = cdsPos - codingLength;
456
457                    logger.debug(" <-  Exon        : " + format(start+1) + " - " + format(cdsEnd) + " | " + format(cdsEnd - start) + " | " + codingLength + " | " + (codingLength % 3));
458                    logger.debug("     UTR         : " + format(cdsEnd + 1) + " - " + format(end));
459                    logger.debug( codingLength + " | " + tmp + " | " + format(start+1));
460                    logger.debug(" -> chromPosForward found position in non coding exon:  " + cdsPos + " " + format(start+1) + " " + format(tmp) + " " + format(cdsStart) + " " + codingLength);
461
462                    return new ChromPos((start +1 + tmp),cdsPos%3);
463                }
464                codingLength += (cdsEnd - start-1);
465
466                logger.debug(" <-  Exon        : " + format(start+1) + " - " + format(cdsEnd) + " | " + format(cdsEnd - start) + " | " + codingLength + " | " + (codingLength % 3));
467                logger.debug("     UTR         : " + format(cdsEnd + 1) + " - " + format(end));
468
469
470            } else if (inCoding) {
471                // A standard coding Exon
472                // tests for the maximum length of this coding exon
473                if (codingLength + (end - start -1)  >= cdsPos) {
474
475                    // we are within the range of this exon
476                    int tmp = cdsPos - codingLength ;
477
478                    logger.debug("     Exon        : " + format(start+1) + " - " + format(end) + " | " + format(end - start) + " | " + tmp + " | " + codingLength);
479                    logger.debug(" -> found chr position in coding exon #" + (i+1) + ":  cdsPos:" + format(cdsPos) + " s:" + format(start) + "-" + format(end) + " tmp:" + format(tmp) + " cdsStart:" + format(cdsStart) + " codingLength:" + codingLength);
480
481                    return new ChromPos((start +1 + tmp),cdsPos%3);
482                }
483                // full exon is coding
484                codingLength += (end - start );
485
486                logger.debug("     Exon        : " + format(start+1) + " - " + format(end) + " | " + format(end - start) + " | " + codingLength + " | " + (codingLength % 3));
487            }
488        }
489        return new ChromPos(-1,-1);
490    }
491
492    /**
493     * Get the length of the coding sequence
494     *
495     * @param exonStarts
496     * @param exonEnds
497     * @param cdsStart
498     * @param cdsEnd
499     * @return
500     */
501    public static int getCDSLengthReverse(List<Integer> exonStarts, List<Integer> exonEnds, int cdsStart, int cdsEnd) {
502
503        int codingLength = 0;
504
505        if (cdsEnd < cdsStart) {
506            int tmp = cdsEnd;
507            cdsEnd = cdsStart;
508            cdsStart = tmp;
509        }
510        cdsStart = cdsStart + base;
511
512        // map reverse
513        for (int i = exonStarts.size() - 1; i >= 0; i--) {
514
515            int end = exonStarts.get(i);
516            int start = exonEnds.get(i);
517
518            if (end < start) {
519                int tmp = end;
520                end = start;
521                start = tmp;
522            }
523            start = start + base;
524
525            if ((start < cdsStart && end < cdsStart) || (start > cdsEnd && end > cdsEnd))
526                continue;
527
528            if (start < cdsStart)
529                start = cdsStart;
530
531            if (end > cdsEnd)
532                end = cdsEnd;
533
534            codingLength += (end - start + 1);
535        }
536        return codingLength - 3;
537    }
538
539    /**
540     * Get the length of the coding sequence
541     *
542     * @param exonStarts
543     * @param exonEnds
544     * @param cdsStart
545     * @param cdsEnd
546     * @return
547     */
548    public static int getCDSLengthForward(List<Integer> exonStarts, List<Integer> exonEnds, int cdsStart, int cdsEnd) {
549
550        int codingLength = 0;
551
552        for (int i = 0; i < exonStarts.size(); i++) {
553
554            int start = exonStarts.get(i)+base;
555            int end = exonEnds.get(i);
556
557            if ( (start < cdsStart+base && end < cdsStart) || (start > cdsEnd && end > cdsEnd) )
558                continue;
559
560            if (start < cdsStart+base)
561                start = cdsStart+base;
562
563            if (end > cdsEnd)
564                end = cdsEnd;
565
566            codingLength += (end - start + 1);
567        }
568        return codingLength - 3;
569    }
570
571    /** 
572     * Extracts the exon boundaries in CDS coordinates. (needs to be divided by 3 to get AA positions)
573     *
574     * @param chromPos
575     * @return
576     */
577    public static List<Range<Integer>> getCDSExonRanges(GeneChromosomePosition chromPos){
578        if ( chromPos.getOrientation() == '+')
579            return getCDSExonRangesForward(chromPos,CDS);
580        return getCDSExonRangesReverse(chromPos,CDS);
581    }
582
583    /** Extracts the boundaries of the coding regions in chromosomal coordinates
584     *
585     * @param chromPos
586     * @return
587     */
588    public static List<Range<Integer>> getChromosomalRangesForCDS(GeneChromosomePosition chromPos){
589        if ( chromPos.getOrientation() == '+')
590            return getCDSExonRangesForward(chromPos,CHROMOSOME);
591        return getCDSExonRangesReverse(chromPos,CHROMOSOME);
592    }
593
594    private static List<Range<Integer>> getCDSExonRangesReverse(GeneChromosomePosition chromPos, String responseType) {
595
596        List<Integer> exonStarts = chromPos.getExonStarts();
597        List<Integer> exonEnds   = chromPos.getExonEnds();
598
599        List<Range<Integer>> data = new ArrayList<>();
600        int cdsStart = chromPos.getCdsStart();
601        int cdsEnd   = chromPos.getCdsEnd();
602
603        boolean inCoding = false;
604        int codingLength = 0;
605
606        if (cdsEnd < cdsStart) {
607            int tmp = cdsEnd;
608            cdsEnd = cdsStart;
609            cdsStart = tmp;
610        }
611
612        java.lang.StringBuffer s =null;
613
614        boolean debug = logger.isDebugEnabled();
615
616        if ( debug)
617            s = new StringBuffer();
618
619                //int lengthExons = 0;
620
621        // map reverse
622        for (int i = exonStarts.size() - 1; i >= 0; i--) {
623
624            int end = exonStarts.get(i);
625            int start = exonEnds.get(i);
626
627            if (end < start) {
628                int tmp = end;
629                end = start;
630                start = tmp;
631            }
632            //lengthExons += end - start;
633            //s.append("Reverse exon: " + end + " - " + start + " | " + (end - start));
634            //s.append(newline);
635
636            if (start <= cdsEnd && end >= cdsEnd) {
637                inCoding = true;
638
639                int tmpstart = start;
640                if (start < cdsStart) {
641                    tmpstart = cdsStart;
642                }
643                codingLength += (cdsEnd - tmpstart);
644                if ( debug ) {
645                    s.append("     UTR         :").append(format(cdsEnd + 1)).append(" | ").append(format(end));
646                    s.append(newline);
647                    if (tmpstart == start)
648                        s.append(" ->  ");
649                    else
650                        s.append(" <-> ");
651                    s.append("Exon        :").append(format(tmpstart + 1)).append(" - ").append(format(cdsEnd)).append(" | ").append(cdsEnd - tmpstart).append(" | ").append(codingLength).append(" | ").append(codingLength % 3);
652                    s.append(newline);
653                    // single exon with UTR on both ends
654                    if (tmpstart != start)
655                        s.append("     UTR         :").append(format(cdsStart)).append(" - ").append(format(start + 1));
656                    s.append(newline);
657                }
658
659
660                Range<Integer> r ;
661                if ( responseType.equals(CDS))
662                    r = Range.closed(0,codingLength);
663                else
664                    r = Range.closed(tmpstart,cdsEnd);
665
666                data.add(r);
667
668            } else if (start <= cdsStart && end >= cdsStart) {
669                inCoding = false;
670
671                Range<Integer> r;
672                if ( responseType.equals(CDS))
673                    r = Range.closed(codingLength,codingLength+(end-cdsStart));
674                else
675                    r = Range.closed(cdsStart+1,end);
676
677                data.add(r);
678
679                codingLength += (end - cdsStart);
680                if  (debug) {
681                    s.append(" <-  Exon        : " + format(cdsStart + 1) + " - " + format(end) + " | " + (end - cdsStart) + " | " + codingLength + " | " + (codingLength % 3));
682                    s.append(newline);
683                    s.append("     UTR         : ").append(format(start + 1)).append(" - ").append(format(cdsStart));
684                    s.append(newline);
685                }
686            } else if (inCoding) {
687                // full exon is coding
688                Range<Integer> r;
689                if ( responseType.equals(CDS))
690                    r = Range.closed(codingLength,codingLength+(end-start));
691                else
692                    r = Range.closed(start,end);
693                data.add(r);
694
695                codingLength += (end - start);
696                if  (debug) {
697                    s.append("     Exon        : " + format(start + 1) + " - " + format(end) + " | " + (end - start) + " | " + codingLength + " | " + (codingLength % 3));
698                    s.append(newline);
699                }
700            } else {
701                // e.g. see UBQLN3
702                if ( debug ) {
703                    s.append(" no translation! UTR: " + format(start) + " - " + format(end));
704                    s.append(newline);
705                }
706            }
707        }
708        if ( debug ) {
709            s.append("CDS length: ").append(Integer.toString(codingLength - 3));
710            s.append(newline);
711            logger.debug(s.toString());
712        }
713
714        return data;
715    }
716
717    private static List<Range<Integer>> getCDSExonRangesForward(GeneChromosomePosition chromPos, String responseType) {
718
719        List<Range<Integer>> data = new ArrayList<>();
720        List<Integer> exonStarts = chromPos.getExonStarts();
721        List<Integer> exonEnds   = chromPos.getExonEnds();
722
723        int cdsStart = chromPos.getCdsStart();
724        int cdsEnd   = chromPos.getCdsEnd();
725
726        boolean inCoding = false;
727        int codingLength = 0;
728
729        for (int i = 0; i < exonStarts.size(); i++) {
730
731            int start = exonStarts.get(i);
732            int end = exonEnds.get(i);
733
734            if (start <= cdsStart && end >= cdsStart) {
735
736                inCoding = true;
737                codingLength += (end - cdsStart);
738
739                Range<Integer> r;
740                if ( responseType.equals(CDS))
741                    r = Range.closed(0,codingLength);
742                else
743                    r = Range.closed(cdsStart,end);
744                data.add(r);
745
746            } else if (start <= cdsEnd && end >= cdsEnd) {
747                //logger.debug(" <-- CDS end at: " + cdsEnd );
748                inCoding = false;
749
750                Range<Integer> r;
751                if ( responseType.equals(CDS))
752                    r = Range.closed(codingLength,codingLength+(cdsEnd-start));
753                else
754                    r = Range.closed(start,cdsEnd);
755                data.add(r);
756                codingLength += (cdsEnd - start);
757
758            } else if (inCoding) {
759                // full exon is coding
760                Range<Integer> r;
761                if ( responseType.equals(CDS))
762                    r = Range.closed(codingLength,codingLength+(end-start));
763                else
764                    r = Range.closed(start,end);
765                data.add(r);
766                codingLength += (end - start);
767            }
768        }
769        return data;
770    }
771
772    /**
773     * I have a genomic coordinate, where is it on the mRNA
774     *
775     * @param coordinate
776     * @param chromosomePosition
777     * @return
778     */
779    public static int getCDSPosForChromosomeCoordinate(int coordinate, GeneChromosomePosition chromosomePosition) {
780
781        if ( chromosomePosition.getOrientation() == '+')
782                return getCDSPosForward(coordinate,
783                    chromosomePosition.getExonStarts(),
784                    chromosomePosition.getExonEnds(),
785                    chromosomePosition.getCdsStart(),
786                    chromosomePosition.getCdsEnd());
787
788        return getCDSPosReverse(coordinate,
789                chromosomePosition.getExonStarts(),
790                chromosomePosition.getExonEnds(),
791                chromosomePosition.getCdsStart(),
792                chromosomePosition.getCdsEnd());
793    }
794    
795        /** 
796         * Converts the genetic coordinate to the position of the nucleotide on the mRNA sequence for a gene 
797         * living on the forward DNA strand.
798         * 
799         * @param chromPos The genetic coordinate on a chromosome 
800     * @param exonStarts The list holding the genetic coordinates pointing to the start positions of the exons (including UTR regions)  
801     * @param exonEnds The list holding the genetic coordinates pointing to the end positions of the exons (including UTR regions)
802     * @param cdsStart The start position of a coding region
803     * @param cdsEnd The end position of a coding region
804     * 
805     * @return the position of the nucleotide base on the mRNA sequence corresponding to the input genetic coordinate (base 1)
806         * 
807         * @author Yana Valasatava
808         */
809    public static int getCDSPosForward(int chromPos, List<Integer> exonStarts, List<Integer> exonEnds,
810            int cdsStart, int cdsEnd) {
811        
812        // the genetic coordinate is not in a coding region
813        if ( (chromPos < (cdsStart+base) ) || ( chromPos > (cdsEnd+base) ) ) {
814                logger.debug("The "+format(chromPos)+" position is not in a coding region");
815            return -1;
816        }
817        
818        logger.debug("looking for CDS position for " +format(chromPos));
819        
820        // map the genetic coordinates of coding region on a stretch of a reverse strand
821        List<Range<Integer>> cdsRegions = getCDSRegions(exonStarts, exonEnds, cdsStart, cdsEnd);
822        
823        int codingLength = 0;
824        int lengthExon = 0;
825        for (Range<Integer> range : cdsRegions) {
826                
827                    int start = range.lowerEndpoint();
828                    int end = range.upperEndpoint();
829                    
830                    lengthExon = end - start;
831
832                    if (start+base <= chromPos && end >= chromPos ) {
833                        return codingLength + (chromPos-start);
834                    }
835                else { 
836                        codingLength += lengthExon;
837                }
838        }
839        return -1;
840    }
841    
842        /** 
843         * Converts the genetic coordinate to the position of the nucleotide on the mRNA sequence for a gene 
844         * living on the reverse DNA strand.
845         * 
846         * @param chromPos The genetic coordinate on a chromosome 
847     * @param exonStarts The list holding the genetic coordinates pointing to the start positions of the exons (including UTR regions)  
848     * @param exonEnds The list holding the genetic coordinates pointing to the end positions of the exons (including UTR regions)
849     * @param cdsStart The start position of a coding region
850     * @param cdsEnd The end position of a coding region
851     * 
852     * @return the position of the nucleotide base on the mRNA sequence corresponding to the input genetic coordinate (base 1)
853         * 
854         * @author Yana Valasatava
855         */
856    public static int getCDSPosReverse(int chromPos, List<Integer> exonStarts, List<Integer> exonEnds,
857            int cdsStart, int cdsEnd) {
858        
859        // the genetic coordinate is not in a coding region
860        if ( (chromPos < (cdsStart+base)) || ( chromPos > (cdsEnd+base) ) ) {
861                logger.debug("The "+format(chromPos)+" position is not in a coding region");
862            return -1;
863        }
864        
865        logger.debug("looking for CDS position for " +format(chromPos));
866                
867        // map the genetic coordinate on a stretch of a reverse strand
868        List<Range<Integer>> cdsRegions = getCDSRegions(exonStarts, exonEnds, cdsStart, cdsEnd);
869        
870        int codingLength = 0;
871        int lengthExon = 0;
872        for ( int i=cdsRegions.size()-1; i>=0; i-- ) {
873                
874                    int start = cdsRegions.get(i).lowerEndpoint();
875                    int end = cdsRegions.get(i).upperEndpoint();
876                    
877                    lengthExon = end - start;
878                    // +1 offset to be a base 1
879                    if (start+base <= chromPos && end >= chromPos ) {
880                        return codingLength + (end-chromPos+1);
881                    }
882                else { 
883                        codingLength += lengthExon;
884                }
885        }
886        return -1;
887    }
888    
889    /** 
890     * Extracts the exons boundaries in CDS coordinates corresponding to the forward DNA strand.
891     *
892     * @param origExonStarts The list holding the genetic coordinates pointing to the start positions of the exons (including UTR regions)
893     * @param origExonEnds The list holding the genetic coordinates pointing to the end positions of the exons (including UTR regions)
894     * @param cdsStart The start position of a coding region
895     * @param cdsEnd The end position of a coding region
896     * 
897     * @return the list of genetic positions corresponding to the exons boundaries in CDS coordinates
898     */
899    public static List<Range<Integer>> getCDSRegions(List<Integer> origExonStarts, List<Integer> origExonEnds, int cdsStart, int cdsEnd) {
900        
901        // remove exons that are fully landed in UTRs
902        List<Integer> exonStarts = new ArrayList<Integer>(origExonStarts);
903        List<Integer> exonEnds = new ArrayList<Integer>(origExonEnds);
904        
905        int j=0;
906        for (int i = 0; i < origExonStarts.size(); i++) {
907                if ( ( origExonEnds.get(i) < cdsStart) || ( origExonStarts.get(i) > cdsEnd) ) {
908                        exonStarts.remove(j);
909                        exonEnds.remove(j);
910                }
911                else {
912                        j++;
913                }
914        }
915        
916        // remove untranslated regions from exons
917        int nExons = exonStarts.size();
918        exonStarts.remove(0);
919        exonStarts.add(0, cdsStart);
920        exonEnds.remove(nExons-1);
921        exonEnds.add(cdsEnd);
922        
923        List<Range<Integer>> cdsRegion = new ArrayList<Range<Integer>>();
924        for ( int i=0; i<nExons; i++ ) {
925                Range<Integer> r = Range.closed(exonStarts.get(i), exonEnds.get(i));
926                cdsRegion.add(r);
927        }
928                return cdsRegion;
929    }
930    
931    /** 
932     * Extracts the DNA sequence transcribed from the input genetic coordinates.
933     *
934     * @param twoBitFacade the facade that provide an access to a 2bit file
935     * @param gcp The container with chromosomal positions
936     * 
937     * @return the DNA sequence transcribed from the input genetic coordinates
938     */
939    public static DNASequence getTranscriptDNASequence(TwoBitFacade twoBitFacade, GeneChromosomePosition gcp) throws Exception {
940        return getTranscriptDNASequence(twoBitFacade,gcp.getChromosome(),gcp.getExonStarts(), gcp.getExonEnds(), gcp.getCdsStart(), gcp.getCdsEnd(), gcp.getOrientation());
941    }
942    
943    /** 
944     * Extracts the DNA sequence transcribed from the input genetic coordinates.
945     *
946     * @param chromosome the name of the chromosome
947     * @param exonStarts The list holding the genetic coordinates pointing to the start positions of the exons (including UTR regions)  
948     * @param exonEnds The list holding the genetic coordinates pointing to the end positions of the exons (including UTR regions)
949     * @param cdsStart The start position of a coding region
950     * @param cdsEnd The end position of a coding region
951     * @param orientation The orientation of the strand where the gene is living
952     * 
953     * @return the DNA sequence transcribed from the input genetic coordinates
954     */
955        public static DNASequence getTranscriptDNASequence(TwoBitFacade twoBitFacade, String chromosome, List<Integer> exonStarts, List<Integer> exonEnds, int cdsStart, int cdsEnd, Character orientation) throws Exception {
956
957                List<Range<Integer>> cdsRegion = getCDSRegions(exonStarts, exonEnds, cdsStart, cdsEnd);
958
959                String dnaSequence = "";
960                for (Range<Integer> range : cdsRegion) {
961                        String exonSequence = twoBitFacade.getSequence(chromosome,range.lowerEndpoint(), range.upperEndpoint());
962            dnaSequence += exonSequence;
963                }
964                if (orientation.equals('-')) {
965            dnaSequence = new StringBuilder(dnaSequence).reverse().toString();
966                        DNASequence dna = new DNASequence(dnaSequence);
967                        SequenceView<NucleotideCompound> compliment = dna.getComplement();
968            dnaSequence = compliment.getSequenceAsString();
969                }
970                return new DNASequence(dnaSequence.toUpperCase());
971        }
972}