test/jalview/analysis/AlignmentUtilsTests.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.analysis;
  22
  23 import static org.testng.AssertJUnit.assertEquals;
  24 import static org.testng.AssertJUnit.assertFalse;
  25 import static org.testng.AssertJUnit.assertNotNull;
  26 import static org.testng.AssertJUnit.assertNull;
  27 import static org.testng.AssertJUnit.assertSame;
  28 import static org.testng.AssertJUnit.assertTrue;
  29
  30 import jalview.analysis.AlignmentUtils.DnaVariant;
  31 import jalview.datamodel.AlignedCodonFrame;
  32 import jalview.datamodel.Alignment;
  33 import jalview.datamodel.AlignmentAnnotation;
  34 import jalview.datamodel.AlignmentI;
  35 import jalview.datamodel.Annotation;
  36 import jalview.datamodel.DBRefEntry;
  37 import jalview.datamodel.GeneLociI;
  38 import jalview.datamodel.Mapping;
  39 import jalview.datamodel.SearchResultMatchI;
  40 import jalview.datamodel.SearchResultsI;
  41 import jalview.datamodel.Sequence;
  42 import jalview.datamodel.SequenceFeature;
  43 import jalview.datamodel.SequenceI;
  44 import jalview.datamodel.features.SequenceFeatures;
  45 import jalview.gui.JvOptionPane;
  46 import jalview.io.AppletFormatAdapter;
  47 import jalview.io.DataSourceType;
  48 import jalview.io.FileFormat;
  49 import jalview.io.FileFormatI;
  50 import jalview.io.FormatAdapter;
  51 import jalview.io.gff.SequenceOntologyI;
  52 import jalview.util.MapList;
  53 import jalview.util.MappingUtils;
  54 import jalview.ws.params.InvalidArgumentException;
  55
  56 import java.io.IOException;
  57 import java.util.ArrayList;
  58 import java.util.Arrays;
  59 import java.util.LinkedHashMap;
  60 import java.util.List;
  61 import java.util.Map;
  62 import java.util.TreeMap;
  63
  64 import org.testng.annotations.BeforeClass;
  65 import org.testng.annotations.Test;
  66
  67 public class AlignmentUtilsTests
  68 {
  69   private static Sequence ts = new Sequence("short",
  70           "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklm");
  71
  72   @BeforeClass(alwaysRun = true)
  73   public void setUpJvOptionPane()
  74   {
  75     JvOptionPane.setInteractiveMode(false);
  76     JvOptionPane.setMockResponse(JvOptionPane.CANCEL_OPTION);
  77   }
  78
  79   @Test(groups = { "Functional" })
  80   public void testExpandContext()
  81   {
  82     AlignmentI al = new Alignment(new Sequence[] {});
  83     for (int i = 4; i < 14; i += 2)
  84     {
  85       SequenceI s1 = ts.deriveSequence().getSubSequence(i, i + 7);
  86       al.addSequence(s1);
  87     }
  88     System.out.println(new AppletFormatAdapter()
  89             .formatSequences(FileFormat.Clustal, al, true));
  90     for (int flnk = -1; flnk < 25; flnk++)
  91     {
  92       AlignmentI exp = AlignmentUtils.expandContext(al, flnk);
  93       System.out.println("\nFlank size: " + flnk);
  94       System.out.println(new AppletFormatAdapter()
  95               .formatSequences(FileFormat.Clustal, exp, true));
  96       if (flnk == -1)
  97       {
  98         /*
  99          * Full expansion to complete sequences
 100          */
 101         for (SequenceI sq : exp.getSequences())
 102         {
 103           String ung = sq.getSequenceAsString().replaceAll("-+", "");
 104           final String errorMsg = "Flanking sequence not the same as original dataset sequence.\n"
 105                   + ung + "\n"
 106                   + sq.getDatasetSequence().getSequenceAsString();
 107           assertTrue(errorMsg, ung.equalsIgnoreCase(
 108                   sq.getDatasetSequence().getSequenceAsString()));
 109         }
 110       }
 111       else if (flnk == 24)
 112       {
 113         /*
 114          * Last sequence is fully expanded, others have leading gaps to match
 115          */
 116         assertTrue(exp.getSequenceAt(4).getSequenceAsString()
 117                 .startsWith("abc"));
 118         assertTrue(exp.getSequenceAt(3).getSequenceAsString()
 119                 .startsWith("--abc"));
 120         assertTrue(exp.getSequenceAt(2).getSequenceAsString()
 121                 .startsWith("----abc"));
 122         assertTrue(exp.getSequenceAt(1).getSequenceAsString()
 123                 .startsWith("------abc"));
 124         assertTrue(exp.getSequenceAt(0).getSequenceAsString()
 125                 .startsWith("--------abc"));
 126       }
 127     }
 128   }
 129
 130   /**
 131    * Test that annotations are correctly adjusted by expandContext
 132    */
 133   @Test(groups = { "Functional" })
 134   public void testExpandContext_annotation()
 135   {
 136     AlignmentI al = new Alignment(new Sequence[] {});
 137     SequenceI ds = new Sequence("Seq1", "ABCDEFGHI");
 138     // subsequence DEF:
 139     SequenceI seq1 = ds.deriveSequence().getSubSequence(3, 6);
 140     al.addSequence(seq1);
 141
 142     /*
 143      * Annotate DEF with 4/5/6 respectively
 144      */
 145     Annotation[] anns = new Annotation[] { new Annotation(4),
 146         new Annotation(5), new Annotation(6) };
 147     AlignmentAnnotation ann = new AlignmentAnnotation("SS",
 148             "secondary structure", anns);
 149     seq1.addAlignmentAnnotation(ann);
 150
 151     /*
 152      * The annotations array should match aligned positions
 153      */
 154     assertEquals(3, ann.annotations.length);
 155     assertEquals(4, ann.annotations[0].value, 0.001);
 156     assertEquals(5, ann.annotations[1].value, 0.001);
 157     assertEquals(6, ann.annotations[2].value, 0.001);
 158
 159     /*
 160      * Check annotation to sequence position mappings before expanding the
 161      * sequence; these are set up in Sequence.addAlignmentAnnotation ->
 162      * Annotation.setSequenceRef -> createSequenceMappings
 163      */
 164     assertNull(ann.getAnnotationForPosition(1));
 165     assertNull(ann.getAnnotationForPosition(2));
 166     assertNull(ann.getAnnotationForPosition(3));
 167     assertEquals(4, ann.getAnnotationForPosition(4).value, 0.001);
 168     assertEquals(5, ann.getAnnotationForPosition(5).value, 0.001);
 169     assertEquals(6, ann.getAnnotationForPosition(6).value, 0.001);
 170     assertNull(ann.getAnnotationForPosition(7));
 171     assertNull(ann.getAnnotationForPosition(8));
 172     assertNull(ann.getAnnotationForPosition(9));
 173
 174     /*
 175      * Expand the subsequence to the full sequence abcDEFghi
 176      */
 177     AlignmentI expanded = AlignmentUtils.expandContext(al, -1);
 178     assertEquals("abcDEFghi",
 179             expanded.getSequenceAt(0).getSequenceAsString());
 180
 181     /*
 182      * Confirm the alignment and sequence have the same SS annotation,
 183      * referencing the expanded sequence
 184      */
 185     ann = expanded.getSequenceAt(0).getAnnotation()[0];
 186     assertSame(ann, expanded.getAlignmentAnnotation()[0]);
 187     assertSame(expanded.getSequenceAt(0), ann.sequenceRef);
 188
 189     /*
 190      * The annotations array should have null values except for annotated
 191      * positions
 192      */
 193     assertNull(ann.annotations[0]);
 194     assertNull(ann.annotations[1]);
 195     assertNull(ann.annotations[2]);
 196     assertEquals(4, ann.annotations[3].value, 0.001);
 197     assertEquals(5, ann.annotations[4].value, 0.001);
 198     assertEquals(6, ann.annotations[5].value, 0.001);
 199     assertNull(ann.annotations[6]);
 200     assertNull(ann.annotations[7]);
 201     assertNull(ann.annotations[8]);
 202
 203     /*
 204      * sequence position mappings should be unchanged
 205      */
 206     assertNull(ann.getAnnotationForPosition(1));
 207     assertNull(ann.getAnnotationForPosition(2));
 208     assertNull(ann.getAnnotationForPosition(3));
 209     assertEquals(4, ann.getAnnotationForPosition(4).value, 0.001);
 210     assertEquals(5, ann.getAnnotationForPosition(5).value, 0.001);
 211     assertEquals(6, ann.getAnnotationForPosition(6).value, 0.001);
 212     assertNull(ann.getAnnotationForPosition(7));
 213     assertNull(ann.getAnnotationForPosition(8));
 214     assertNull(ann.getAnnotationForPosition(9));
 215   }
 216
 217   /**
 218    * Test method that returns a map of lists of sequences by sequence name.
 219    *
 220    * @throws IOException
 221    */
 222   @Test(groups = { "Functional" })
 223   public void testGetSequencesByName() throws IOException
 224   {
 225     final String data = ">Seq1Name\nKQYL\n" + ">Seq2Name\nRFPW\n"
 226             + ">Seq1Name\nABCD\n";
 227     AlignmentI al = loadAlignment(data, FileFormat.Fasta);
 228     Map<String, List<SequenceI>> map = AlignmentUtils
 229             .getSequencesByName(al);
 230     assertEquals(2, map.keySet().size());
 231     assertEquals(2, map.get("Seq1Name").size());
 232     assertEquals("KQYL", map.get("Seq1Name").get(0).getSequenceAsString());
 233     assertEquals("ABCD", map.get("Seq1Name").get(1).getSequenceAsString());
 234     assertEquals(1, map.get("Seq2Name").size());
 235     assertEquals("RFPW", map.get("Seq2Name").get(0).getSequenceAsString());
 236   }
 237
 238   /**
 239    * Helper method to load an alignment and ensure dataset sequences are set up.
 240    *
 241    * @param data
 242    * @param format
 243    *          TODO
 244    * @return
 245    * @throws IOException
 246    */
 247   protected AlignmentI loadAlignment(final String data, FileFormatI format)
 248           throws IOException
 249   {
 250     AlignmentI a = new FormatAdapter().readFile(data, DataSourceType.PASTE,
 251             format);
 252     a.setDataset(null);
 253     return a;
 254   }
 255
 256   /**
 257    * Test mapping of protein to cDNA, for the case where we have no sequence
 258    * cross-references, so mappings are made first-served 1-1 where sequences
 259    * translate.
 260    *
 261    * @throws IOException
 262    */
 263   @Test(groups = { "Functional" })
 264   public void testMapProteinAlignmentToCdna_noXrefs() throws IOException
 265   {
 266     List<SequenceI> protseqs = new ArrayList<>();
 267     protseqs.add(new Sequence("UNIPROT|V12345", "EIQ"));
 268     protseqs.add(new Sequence("UNIPROT|V12346", "EIQ"));
 269     protseqs.add(new Sequence("UNIPROT|V12347", "SAR"));
 270     AlignmentI protein = new Alignment(protseqs.toArray(new SequenceI[3]));
 271     protein.setDataset(null);
 272
 273     List<SequenceI> dnaseqs = new ArrayList<>();
 274     dnaseqs.add(new Sequence("EMBL|A11111", "TCAGCACGC")); // = SAR
 275     dnaseqs.add(new Sequence("EMBL|A22222", "GAGATACAA")); // = EIQ
 276     dnaseqs.add(new Sequence("EMBL|A33333", "GAAATCCAG")); // = EIQ
 277     dnaseqs.add(new Sequence("EMBL|A44444", "GAAATTCAG")); // = EIQ
 278     AlignmentI cdna = new Alignment(dnaseqs.toArray(new SequenceI[4]));
 279     cdna.setDataset(null);
 280
 281     assertTrue(AlignmentUtils.mapProteinAlignmentToCdna(protein, cdna));
 282
 283     // 3 mappings made, each from 1 to 1 sequence
 284     assertEquals(3, protein.getCodonFrames().size());
 285     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(0)).size());
 286     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(1)).size());
 287     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(2)).size());
 288
 289     // V12345 mapped to A22222
 290     AlignedCodonFrame acf = protein.getCodonFrame(protein.getSequenceAt(0))
 291             .get(0);
 292     assertEquals(1, acf.getdnaSeqs().length);
 293     assertEquals(cdna.getSequenceAt(1).getDatasetSequence(),
 294             acf.getdnaSeqs()[0]);
 295     Mapping[] protMappings = acf.getProtMappings();
 296     assertEquals(1, protMappings.length);
 297     MapList mapList = protMappings[0].getMap();
 298     assertEquals(3, mapList.getFromRatio());
 299     assertEquals(1, mapList.getToRatio());
 300     assertTrue(
 301             Arrays.equals(new int[]
 302             { 1, 9 }, mapList.getFromRanges().get(0)));
 303     assertEquals(1, mapList.getFromRanges().size());
 304     assertTrue(
 305             Arrays.equals(new int[]
 306             { 1, 3 }, mapList.getToRanges().get(0)));
 307     assertEquals(1, mapList.getToRanges().size());
 308
 309     // V12346 mapped to A33333
 310     acf = protein.getCodonFrame(protein.getSequenceAt(1)).get(0);
 311     assertEquals(1, acf.getdnaSeqs().length);
 312     assertEquals(cdna.getSequenceAt(2).getDatasetSequence(),
 313             acf.getdnaSeqs()[0]);
 314
 315     // V12347 mapped to A11111
 316     acf = protein.getCodonFrame(protein.getSequenceAt(2)).get(0);
 317     assertEquals(1, acf.getdnaSeqs().length);
 318     assertEquals(cdna.getSequenceAt(0).getDatasetSequence(),
 319             acf.getdnaSeqs()[0]);
 320
 321     // no mapping involving the 'extra' A44444
 322     assertTrue(protein.getCodonFrame(cdna.getSequenceAt(3)).isEmpty());
 323   }
 324
 325   /**
 326    * Test for the alignSequenceAs method that takes two sequences and a mapping.
 327    */
 328   @Test(groups = { "Functional" })
 329   public void testAlignSequenceAs_withMapping_noIntrons()
 330   {
 331     MapList map = new MapList(new int[] { 1, 6 }, new int[] { 1, 2 }, 3, 1);
 332
 333     /*
 334      * No existing gaps in dna:
 335      */
 336     checkAlignSequenceAs("GGGAAA", "-A-L-", false, false, map,
 337             "---GGG---AAA");
 338
 339     /*
 340      * Now introduce gaps in dna but ignore them when realigning.
 341      */
 342     checkAlignSequenceAs("-G-G-G-A-A-A-", "-A-L-", false, false, map,
 343             "---GGG---AAA");
 344
 345     /*
 346      * Now include gaps in dna when realigning. First retaining 'mapped' gaps
 347      * only, i.e. those within the exon region.
 348      */
 349     checkAlignSequenceAs("-G-G--G-A--A-A-", "-A-L-", true, false, map,
 350             "---G-G--G---A--A-A");
 351
 352     /*
 353      * Include all gaps in dna when realigning (within and without the exon
 354      * region). The leading gap, and the gaps between codons, are subsumed by
 355      * the protein alignment gap.
 356      */
 357     checkAlignSequenceAs("-G-GG--AA-A---", "-A-L-", true, true, map,
 358             "---G-GG---AA-A---");
 359
 360     /*
 361      * Include only unmapped gaps in dna when realigning (outside the exon
 362      * region). The leading gap, and the gaps between codons, are subsumed by
 363      * the protein alignment gap.
 364      */
 365     checkAlignSequenceAs("-G-GG--AA-A-", "-A-L-", false, true, map,
 366             "---GGG---AAA---");
 367   }
 368
 369   /**
 370    * Test for the alignSequenceAs method that takes two sequences and a mapping.
 371    */
 372   @Test(groups = { "Functional" })
 373   public void testAlignSequenceAs_withMapping_withIntrons()
 374   {
 375     /*
 376      * Exons at codon 2 (AAA) and 4 (TTT)
 377      */
 378     MapList map = new MapList(new int[] { 4, 6, 10, 12 },
 379             new int[]
 380             { 1, 2 }, 3, 1);
 381
 382     /*
 383      * Simple case: no gaps in dna
 384      */
 385     checkAlignSequenceAs("GGGAAACCCTTTGGG", "--A-L-", false, false, map,
 386             "GGG---AAACCCTTTGGG");
 387
 388     /*
 389      * Add gaps to dna - but ignore when realigning.
 390      */
 391     checkAlignSequenceAs("-G-G-G--A--A---AC-CC-T-TT-GG-G-", "--A-L-", false,
 392             false, map, "GGG---AAACCCTTTGGG");
 393
 394     /*
 395      * Add gaps to dna - include within exons only when realigning.
 396      */
 397     checkAlignSequenceAs("-G-G-G--A--A---A-C-CC-T-TT-GG-G-", "--A-L-", true,
 398             false, map, "GGG---A--A---ACCCT-TTGGG");
 399
 400     /*
 401      * Include gaps outside exons only when realigning.
 402      */
 403     checkAlignSequenceAs("-G-G-G--A--A---A-C-CC-T-TT-GG-G-", "--A-L-",
 404             false, true, map, "-G-G-GAAAC-CCTTT-GG-G-");
 405
 406     /*
 407      * Include gaps following first intron if we are 'preserving mapped gaps'
 408      */
 409     checkAlignSequenceAs("-G-G-G--A--A---A-C-CC-T-TT-GG-G-", "--A-L-", true,
 410             true, map, "-G-G-G--A--A---A-C-CC-T-TT-GG-G-");
 411
 412     /*
 413      * Include all gaps in dna when realigning.
 414      */
 415     checkAlignSequenceAs("-G-G-G--A--A---A-C-CC-T-TT-GG-G-", "--A-L-", true,
 416             true, map, "-G-G-G--A--A---A-C-CC-T-TT-GG-G-");
 417   }
 418
 419   /**
 420    * Test for the case where not all of the protein sequence is mapped to cDNA.
 421    */
 422   @Test(groups = { "Functional" })
 423   public void testAlignSequenceAs_withMapping_withUnmappedProtein()
 424   {
 425     /*
 426      * Exons at codon 2 (AAA) and 4 (TTT) mapped to A and P
 427      */
 428     final MapList map = new MapList(new int[] { 4, 6, 10, 12 },
 429             new int[]
 430             { 1, 1, 3, 3 }, 3, 1);
 431
 432     /*
 433      * -L- 'aligns' ccc------
 434      */
 435     checkAlignSequenceAs("gggAAAcccTTTggg", "-A-L-P-", false, false, map,
 436             "gggAAAccc------TTTggg");
 437   }
 438
 439   /**
 440    * Helper method that performs and verifies the method under test.
 441    *
 442    * @param alignee
 443    *          the sequence to be realigned
 444    * @param alignModel
 445    *          the sequence whose alignment is to be copied
 446    * @param preserveMappedGaps
 447    * @param preserveUnmappedGaps
 448    * @param map
 449    * @param expected
 450    */
 451   protected void checkAlignSequenceAs(final String alignee,
 452           final String alignModel, final boolean preserveMappedGaps,
 453           final boolean preserveUnmappedGaps, MapList map,
 454           final String expected)
 455   {
 456     SequenceI alignMe = new Sequence("Seq1", alignee);
 457     alignMe.createDatasetSequence();
 458     SequenceI alignFrom = new Sequence("Seq2", alignModel);
 459     alignFrom.createDatasetSequence();
 460     AlignedCodonFrame acf = new AlignedCodonFrame();
 461     acf.addMap(alignMe.getDatasetSequence(), alignFrom.getDatasetSequence(),
 462             map);
 463
 464     AlignmentUtils.alignSequenceAs(alignMe, alignFrom, acf, "---", '-',
 465             preserveMappedGaps, preserveUnmappedGaps);
 466     assertEquals(expected, alignMe.getSequenceAsString());
 467   }
 468
 469   /**
 470    * Test for the alignSequenceAs method where we preserve gaps in introns only.
 471    */
 472   @Test(groups = { "Functional" })
 473   public void testAlignSequenceAs_keepIntronGapsOnly()
 474   {
 475
 476     /*
 477      * Intron GGGAAA followed by exon CCCTTT
 478      */
 479     MapList map = new MapList(new int[] { 7, 12 }, new int[] { 1, 2 }, 3,
 480             1);
 481
 482     checkAlignSequenceAs("GG-G-AA-A-C-CC-T-TT", "AL", false, true, map,
 483             "GG-G-AA-ACCCTTT");
 484   }
 485
 486   /**
 487    * Test the method that realigns protein to match mapped codon alignment.
 488    */
 489   @Test(groups = { "Functional" })
 490   public void testAlignProteinAsDna()
 491   {
 492     // seq1 codons are [1,2,3] [4,5,6] [7,8,9] [10,11,12]
 493     SequenceI dna1 = new Sequence("Seq1", "TGCCATTACCAG-");
 494     // seq2 codons are [1,3,4] [5,6,7] [8,9,10] [11,12,13]
 495     SequenceI dna2 = new Sequence("Seq2", "T-GCCATTACCAG");
 496     // seq3 codons are [1,2,3] [4,5,7] [8,9,10] [11,12,13]
 497     SequenceI dna3 = new Sequence("Seq3", "TGCCA-TTACCAG");
 498     AlignmentI dna = new Alignment(new SequenceI[] { dna1, dna2, dna3 });
 499     dna.setDataset(null);
 500
 501     // protein alignment will be realigned like dna
 502     SequenceI prot1 = new Sequence("Seq1", "CHYQ");
 503     SequenceI prot2 = new Sequence("Seq2", "CHYQ");
 504     SequenceI prot3 = new Sequence("Seq3", "CHYQ");
 505     SequenceI prot4 = new Sequence("Seq4", "R-QSV"); // unmapped, unchanged
 506     AlignmentI protein = new Alignment(
 507             new SequenceI[]
 508             { prot1, prot2, prot3, prot4 });
 509     protein.setDataset(null);
 510
 511     MapList map = new MapList(new int[] { 1, 12 }, new int[] { 1, 4 }, 3,
 512             1);
 513     AlignedCodonFrame acf = new AlignedCodonFrame();
 514     acf.addMap(dna1.getDatasetSequence(), prot1.getDatasetSequence(), map);
 515     acf.addMap(dna2.getDatasetSequence(), prot2.getDatasetSequence(), map);
 516     acf.addMap(dna3.getDatasetSequence(), prot3.getDatasetSequence(), map);
 517     ArrayList<AlignedCodonFrame> acfs = new ArrayList<>();
 518     acfs.add(acf);
 519     protein.setCodonFrames(acfs);
 520
 521     /*
 522      * Translated codon order is [1,2,3] [1,3,4] [4,5,6] [4,5,7] [5,6,7] [7,8,9]
 523      * [8,9,10] [10,11,12] [11,12,13]
 524      */
 525     AlignmentUtils.alignProteinAsDna(protein, dna);
 526     assertEquals("C-H--Y-Q-", prot1.getSequenceAsString());
 527     assertEquals("-C--H-Y-Q", prot2.getSequenceAsString());
 528     assertEquals("C--H--Y-Q", prot3.getSequenceAsString());
 529     assertEquals("R-QSV", prot4.getSequenceAsString());
 530   }
 531
 532   /**
 533    * Test the method that tests whether a CDNA sequence translates to a protein
 534    * sequence
 535    */
 536   @Test(groups = { "Functional" })
 537   public void testTranslatesAs()
 538   {
 539     // null arguments check
 540     assertFalse(AlignmentUtils.translatesAs(null, 0, null));
 541     assertFalse(AlignmentUtils.translatesAs(new char[] { 't' }, 0, null));
 542     assertFalse(AlignmentUtils.translatesAs(null, 0, new char[] { 'a' }));
 543
 544     // straight translation
 545     assertTrue(AlignmentUtils.translatesAs("tttcccaaaggg".toCharArray(), 0,
 546             "FPKG".toCharArray()));
 547     // with extra start codon (not in protein)
 548     assertTrue(AlignmentUtils.translatesAs("atgtttcccaaaggg".toCharArray(),
 549             3, "FPKG".toCharArray()));
 550     // with stop codon1 (not in protein)
 551     assertTrue(AlignmentUtils.translatesAs("tttcccaaagggtaa".toCharArray(),
 552             0, "FPKG".toCharArray()));
 553     // with stop codon1 (in protein as *)
 554     assertTrue(AlignmentUtils.translatesAs("tttcccaaagggtaa".toCharArray(),
 555             0, "FPKG*".toCharArray()));
 556     // with stop codon2 (not in protein)
 557     assertTrue(AlignmentUtils.translatesAs("tttcccaaagggtag".toCharArray(),
 558             0, "FPKG".toCharArray()));
 559     // with stop codon3 (not in protein)
 560     assertTrue(AlignmentUtils.translatesAs("tttcccaaagggtga".toCharArray(),
 561             0, "FPKG".toCharArray()));
 562     // with start and stop codon1
 563     assertTrue(AlignmentUtils.translatesAs(
 564             "atgtttcccaaagggtaa".toCharArray(), 3, "FPKG".toCharArray()));
 565     // with start and stop codon1 (in protein as *)
 566     assertTrue(AlignmentUtils.translatesAs(
 567             "atgtttcccaaagggtaa".toCharArray(), 3, "FPKG*".toCharArray()));
 568     // with start and stop codon2
 569     assertTrue(AlignmentUtils.translatesAs(
 570             "atgtttcccaaagggtag".toCharArray(), 3, "FPKG".toCharArray()));
 571     // with start and stop codon3
 572     assertTrue(AlignmentUtils.translatesAs(
 573             "atgtttcccaaagggtga".toCharArray(), 3, "FPKG".toCharArray()));
 574
 575     // with embedded stop codons
 576     assertTrue(AlignmentUtils.translatesAs(
 577             "atgtttTAGcccaaaTAAgggtga".toCharArray(), 3,
 578             "F*PK*G".toCharArray()));
 579
 580     // wrong protein
 581     assertFalse(AlignmentUtils.translatesAs("tttcccaaaggg".toCharArray(), 0,
 582             "FPMG".toCharArray()));
 583
 584     // truncated dna
 585     assertFalse(AlignmentUtils.translatesAs("tttcccaaagg".toCharArray(), 0,
 586             "FPKG".toCharArray()));
 587
 588     // truncated protein
 589     assertFalse(AlignmentUtils.translatesAs("tttcccaaaggg".toCharArray(), 0,
 590             "FPK".toCharArray()));
 591
 592     // overlong dna (doesn't end in stop codon)
 593     assertFalse(AlignmentUtils.translatesAs("tttcccaaagggttt".toCharArray(),
 594             0, "FPKG".toCharArray()));
 595
 596     // dna + stop codon + more
 597     assertFalse(AlignmentUtils.translatesAs(
 598             "tttcccaaagggttaga".toCharArray(), 0, "FPKG".toCharArray()));
 599
 600     // overlong protein
 601     assertFalse(AlignmentUtils.translatesAs("tttcccaaaggg".toCharArray(), 0,
 602             "FPKGQ".toCharArray()));
 603   }
 604
 605   /**
 606    * Test mapping of protein to cDNA, for cases where the cDNA has start and/or
 607    * stop codons in addition to the protein coding sequence.
 608    *
 609    * @throws IOException
 610    */
 611   @Test(groups = { "Functional" })
 612   public void testMapProteinAlignmentToCdna_withStartAndStopCodons()
 613           throws IOException
 614   {
 615     List<SequenceI> protseqs = new ArrayList<>();
 616     protseqs.add(new Sequence("UNIPROT|V12345", "EIQ"));
 617     protseqs.add(new Sequence("UNIPROT|V12346", "EIQ"));
 618     protseqs.add(new Sequence("UNIPROT|V12347", "SAR"));
 619     AlignmentI protein = new Alignment(protseqs.toArray(new SequenceI[3]));
 620     protein.setDataset(null);
 621
 622     List<SequenceI> dnaseqs = new ArrayList<>();
 623     // start + SAR:
 624     dnaseqs.add(new Sequence("EMBL|A11111", "ATGTCAGCACGC"));
 625     // = EIQ + stop
 626     dnaseqs.add(new Sequence("EMBL|A22222", "GAGATACAATAA"));
 627     // = start +EIQ + stop
 628     dnaseqs.add(new Sequence("EMBL|A33333", "ATGGAAATCCAGTAG"));
 629     dnaseqs.add(new Sequence("EMBL|A44444", "GAAATTCAG"));
 630     AlignmentI cdna = new Alignment(dnaseqs.toArray(new SequenceI[4]));
 631     cdna.setDataset(null);
 632
 633     assertTrue(AlignmentUtils.mapProteinAlignmentToCdna(protein, cdna));
 634
 635     // 3 mappings made, each from 1 to 1 sequence
 636     assertEquals(3, protein.getCodonFrames().size());
 637     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(0)).size());
 638     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(1)).size());
 639     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(2)).size());
 640
 641     // V12345 mapped from A22222
 642     AlignedCodonFrame acf = protein.getCodonFrame(protein.getSequenceAt(0))
 643             .get(0);
 644     assertEquals(1, acf.getdnaSeqs().length);
 645     assertEquals(cdna.getSequenceAt(1).getDatasetSequence(),
 646             acf.getdnaSeqs()[0]);
 647     Mapping[] protMappings = acf.getProtMappings();
 648     assertEquals(1, protMappings.length);
 649     MapList mapList = protMappings[0].getMap();
 650     assertEquals(3, mapList.getFromRatio());
 651     assertEquals(1, mapList.getToRatio());
 652     assertTrue(
 653             Arrays.equals(new int[]
 654             { 1, 9 }, mapList.getFromRanges().get(0)));
 655     assertEquals(1, mapList.getFromRanges().size());
 656     assertTrue(
 657             Arrays.equals(new int[]
 658             { 1, 3 }, mapList.getToRanges().get(0)));
 659     assertEquals(1, mapList.getToRanges().size());
 660
 661     // V12346 mapped from A33333 starting position 4
 662     acf = protein.getCodonFrame(protein.getSequenceAt(1)).get(0);
 663     assertEquals(1, acf.getdnaSeqs().length);
 664     assertEquals(cdna.getSequenceAt(2).getDatasetSequence(),
 665             acf.getdnaSeqs()[0]);
 666     protMappings = acf.getProtMappings();
 667     assertEquals(1, protMappings.length);
 668     mapList = protMappings[0].getMap();
 669     assertEquals(3, mapList.getFromRatio());
 670     assertEquals(1, mapList.getToRatio());
 671     assertTrue(
 672             Arrays.equals(new int[]
 673             { 4, 12 }, mapList.getFromRanges().get(0)));
 674     assertEquals(1, mapList.getFromRanges().size());
 675     assertTrue(
 676             Arrays.equals(new int[]
 677             { 1, 3 }, mapList.getToRanges().get(0)));
 678     assertEquals(1, mapList.getToRanges().size());
 679
 680     // V12347 mapped to A11111 starting position 4
 681     acf = protein.getCodonFrame(protein.getSequenceAt(2)).get(0);
 682     assertEquals(1, acf.getdnaSeqs().length);
 683     assertEquals(cdna.getSequenceAt(0).getDatasetSequence(),
 684             acf.getdnaSeqs()[0]);
 685     protMappings = acf.getProtMappings();
 686     assertEquals(1, protMappings.length);
 687     mapList = protMappings[0].getMap();
 688     assertEquals(3, mapList.getFromRatio());
 689     assertEquals(1, mapList.getToRatio());
 690     assertTrue(
 691             Arrays.equals(new int[]
 692             { 4, 12 }, mapList.getFromRanges().get(0)));
 693     assertEquals(1, mapList.getFromRanges().size());
 694     assertTrue(
 695             Arrays.equals(new int[]
 696             { 1, 3 }, mapList.getToRanges().get(0)));
 697     assertEquals(1, mapList.getToRanges().size());
 698
 699     // no mapping involving the 'extra' A44444
 700     assertTrue(protein.getCodonFrame(cdna.getSequenceAt(3)).isEmpty());
 701   }
 702
 703   /**
 704    * Test mapping of protein to cDNA, for the case where we have some sequence
 705    * cross-references. Verify that 1-to-many mappings are made where
 706    * cross-references exist and sequences are mappable.
 707    *
 708    * @throws IOException
 709    */
 710   @Test(groups = { "Functional" })
 711   public void testMapProteinAlignmentToCdna_withXrefs() throws IOException
 712   {
 713     List<SequenceI> protseqs = new ArrayList<>();
 714     protseqs.add(new Sequence("UNIPROT|V12345", "EIQ"));
 715     protseqs.add(new Sequence("UNIPROT|V12346", "EIQ"));
 716     protseqs.add(new Sequence("UNIPROT|V12347", "SAR"));
 717     AlignmentI protein = new Alignment(protseqs.toArray(new SequenceI[3]));
 718     protein.setDataset(null);
 719
 720     List<SequenceI> dnaseqs = new ArrayList<>();
 721     dnaseqs.add(new Sequence("EMBL|A11111", "TCAGCACGC")); // = SAR
 722     dnaseqs.add(new Sequence("EMBL|A22222", "ATGGAGATACAA")); // = start + EIQ
 723     dnaseqs.add(new Sequence("EMBL|A33333", "GAAATCCAG")); // = EIQ
 724     dnaseqs.add(new Sequence("EMBL|A44444", "GAAATTCAG")); // = EIQ
 725     dnaseqs.add(new Sequence("EMBL|A55555", "GAGATTCAG")); // = EIQ
 726     AlignmentI cdna = new Alignment(dnaseqs.toArray(new SequenceI[5]));
 727     cdna.setDataset(null);
 728
 729     // Xref A22222 to V12345 (should get mapped)
 730     dnaseqs.get(1).addDBRef(new DBRefEntry("UNIPROT", "1", "V12345"));
 731     // Xref V12345 to A44444 (should get mapped)
 732     protseqs.get(0).addDBRef(new DBRefEntry("EMBL", "1", "A44444"));
 733     // Xref A33333 to V12347 (sequence mismatch - should not get mapped)
 734     dnaseqs.get(2).addDBRef(new DBRefEntry("UNIPROT", "1", "V12347"));
 735     // as V12345 is mapped to A22222 and A44444, this leaves V12346 unmapped.
 736     // it should get paired up with the unmapped A33333
 737     // A11111 should be mapped to V12347
 738     // A55555 is spare and has no xref so is not mapped
 739
 740     assertTrue(AlignmentUtils.mapProteinAlignmentToCdna(protein, cdna));
 741
 742     // 4 protein mappings made for 3 proteins, 2 to V12345, 1 each to V12346/7
 743     assertEquals(3, protein.getCodonFrames().size());
 744     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(0)).size());
 745     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(1)).size());
 746     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(2)).size());
 747
 748     // one mapping for each of the first 4 cDNA sequences
 749     assertEquals(1, protein.getCodonFrame(cdna.getSequenceAt(0)).size());
 750     assertEquals(1, protein.getCodonFrame(cdna.getSequenceAt(1)).size());
 751     assertEquals(1, protein.getCodonFrame(cdna.getSequenceAt(2)).size());
 752     assertEquals(1, protein.getCodonFrame(cdna.getSequenceAt(3)).size());
 753
 754     // V12345 mapped to A22222 and A44444
 755     AlignedCodonFrame acf = protein.getCodonFrame(protein.getSequenceAt(0))
 756             .get(0);
 757     assertEquals(2, acf.getdnaSeqs().length);
 758     assertEquals(cdna.getSequenceAt(1).getDatasetSequence(),
 759             acf.getdnaSeqs()[0]);
 760     assertEquals(cdna.getSequenceAt(3).getDatasetSequence(),
 761             acf.getdnaSeqs()[1]);
 762
 763     // V12346 mapped to A33333
 764     acf = protein.getCodonFrame(protein.getSequenceAt(1)).get(0);
 765     assertEquals(1, acf.getdnaSeqs().length);
 766     assertEquals(cdna.getSequenceAt(2).getDatasetSequence(),
 767             acf.getdnaSeqs()[0]);
 768
 769     // V12347 mapped to A11111
 770     acf = protein.getCodonFrame(protein.getSequenceAt(2)).get(0);
 771     assertEquals(1, acf.getdnaSeqs().length);
 772     assertEquals(cdna.getSequenceAt(0).getDatasetSequence(),
 773             acf.getdnaSeqs()[0]);
 774
 775     // no mapping involving the 'extra' A55555
 776     assertTrue(protein.getCodonFrame(cdna.getSequenceAt(4)).isEmpty());
 777   }
 778
 779   /**
 780    * Test mapping of protein to cDNA, for the case where we have some sequence
 781    * cross-references. Verify that once we have made an xref mapping we don't
 782    * also map un-xrefd sequeces.
 783    *
 784    * @throws IOException
 785    */
 786   @Test(groups = { "Functional" })
 787   public void testMapProteinAlignmentToCdna_prioritiseXrefs()
 788           throws IOException
 789   {
 790     List<SequenceI> protseqs = new ArrayList<>();
 791     protseqs.add(new Sequence("UNIPROT|V12345", "EIQ"));
 792     protseqs.add(new Sequence("UNIPROT|V12346", "EIQ"));
 793     AlignmentI protein = new Alignment(
 794             protseqs.toArray(new SequenceI[protseqs.size()]));
 795     protein.setDataset(null);
 796
 797     List<SequenceI> dnaseqs = new ArrayList<>();
 798     dnaseqs.add(new Sequence("EMBL|A11111", "GAAATCCAG")); // = EIQ
 799     dnaseqs.add(new Sequence("EMBL|A22222", "GAAATTCAG")); // = EIQ
 800     AlignmentI cdna = new Alignment(
 801             dnaseqs.toArray(new SequenceI[dnaseqs.size()]));
 802     cdna.setDataset(null);
 803
 804     // Xref A22222 to V12345 (should get mapped)
 805     // A11111 should then be mapped to the unmapped V12346
 806     dnaseqs.get(1).addDBRef(new DBRefEntry("UNIPROT", "1", "V12345"));
 807
 808     assertTrue(AlignmentUtils.mapProteinAlignmentToCdna(protein, cdna));
 809
 810     // 2 protein mappings made
 811     assertEquals(2, protein.getCodonFrames().size());
 812     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(0)).size());
 813     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(1)).size());
 814
 815     // one mapping for each of the cDNA sequences
 816     assertEquals(1, protein.getCodonFrame(cdna.getSequenceAt(0)).size());
 817     assertEquals(1, protein.getCodonFrame(cdna.getSequenceAt(1)).size());
 818
 819     // V12345 mapped to A22222
 820     AlignedCodonFrame acf = protein.getCodonFrame(protein.getSequenceAt(0))
 821             .get(0);
 822     assertEquals(1, acf.getdnaSeqs().length);
 823     assertEquals(cdna.getSequenceAt(1).getDatasetSequence(),
 824             acf.getdnaSeqs()[0]);
 825
 826     // V12346 mapped to A11111
 827     acf = protein.getCodonFrame(protein.getSequenceAt(1)).get(0);
 828     assertEquals(1, acf.getdnaSeqs().length);
 829     assertEquals(cdna.getSequenceAt(0).getDatasetSequence(),
 830             acf.getdnaSeqs()[0]);
 831   }
 832
 833   /**
 834    * Test the method that shows or hides sequence annotations by type(s) and
 835    * selection group.
 836    */
 837   @Test(groups = { "Functional" })
 838   public void testShowOrHideSequenceAnnotations()
 839   {
 840     SequenceI seq1 = new Sequence("Seq1", "AAA");
 841     SequenceI seq2 = new Sequence("Seq2", "BBB");
 842     SequenceI seq3 = new Sequence("Seq3", "CCC");
 843     Annotation[] anns = new Annotation[] { new Annotation(2f) };
 844     AlignmentAnnotation ann1 = new AlignmentAnnotation("Structure", "ann1",
 845             anns);
 846     ann1.setSequenceRef(seq1);
 847     AlignmentAnnotation ann2 = new AlignmentAnnotation("Structure", "ann2",
 848             anns);
 849     ann2.setSequenceRef(seq2);
 850     AlignmentAnnotation ann3 = new AlignmentAnnotation("Structure", "ann3",
 851             anns);
 852     AlignmentAnnotation ann4 = new AlignmentAnnotation("Temp", "ann4",
 853             anns);
 854     ann4.setSequenceRef(seq1);
 855     AlignmentAnnotation ann5 = new AlignmentAnnotation("Temp", "ann5",
 856             anns);
 857     ann5.setSequenceRef(seq2);
 858     AlignmentAnnotation ann6 = new AlignmentAnnotation("Temp", "ann6",
 859             anns);
 860     AlignmentI al = new Alignment(new SequenceI[] { seq1, seq2, seq3 });
 861     al.addAnnotation(ann1); // Structure for Seq1
 862     al.addAnnotation(ann2); // Structure for Seq2
 863     al.addAnnotation(ann3); // Structure for no sequence
 864     al.addAnnotation(ann4); // Temp for seq1
 865     al.addAnnotation(ann5); // Temp for seq2
 866     al.addAnnotation(ann6); // Temp for no sequence
 867     List<String> types = new ArrayList<>();
 868     List<SequenceI> scope = new ArrayList<>();
 869
 870     /*
 871      * Set all sequence related Structure to hidden (ann1, ann2)
 872      */
 873     types.add("Structure");
 874     AlignmentUtils.showOrHideSequenceAnnotations(al, types, null, false,
 875             false);
 876     assertFalse(ann1.visible);
 877     assertFalse(ann2.visible);
 878     assertTrue(ann3.visible); // not sequence-related, not affected
 879     assertTrue(ann4.visible); // not Structure, not affected
 880     assertTrue(ann5.visible); // "
 881     assertTrue(ann6.visible); // not sequence-related, not affected
 882
 883     /*
 884      * Set Temp in {seq1, seq3} to hidden
 885      */
 886     types.clear();
 887     types.add("Temp");
 888     scope.add(seq1);
 889     scope.add(seq3);
 890     AlignmentUtils.showOrHideSequenceAnnotations(al, types, scope, false,
 891             false);
 892     assertFalse(ann1.visible); // unchanged
 893     assertFalse(ann2.visible); // unchanged
 894     assertTrue(ann3.visible); // not sequence-related, not affected
 895     assertFalse(ann4.visible); // Temp for seq1 hidden
 896     assertTrue(ann5.visible); // not in scope, not affected
 897     assertTrue(ann6.visible); // not sequence-related, not affected
 898
 899     /*
 900      * Set Temp in all sequences to hidden
 901      */
 902     types.clear();
 903     types.add("Temp");
 904     scope.add(seq1);
 905     scope.add(seq3);
 906     AlignmentUtils.showOrHideSequenceAnnotations(al, types, null, false,
 907             false);
 908     assertFalse(ann1.visible); // unchanged
 909     assertFalse(ann2.visible); // unchanged
 910     assertTrue(ann3.visible); // not sequence-related, not affected
 911     assertFalse(ann4.visible); // Temp for seq1 hidden
 912     assertFalse(ann5.visible); // Temp for seq2 hidden
 913     assertTrue(ann6.visible); // not sequence-related, not affected
 914
 915     /*
 916      * Set all types in {seq1, seq3} to visible
 917      */
 918     types.clear();
 919     scope.clear();
 920     scope.add(seq1);
 921     scope.add(seq3);
 922     AlignmentUtils.showOrHideSequenceAnnotations(al, types, scope, true,
 923             true);
 924     assertTrue(ann1.visible); // Structure for seq1 set visible
 925     assertFalse(ann2.visible); // not in scope, unchanged
 926     assertTrue(ann3.visible); // not sequence-related, not affected
 927     assertTrue(ann4.visible); // Temp for seq1 set visible
 928     assertFalse(ann5.visible); // not in scope, unchanged
 929     assertTrue(ann6.visible); // not sequence-related, not affected
 930
 931     /*
 932      * Set all types in all scope to hidden
 933      */
 934     AlignmentUtils.showOrHideSequenceAnnotations(al, types, null, true,
 935             false);
 936     assertFalse(ann1.visible);
 937     assertFalse(ann2.visible);
 938     assertTrue(ann3.visible); // not sequence-related, not affected
 939     assertFalse(ann4.visible);
 940     assertFalse(ann5.visible);
 941     assertTrue(ann6.visible); // not sequence-related, not affected
 942   }
 943
 944   /**
 945    * Tests for the method that checks if one sequence cross-references another
 946    */
 947   @Test(groups = { "Functional" })
 948   public void testHasCrossRef()
 949   {
 950     assertFalse(AlignmentUtils.hasCrossRef(null, null));
 951     SequenceI seq1 = new Sequence("EMBL|A12345", "ABCDEF");
 952     assertFalse(AlignmentUtils.hasCrossRef(seq1, null));
 953     assertFalse(AlignmentUtils.hasCrossRef(null, seq1));
 954     SequenceI seq2 = new Sequence("UNIPROT|V20192", "ABCDEF");
 955     assertFalse(AlignmentUtils.hasCrossRef(seq1, seq2));
 956
 957     // different ref
 958     seq1.addDBRef(new DBRefEntry("UNIPROT", "1", "v20193"));
 959     assertFalse(AlignmentUtils.hasCrossRef(seq1, seq2));
 960
 961     // case-insensitive; version number is ignored
 962     seq1.addDBRef(new DBRefEntry("UNIPROT", "1", "v20192"));
 963     assertTrue(AlignmentUtils.hasCrossRef(seq1, seq2));
 964
 965     // right case!
 966     seq1.addDBRef(new DBRefEntry("UNIPROT", "1", "V20192"));
 967     assertTrue(AlignmentUtils.hasCrossRef(seq1, seq2));
 968     // test is one-way only
 969     assertFalse(AlignmentUtils.hasCrossRef(seq2, seq1));
 970   }
 971
 972   /**
 973    * Tests for the method that checks if either sequence cross-references the
 974    * other
 975    */
 976   @Test(groups = { "Functional" })
 977   public void testHaveCrossRef()
 978   {
 979     assertFalse(AlignmentUtils.hasCrossRef(null, null));
 980     SequenceI seq1 = new Sequence("EMBL|A12345", "ABCDEF");
 981     assertFalse(AlignmentUtils.haveCrossRef(seq1, null));
 982     assertFalse(AlignmentUtils.haveCrossRef(null, seq1));
 983     SequenceI seq2 = new Sequence("UNIPROT|V20192", "ABCDEF");
 984     assertFalse(AlignmentUtils.haveCrossRef(seq1, seq2));
 985
 986     seq1.addDBRef(new DBRefEntry("UNIPROT", "1", "V20192"));
 987     assertTrue(AlignmentUtils.haveCrossRef(seq1, seq2));
 988     // next is true for haveCrossRef, false for hasCrossRef
 989     assertTrue(AlignmentUtils.haveCrossRef(seq2, seq1));
 990
 991     // now the other way round
 992     seq1.setDBRefs(null);
 993     seq2.addDBRef(new DBRefEntry("EMBL", "1", "A12345"));
 994     assertTrue(AlignmentUtils.haveCrossRef(seq1, seq2));
 995     assertTrue(AlignmentUtils.haveCrossRef(seq2, seq1));
 996
 997     // now both ways
 998     seq1.addDBRef(new DBRefEntry("UNIPROT", "1", "V20192"));
 999     assertTrue(AlignmentUtils.haveCrossRef(seq1, seq2));
1000     assertTrue(AlignmentUtils.haveCrossRef(seq2, seq1));
1001   }
1002
1003   /**
1004    * Test the method that extracts the cds-only part of a dna alignment.
1005    */
1006   @Test(groups = { "Functional" })
1007   public void testMakeCdsAlignment()
1008   {
1009     /*
1010      * scenario:
1011      *     dna1 --> [4, 6] [10,12]        --> pep1
1012      *     dna2 --> [1, 3] [7, 9] [13,15] --> pep2
1013      */
1014     SequenceI dna1 = new Sequence("dna1", "aaaGGGcccTTTaaa");
1015     SequenceI dna2 = new Sequence("dna2", "GGGcccTTTaaaCCC");
1016     SequenceI pep1 = new Sequence("pep1", "GF");
1017     SequenceI pep2 = new Sequence("pep2", "GFP");
1018     pep1.addDBRef(new DBRefEntry("UNIPROT", "0", "pep1"));
1019     pep2.addDBRef(new DBRefEntry("UNIPROT", "0", "pep2"));
1020     dna1.createDatasetSequence();
1021     dna2.createDatasetSequence();
1022     pep1.createDatasetSequence();
1023     pep2.createDatasetSequence();
1024     AlignmentI dna = new Alignment(new SequenceI[] { dna1, dna2 });
1025     dna.setDataset(null);
1026
1027     /*
1028      * put a variant feature on dna2 base 8
1029      * - should transfer to cds2 base 5
1030      */
1031     dna2.addSequenceFeature(
1032             new SequenceFeature("variant", "hgmd", 8, 8, 0f, null));
1033
1034     /*
1035      * need a sourceDbRef if we are to construct dbrefs to the CDS
1036      * sequence from the dna contig sequences
1037      */
1038     DBRefEntry dbref = new DBRefEntry("ENSEMBL", "0", "dna1");
1039     dna1.getDatasetSequence().addDBRef(dbref);
1040     org.testng.Assert.assertEquals(dbref, dna1.getPrimaryDBRefs().get(0));
1041     dbref = new DBRefEntry("ENSEMBL", "0", "dna2");
1042     dna2.getDatasetSequence().addDBRef(dbref);
1043     org.testng.Assert.assertEquals(dbref, dna2.getPrimaryDBRefs().get(0));
1044
1045     /*
1046      * CDS sequences are 'discovered' from dna-to-protein mappings on the alignment
1047      * dataset (e.g. added from dbrefs by CrossRef.findXrefSequences)
1048      */
1049     MapList mapfordna1 = new MapList(new int[] { 4, 6, 10, 12 },
1050             new int[]
1051             { 1, 2 }, 3, 1);
1052     AlignedCodonFrame acf = new AlignedCodonFrame();
1053     acf.addMap(dna1.getDatasetSequence(), pep1.getDatasetSequence(),
1054             mapfordna1);
1055     dna.addCodonFrame(acf);
1056     MapList mapfordna2 = new MapList(new int[] { 1, 3, 7, 9, 13, 15 },
1057             new int[]
1058             { 1, 3 }, 3, 1);
1059     acf = new AlignedCodonFrame();
1060     acf.addMap(dna2.getDatasetSequence(), pep2.getDatasetSequence(),
1061             mapfordna2);
1062     dna.addCodonFrame(acf);
1063
1064     /*
1065      * In this case, mappings originally came from matching Uniprot accessions
1066      * - so need an xref on dna involving those regions.
1067      * These are normally constructed from CDS annotation
1068      */
1069     DBRefEntry dna1xref = new DBRefEntry("UNIPROT", "ENSEMBL", "pep1",
1070             new Mapping(mapfordna1));
1071     dna1.addDBRef(dna1xref);
1072     assertEquals(2, dna1.getDBRefs().size()); // to self and to pep1
1073     DBRefEntry dna2xref = new DBRefEntry("UNIPROT", "ENSEMBL", "pep2",
1074             new Mapping(mapfordna2));
1075     dna2.addDBRef(dna2xref);
1076     assertEquals(2, dna2.getDBRefs().size()); // to self and to pep2
1077
1078     /*
1079      * execute method under test:
1080      */
1081     AlignmentI cds = AlignmentUtils
1082             .makeCdsAlignment(new SequenceI[]
1083             { dna1, dna2 }, dna.getDataset(), null);
1084
1085     /*
1086      * verify cds sequences
1087      */
1088     assertEquals(2, cds.getSequences().size());
1089     assertEquals("GGGTTT", cds.getSequenceAt(0).getSequenceAsString());
1090     assertEquals("GGGTTTCCC", cds.getSequenceAt(1).getSequenceAsString());
1091
1092     /*
1093      * verify shared, extended alignment dataset
1094      */
1095     assertSame(dna.getDataset(), cds.getDataset());
1096     SequenceI cds1Dss = cds.getSequenceAt(0).getDatasetSequence();
1097     SequenceI cds2Dss = cds.getSequenceAt(1).getDatasetSequence();
1098     assertTrue(dna.getDataset().getSequences().contains(cds1Dss));
1099     assertTrue(dna.getDataset().getSequences().contains(cds2Dss));
1100
1101     /*
1102      * verify CDS has a dbref with mapping to peptide
1103      */
1104     assertNotNull(cds1Dss.getDBRefs());
1105     assertEquals(2, cds1Dss.getDBRefs().size());
1106     dbref = cds1Dss.getDBRefs().get(0);
1107     assertEquals(dna1xref.getSource(), dbref.getSource());
1108     // version is via ensembl's primary ref
1109     assertEquals(dna1xref.getVersion(), dbref.getVersion());
1110     assertEquals(dna1xref.getAccessionId(), dbref.getAccessionId());
1111     assertNotNull(dbref.getMap());
1112     assertSame(pep1.getDatasetSequence(), dbref.getMap().getTo());
1113     MapList cdsMapping = new MapList(new int[] { 1, 6 }, new int[] { 1, 2 },
1114             3, 1);
1115     assertEquals(cdsMapping, dbref.getMap().getMap());
1116
1117     /*
1118      * verify peptide has added a dbref with reverse mapping to CDS
1119      */
1120     assertNotNull(pep1.getDBRefs());
1121     // FIXME pep1.getDBRefs() is 1 - is that the correct behaviour ?
1122     assertEquals(2, pep1.getDBRefs().size());
1123     dbref = pep1.getDBRefs().get(1);
1124     assertEquals("ENSEMBL", dbref.getSource());
1125     assertEquals("0", dbref.getVersion());
1126     assertEquals("CDS|dna1", dbref.getAccessionId());
1127     assertNotNull(dbref.getMap());
1128     assertSame(cds1Dss, dbref.getMap().getTo());
1129     assertEquals(cdsMapping.getInverse(), dbref.getMap().getMap());
1130
1131     /*
1132      * verify cDNA has added a dbref with mapping to CDS
1133      */
1134     assertEquals(3, dna1.getDBRefs().size());
1135     DBRefEntry dbRefEntry = dna1.getDBRefs().get(2);
1136     assertSame(cds1Dss, dbRefEntry.getMap().getTo());
1137     MapList dnaToCdsMapping = new MapList(new int[] { 4, 6, 10, 12 },
1138             new int[]
1139             { 1, 6 }, 1, 1);
1140     assertEquals(dnaToCdsMapping, dbRefEntry.getMap().getMap());
1141     assertEquals(3, dna2.getDBRefs().size());
1142     dbRefEntry = dna2.getDBRefs().get(2);
1143     assertSame(cds2Dss, dbRefEntry.getMap().getTo());
1144     dnaToCdsMapping = new MapList(new int[] { 1, 3, 7, 9, 13, 15 },
1145             new int[]
1146             { 1, 9 }, 1, 1);
1147     assertEquals(dnaToCdsMapping, dbRefEntry.getMap().getMap());
1148
1149     /*
1150      * verify CDS has added a dbref with mapping to cDNA
1151      */
1152     assertEquals(2, cds1Dss.getDBRefs().size());
1153     dbRefEntry = cds1Dss.getDBRefs().get(1);
1154     assertSame(dna1.getDatasetSequence(), dbRefEntry.getMap().getTo());
1155     MapList cdsToDnaMapping = new MapList(new int[] { 1, 6 },
1156             new int[]
1157             { 4, 6, 10, 12 }, 1, 1);
1158     assertEquals(cdsToDnaMapping, dbRefEntry.getMap().getMap());
1159     assertEquals(2, cds2Dss.getDBRefs().size());
1160     dbRefEntry = cds2Dss.getDBRefs().get(1);
1161     assertSame(dna2.getDatasetSequence(), dbRefEntry.getMap().getTo());
1162     cdsToDnaMapping = new MapList(new int[] { 1, 9 },
1163             new int[]
1164             { 1, 3, 7, 9, 13, 15 }, 1, 1);
1165     assertEquals(cdsToDnaMapping, dbRefEntry.getMap().getMap());
1166
1167     /*
1168      * Verify mappings from CDS to peptide, cDNA to CDS, and cDNA to peptide
1169      * the mappings are on the shared alignment dataset
1170      * 6 mappings, 2*(DNA->CDS), 2*(DNA->Pep), 2*(CDS->Pep)
1171      */
1172     List<AlignedCodonFrame> cdsMappings = cds.getDataset().getCodonFrames();
1173     assertEquals(6, cdsMappings.size());
1174
1175     /*
1176      * verify that mapping sets for dna and cds alignments are different
1177      * [not current behaviour - all mappings are on the alignment dataset]
1178      */
1179     // select -> subselect type to test.
1180     // Assert.assertNotSame(dna.getCodonFrames(), cds.getCodonFrames());
1181     // assertEquals(4, dna.getCodonFrames().size());
1182     // assertEquals(4, cds.getCodonFrames().size());
1183
1184     /*
1185      * Two mappings involve pep1 (dna to pep1, cds to pep1)
1186      * Mapping from pep1 to GGGTTT in first new exon sequence
1187      */
1188     List<AlignedCodonFrame> pep1Mappings = MappingUtils
1189             .findMappingsForSequence(pep1, cdsMappings);
1190     assertEquals(2, pep1Mappings.size());
1191     List<AlignedCodonFrame> mappings = MappingUtils
1192             .findMappingsForSequence(cds.getSequenceAt(0), pep1Mappings);
1193     assertEquals(1, mappings.size());
1194
1195     // map G to GGG
1196     SearchResultsI sr = MappingUtils.buildSearchResults(pep1, 1, mappings);
1197     assertEquals(1, sr.getResults().size());
1198     SearchResultMatchI m = sr.getResults().get(0);
1199     assertSame(cds1Dss, m.getSequence());
1200     assertEquals(1, m.getStart());
1201     assertEquals(3, m.getEnd());
1202     // map F to TTT
1203     sr = MappingUtils.buildSearchResults(pep1, 2, mappings);
1204     m = sr.getResults().get(0);
1205     assertSame(cds1Dss, m.getSequence());
1206     assertEquals(4, m.getStart());
1207     assertEquals(6, m.getEnd());
1208
1209     /*
1210      * Two mappings involve pep2 (dna to pep2, cds to pep2)
1211      * Verify mapping from pep2 to GGGTTTCCC in second new exon sequence
1212      */
1213     List<AlignedCodonFrame> pep2Mappings = MappingUtils
1214             .findMappingsForSequence(pep2, cdsMappings);
1215     assertEquals(2, pep2Mappings.size());
1216     mappings = MappingUtils.findMappingsForSequence(cds.getSequenceAt(1),
1217             pep2Mappings);
1218     assertEquals(1, mappings.size());
1219     // map G to GGG
1220     sr = MappingUtils.buildSearchResults(pep2, 1, mappings);
1221     assertEquals(1, sr.getResults().size());
1222     m = sr.getResults().get(0);
1223     assertSame(cds2Dss, m.getSequence());
1224     assertEquals(1, m.getStart());
1225     assertEquals(3, m.getEnd());
1226     // map F to TTT
1227     sr = MappingUtils.buildSearchResults(pep2, 2, mappings);
1228     m = sr.getResults().get(0);
1229     assertSame(cds2Dss, m.getSequence());
1230     assertEquals(4, m.getStart());
1231     assertEquals(6, m.getEnd());
1232     // map P to CCC
1233     sr = MappingUtils.buildSearchResults(pep2, 3, mappings);
1234     m = sr.getResults().get(0);
1235     assertSame(cds2Dss, m.getSequence());
1236     assertEquals(7, m.getStart());
1237     assertEquals(9, m.getEnd());
1238
1239     /*
1240      * check cds2 acquired a variant feature in position 5
1241      */
1242     List<SequenceFeature> sfs = cds2Dss.getSequenceFeatures();
1243     assertNotNull(sfs);
1244     assertEquals(1, sfs.size());
1245     assertEquals("variant", sfs.get(0).type);
1246     assertEquals(5, sfs.get(0).begin);
1247     assertEquals(5, sfs.get(0).end);
1248   }
1249
1250   /**
1251    * Test the method that makes a cds-only alignment from a DNA sequence and its
1252    * product mappings, for the case where there are multiple exon mappings to
1253    * different protein products.
1254    */
1255   @Test(groups = { "Functional" })
1256   public void testMakeCdsAlignment_multipleProteins()
1257   {
1258     SequenceI dna1 = new Sequence("dna1", "aaaGGGcccTTTaaa");
1259     SequenceI pep1 = new Sequence("pep1", "GF"); // GGGTTT
1260     SequenceI pep2 = new Sequence("pep2", "KP"); // aaaccc
1261     SequenceI pep3 = new Sequence("pep3", "KF"); // aaaTTT
1262     dna1.createDatasetSequence();
1263     pep1.createDatasetSequence();
1264     pep2.createDatasetSequence();
1265     pep3.createDatasetSequence();
1266     pep1.getDatasetSequence()
1267             .addDBRef(new DBRefEntry("EMBLCDS", "2", "A12345"));
1268     pep2.getDatasetSequence()
1269             .addDBRef(new DBRefEntry("EMBLCDS", "3", "A12346"));
1270     pep3.getDatasetSequence()
1271             .addDBRef(new DBRefEntry("EMBLCDS", "4", "A12347"));
1272
1273     /*
1274      * Create the CDS alignment
1275      */
1276     AlignmentI dna = new Alignment(new SequenceI[] { dna1 });
1277     dna.setDataset(null);
1278
1279     /*
1280      * Make the mappings from dna to protein
1281      */
1282     // map ...GGG...TTT to GF
1283     MapList map = new MapList(new int[] { 4, 6, 10, 12 },
1284             new int[]
1285             { 1, 2 }, 3, 1);
1286     AlignedCodonFrame acf = new AlignedCodonFrame();
1287     acf.addMap(dna1.getDatasetSequence(), pep1.getDatasetSequence(), map);
1288     dna.addCodonFrame(acf);
1289
1290     // map aaa...ccc to KP
1291     map = new MapList(new int[] { 1, 3, 7, 9 }, new int[] { 1, 2 }, 3, 1);
1292     acf = new AlignedCodonFrame();
1293     acf.addMap(dna1.getDatasetSequence(), pep2.getDatasetSequence(), map);
1294     dna.addCodonFrame(acf);
1295
1296     // map aaa......TTT to KF
1297     map = new MapList(new int[] { 1, 3, 10, 12 }, new int[] { 1, 2 }, 3, 1);
1298     acf = new AlignedCodonFrame();
1299     acf.addMap(dna1.getDatasetSequence(), pep3.getDatasetSequence(), map);
1300     dna.addCodonFrame(acf);
1301
1302     /*
1303      * execute method under test
1304      */
1305     AlignmentI cdsal = AlignmentUtils
1306             .makeCdsAlignment(new SequenceI[]
1307             { dna1 }, dna.getDataset(), null);
1308
1309     /*
1310      * Verify we have 3 cds sequences, mapped to pep1/2/3 respectively
1311      */
1312     List<SequenceI> cds = cdsal.getSequences();
1313     assertEquals(3, cds.size());
1314
1315     /*
1316      * verify shared, extended alignment dataset
1317      */
1318     assertSame(cdsal.getDataset(), dna.getDataset());
1319     assertTrue(dna.getDataset().getSequences()
1320             .contains(cds.get(0).getDatasetSequence()));
1321     assertTrue(dna.getDataset().getSequences()
1322             .contains(cds.get(1).getDatasetSequence()));
1323     assertTrue(dna.getDataset().getSequences()
1324             .contains(cds.get(2).getDatasetSequence()));
1325
1326     /*
1327      * verify aligned cds sequences and their xrefs
1328      */
1329     SequenceI cdsSeq = cds.get(0);
1330     assertEquals("GGGTTT", cdsSeq.getSequenceAsString());
1331     // assertEquals("dna1|A12345", cdsSeq.getName());
1332     assertEquals("CDS|dna1", cdsSeq.getName());
1333     // assertEquals(1, cdsSeq.getDBRefs().length);
1334     // DBRefEntry cdsRef = cdsSeq.getDBRefs()[0];
1335     // assertEquals("EMBLCDS", cdsRef.getSource());
1336     // assertEquals("2", cdsRef.getVersion());
1337     // assertEquals("A12345", cdsRef.getAccessionId());
1338
1339     cdsSeq = cds.get(1);
1340     assertEquals("aaaccc", cdsSeq.getSequenceAsString());
1341     // assertEquals("dna1|A12346", cdsSeq.getName());
1342     assertEquals("CDS|dna1", cdsSeq.getName());
1343     // assertEquals(1, cdsSeq.getDBRefs().length);
1344     // cdsRef = cdsSeq.getDBRefs()[0];
1345     // assertEquals("EMBLCDS", cdsRef.getSource());
1346     // assertEquals("3", cdsRef.getVersion());
1347     // assertEquals("A12346", cdsRef.getAccessionId());
1348
1349     cdsSeq = cds.get(2);
1350     assertEquals("aaaTTT", cdsSeq.getSequenceAsString());
1351     // assertEquals("dna1|A12347", cdsSeq.getName());
1352     assertEquals("CDS|dna1", cdsSeq.getName());
1353     // assertEquals(1, cdsSeq.getDBRefs().length);
1354     // cdsRef = cdsSeq.getDBRefs()[0];
1355     // assertEquals("EMBLCDS", cdsRef.getSource());
1356     // assertEquals("4", cdsRef.getVersion());
1357     // assertEquals("A12347", cdsRef.getAccessionId());
1358
1359     /*
1360      * Verify there are mappings from each cds sequence to its protein product
1361      * and also to its dna source
1362      */
1363     List<AlignedCodonFrame> newMappings = cdsal.getCodonFrames();
1364
1365     /*
1366      * 6 mappings involve dna1 (to pep1/2/3, cds1/2/3)
1367      */
1368     List<AlignedCodonFrame> dnaMappings = MappingUtils
1369             .findMappingsForSequence(dna1, newMappings);
1370     assertEquals(6, dnaMappings.size());
1371
1372     /*
1373      * dna1 to pep1
1374      */
1375     List<AlignedCodonFrame> mappings = MappingUtils
1376             .findMappingsForSequence(pep1, dnaMappings);
1377     assertEquals(1, mappings.size());
1378     assertEquals(1, mappings.get(0).getMappings().size());
1379     assertSame(pep1.getDatasetSequence(),
1380             mappings.get(0).getMappings().get(0).getMapping().getTo());
1381
1382     /*
1383      * dna1 to cds1
1384      */
1385     List<AlignedCodonFrame> dnaToCds1Mappings = MappingUtils
1386             .findMappingsForSequence(cds.get(0), dnaMappings);
1387     Mapping mapping = dnaToCds1Mappings.get(0).getMappings().get(0)
1388             .getMapping();
1389     assertSame(cds.get(0).getDatasetSequence(), mapping.getTo());
1390     assertEquals("G(1) in CDS should map to G(4) in DNA", 4,
1391             mapping.getMap().getToPosition(1));
1392
1393     /*
1394      * dna1 to pep2
1395      */
1396     mappings = MappingUtils.findMappingsForSequence(pep2, dnaMappings);
1397     assertEquals(1, mappings.size());
1398     assertEquals(1, mappings.get(0).getMappings().size());
1399     assertSame(pep2.getDatasetSequence(),
1400             mappings.get(0).getMappings().get(0).getMapping().getTo());
1401
1402     /*
1403      * dna1 to cds2
1404      */
1405     List<AlignedCodonFrame> dnaToCds2Mappings = MappingUtils
1406             .findMappingsForSequence(cds.get(1), dnaMappings);
1407     mapping = dnaToCds2Mappings.get(0).getMappings().get(0).getMapping();
1408     assertSame(cds.get(1).getDatasetSequence(), mapping.getTo());
1409     assertEquals("c(4) in CDS should map to c(7) in DNA", 7,
1410             mapping.getMap().getToPosition(4));
1411
1412     /*
1413      * dna1 to pep3
1414      */
1415     mappings = MappingUtils.findMappingsForSequence(pep3, dnaMappings);
1416     assertEquals(1, mappings.size());
1417     assertEquals(1, mappings.get(0).getMappings().size());
1418     assertSame(pep3.getDatasetSequence(),
1419             mappings.get(0).getMappings().get(0).getMapping().getTo());
1420
1421     /*
1422      * dna1 to cds3
1423      */
1424     List<AlignedCodonFrame> dnaToCds3Mappings = MappingUtils
1425             .findMappingsForSequence(cds.get(2), dnaMappings);
1426     mapping = dnaToCds3Mappings.get(0).getMappings().get(0).getMapping();
1427     assertSame(cds.get(2).getDatasetSequence(), mapping.getTo());
1428     assertEquals("T(4) in CDS should map to T(10) in DNA", 10,
1429             mapping.getMap().getToPosition(4));
1430   }
1431
1432   @Test(groups = { "Functional" })
1433   public void testIsMappable()
1434   {
1435     SequenceI dna1 = new Sequence("dna1", "cgCAGtgGT");
1436     SequenceI aa1 = new Sequence("aa1", "RSG");
1437     AlignmentI al1 = new Alignment(new SequenceI[] { dna1 });
1438     AlignmentI al2 = new Alignment(new SequenceI[] { aa1 });
1439
1440     assertFalse(AlignmentUtils.isMappable(null, null));
1441     assertFalse(AlignmentUtils.isMappable(al1, null));
1442     assertFalse(AlignmentUtils.isMappable(null, al1));
1443     assertFalse(AlignmentUtils.isMappable(al1, al1));
1444     assertFalse(AlignmentUtils.isMappable(al2, al2));
1445
1446     assertTrue(AlignmentUtils.isMappable(al1, al2));
1447     assertTrue(AlignmentUtils.isMappable(al2, al1));
1448   }
1449
1450   /**
1451    * Test creating a mapping when the sequences involved do not start at residue
1452    * 1
1453    *
1454    * @throws IOException
1455    */
1456   @Test(groups = { "Functional" })
1457   public void testMapCdnaToProtein_forSubsequence() throws IOException
1458   {
1459     SequenceI prot = new Sequence("UNIPROT|V12345", "E-I--Q", 10, 12);
1460     prot.createDatasetSequence();
1461
1462     SequenceI dna = new Sequence("EMBL|A33333", "GAA--AT-C-CAG", 40, 48);
1463     dna.createDatasetSequence();
1464
1465     MapList map = AlignmentUtils.mapCdnaToProtein(prot, dna);
1466     assertEquals(10, map.getToLowest());
1467     assertEquals(12, map.getToHighest());
1468     assertEquals(40, map.getFromLowest());
1469     assertEquals(48, map.getFromHighest());
1470   }
1471
1472   /**
1473    * Test for the alignSequenceAs method where we have protein mapped to protein
1474    */
1475   @Test(groups = { "Functional" })
1476   public void testAlignSequenceAs_mappedProteinProtein()
1477   {
1478
1479     SequenceI alignMe = new Sequence("Match", "MGAASEV");
1480     alignMe.createDatasetSequence();
1481     SequenceI alignFrom = new Sequence("Query", "LQTGYMGAASEVMFSPTRR");
1482     alignFrom.createDatasetSequence();
1483
1484     AlignedCodonFrame acf = new AlignedCodonFrame();
1485     // this is like a domain or motif match of part of a peptide sequence
1486     MapList map = new MapList(new int[] { 6, 12 }, new int[] { 1, 7 }, 1,
1487             1);
1488     acf.addMap(alignFrom.getDatasetSequence(), alignMe.getDatasetSequence(),
1489             map);
1490
1491     AlignmentUtils.alignSequenceAs(alignMe, alignFrom, acf, "-", '-', true,
1492             true);
1493     assertEquals("-----MGAASEV-------", alignMe.getSequenceAsString());
1494   }
1495
1496   /**
1497    * Test for the alignSequenceAs method where there are trailing unmapped
1498    * residues in the model sequence
1499    */
1500   @Test(groups = { "Functional" })
1501   public void testAlignSequenceAs_withTrailingPeptide()
1502   {
1503     // map first 3 codons to KPF; G is a trailing unmapped residue
1504     MapList map = new MapList(new int[] { 1, 9 }, new int[] { 1, 3 }, 3, 1);
1505
1506     checkAlignSequenceAs("AAACCCTTT", "K-PFG", true, true, map,
1507             "AAA---CCCTTT---");
1508   }
1509
1510   /**
1511    * Tests for transferring features between mapped sequences
1512    */
1513   @Test(groups = { "Functional" })
1514   public void testTransferFeatures()
1515   {
1516     SequenceI dna = new Sequence("dna/20-34", "acgTAGcaaGCCcgt");
1517     SequenceI cds = new Sequence("cds/10-15", "TAGGCC");
1518
1519     // no overlap
1520     dna.addSequenceFeature(
1521             new SequenceFeature("type1", "desc1", 1, 2, 1f, null));
1522     // partial overlap - to [1, 1]
1523     dna.addSequenceFeature(
1524             new SequenceFeature("type2", "desc2", 3, 4, 2f, null));
1525     // exact overlap - to [1, 3]
1526     dna.addSequenceFeature(
1527             new SequenceFeature("type3", "desc3", 4, 6, 3f, null));
1528     // spanning overlap - to [2, 5]
1529     dna.addSequenceFeature(
1530             new SequenceFeature("type4", "desc4", 5, 11, 4f, null));
1531     // exactly overlaps whole mapped range [1, 6]
1532     dna.addSequenceFeature(
1533             new SequenceFeature("type5", "desc5", 4, 12, 5f, null));
1534     // no overlap (internal)
1535     dna.addSequenceFeature(
1536             new SequenceFeature("type6", "desc6", 7, 9, 6f, null));
1537     // no overlap (3' end)
1538     dna.addSequenceFeature(
1539             new SequenceFeature("type7", "desc7", 13, 15, 7f, null));
1540     // overlap (3' end) - to [6, 6]
1541     dna.addSequenceFeature(
1542             new SequenceFeature("type8", "desc8", 12, 12, 8f, null));
1543     // extended overlap - to [6, +]
1544     dna.addSequenceFeature(
1545             new SequenceFeature("type9", "desc9", 12, 13, 9f, null));
1546
1547     MapList map = new MapList(new int[] { 4, 6, 10, 12 },
1548             new int[]
1549             { 1, 6 }, 1, 1);
1550
1551     /*
1552      * transferFeatures() will build 'partial overlap' for regions
1553      * that partially overlap 5' or 3' (start or end) of target sequence
1554      */
1555     AlignmentUtils.transferFeatures(dna, cds, map, null);
1556     List<SequenceFeature> sfs = cds.getSequenceFeatures();
1557     assertEquals(6, sfs.size());
1558
1559     SequenceFeature sf = sfs.get(0);
1560     assertEquals("type2", sf.getType());
1561     assertEquals("desc2", sf.getDescription());
1562     assertEquals(2f, sf.getScore());
1563     assertEquals(1, sf.getBegin());
1564     assertEquals(1, sf.getEnd());
1565
1566     sf = sfs.get(1);
1567     assertEquals("type3", sf.getType());
1568     assertEquals("desc3", sf.getDescription());
1569     assertEquals(3f, sf.getScore());
1570     assertEquals(1, sf.getBegin());
1571     assertEquals(3, sf.getEnd());
1572
1573     sf = sfs.get(2);
1574     assertEquals("type4", sf.getType());
1575     assertEquals(2, sf.getBegin());
1576     assertEquals(5, sf.getEnd());
1577
1578     sf = sfs.get(3);
1579     assertEquals("type5", sf.getType());
1580     assertEquals(1, sf.getBegin());
1581     assertEquals(6, sf.getEnd());
1582
1583     sf = sfs.get(4);
1584     assertEquals("type8", sf.getType());
1585     assertEquals(6, sf.getBegin());
1586     assertEquals(6, sf.getEnd());
1587
1588     sf = sfs.get(5);
1589     assertEquals("type9", sf.getType());
1590     assertEquals(6, sf.getBegin());
1591     assertEquals(6, sf.getEnd());
1592   }
1593
1594   /**
1595    * Tests for transferring features between mapped sequences
1596    */
1597   @Test(groups = { "Functional" })
1598   public void testTransferFeatures_withOmit()
1599   {
1600     SequenceI dna = new Sequence("dna/20-34", "acgTAGcaaGCCcgt");
1601     SequenceI cds = new Sequence("cds/10-15", "TAGGCC");
1602
1603     MapList map = new MapList(new int[] { 4, 6, 10, 12 },
1604             new int[]
1605             { 1, 6 }, 1, 1);
1606
1607     // [5, 11] maps to [2, 5]
1608     dna.addSequenceFeature(
1609             new SequenceFeature("type4", "desc4", 5, 11, 4f, null));
1610     // [4, 12] maps to [1, 6]
1611     dna.addSequenceFeature(
1612             new SequenceFeature("type5", "desc5", 4, 12, 5f, null));
1613     // [12, 12] maps to [6, 6]
1614     dna.addSequenceFeature(
1615             new SequenceFeature("type8", "desc8", 12, 12, 8f, null));
1616
1617     // desc4 and desc8 are the 'omit these' varargs
1618     AlignmentUtils.transferFeatures(dna, cds, map, null, "type4", "type8");
1619     List<SequenceFeature> sfs = cds.getSequenceFeatures();
1620     assertEquals(1, sfs.size());
1621
1622     SequenceFeature sf = sfs.get(0);
1623     assertEquals("type5", sf.getType());
1624     assertEquals(1, sf.getBegin());
1625     assertEquals(6, sf.getEnd());
1626   }
1627
1628   /**
1629    * Tests for transferring features between mapped sequences
1630    */
1631   @Test(groups = { "Functional" })
1632   public void testTransferFeatures_withSelect()
1633   {
1634     SequenceI dna = new Sequence("dna/20-34", "acgTAGcaaGCCcgt");
1635     SequenceI cds = new Sequence("cds/10-15", "TAGGCC");
1636
1637     MapList map = new MapList(new int[] { 4, 6, 10, 12 },
1638             new int[]
1639             { 1, 6 }, 1, 1);
1640
1641     // [5, 11] maps to [2, 5]
1642     dna.addSequenceFeature(
1643             new SequenceFeature("type4", "desc4", 5, 11, 4f, null));
1644     // [4, 12] maps to [1, 6]
1645     dna.addSequenceFeature(
1646             new SequenceFeature("type5", "desc5", 4, 12, 5f, null));
1647     // [12, 12] maps to [6, 6]
1648     dna.addSequenceFeature(
1649             new SequenceFeature("type8", "desc8", 12, 12, 8f, null));
1650
1651     // "type5" is the 'select this type' argument
1652     AlignmentUtils.transferFeatures(dna, cds, map, "type5");
1653     List<SequenceFeature> sfs = cds.getSequenceFeatures();
1654     assertEquals(1, sfs.size());
1655
1656     SequenceFeature sf = sfs.get(0);
1657     assertEquals("type5", sf.getType());
1658     assertEquals(1, sf.getBegin());
1659     assertEquals(6, sf.getEnd());
1660   }
1661
1662   /**
1663    * Test the method that extracts the cds-only part of a dna alignment, for the
1664    * case where the cds should be aligned to match its nucleotide sequence.
1665    */
1666   @Test(groups = { "Functional" })
1667   public void testMakeCdsAlignment_alternativeTranscripts()
1668   {
1669     SequenceI dna1 = new Sequence("dna1", "aaaGGGCC-----CTTTaaaGGG");
1670     // alternative transcript of same dna skips CCC codon
1671     SequenceI dna2 = new Sequence("dna2", "aaaGGGCC-----cttTaaaGGG");
1672     // dna3 has no mapping (protein product) so should be ignored here
1673     SequenceI dna3 = new Sequence("dna3", "aaaGGGCCCCCGGGcttTaaaGGG");
1674     SequenceI pep1 = new Sequence("pep1", "GPFG");
1675     SequenceI pep2 = new Sequence("pep2", "GPG");
1676     dna1.createDatasetSequence();
1677     dna2.createDatasetSequence();
1678     dna3.createDatasetSequence();
1679     pep1.createDatasetSequence();
1680     pep2.createDatasetSequence();
1681
1682     AlignmentI dna = new Alignment(new SequenceI[] { dna1, dna2, dna3 });
1683     dna.setDataset(null);
1684
1685     MapList map = new MapList(new int[] { 4, 12, 16, 18 },
1686             new int[]
1687             { 1, 4 }, 3, 1);
1688     AlignedCodonFrame acf = new AlignedCodonFrame();
1689     acf.addMap(dna1.getDatasetSequence(), pep1.getDatasetSequence(), map);
1690     dna.addCodonFrame(acf);
1691     map = new MapList(new int[] { 4, 8, 12, 12, 16, 18 },
1692             new int[]
1693             { 1, 3 }, 3, 1);
1694     acf = new AlignedCodonFrame();
1695     acf.addMap(dna2.getDatasetSequence(), pep2.getDatasetSequence(), map);
1696     dna.addCodonFrame(acf);
1697
1698     AlignmentI cds = AlignmentUtils
1699             .makeCdsAlignment(new SequenceI[]
1700             { dna1, dna2, dna3 }, dna.getDataset(), null);
1701     List<SequenceI> cdsSeqs = cds.getSequences();
1702     assertEquals(2, cdsSeqs.size());
1703     assertEquals("GGGCCCTTTGGG", cdsSeqs.get(0).getSequenceAsString());
1704     assertEquals("GGGCCTGGG", cdsSeqs.get(1).getSequenceAsString());
1705
1706     /*
1707      * verify shared, extended alignment dataset
1708      */
1709     assertSame(dna.getDataset(), cds.getDataset());
1710     assertTrue(dna.getDataset().getSequences()
1711             .contains(cdsSeqs.get(0).getDatasetSequence()));
1712     assertTrue(dna.getDataset().getSequences()
1713             .contains(cdsSeqs.get(1).getDatasetSequence()));
1714
1715     /*
1716      * Verify 6 mappings: dna1 to cds1, cds1 to pep1, dna1 to pep1
1717      * and the same for dna2/cds2/pep2
1718      */
1719     List<AlignedCodonFrame> mappings = cds.getCodonFrames();
1720     assertEquals(6, mappings.size());
1721
1722     /*
1723      * 2 mappings involve pep1
1724      */
1725     List<AlignedCodonFrame> pep1Mappings = MappingUtils
1726             .findMappingsForSequence(pep1, mappings);
1727     assertEquals(2, pep1Mappings.size());
1728
1729     /*
1730      * Get mapping of pep1 to cds1 and verify it
1731      * maps GPFG to 1-3,4-6,7-9,10-12
1732      */
1733     List<AlignedCodonFrame> pep1CdsMappings = MappingUtils
1734             .findMappingsForSequence(cds.getSequenceAt(0), pep1Mappings);
1735     assertEquals(1, pep1CdsMappings.size());
1736     SearchResultsI sr = MappingUtils.buildSearchResults(pep1, 1,
1737             pep1CdsMappings);
1738     assertEquals(1, sr.getResults().size());
1739     SearchResultMatchI m = sr.getResults().get(0);
1740     assertEquals(cds.getSequenceAt(0).getDatasetSequence(),
1741             m.getSequence());
1742     assertEquals(1, m.getStart());
1743     assertEquals(3, m.getEnd());
1744     sr = MappingUtils.buildSearchResults(pep1, 2, pep1CdsMappings);
1745     m = sr.getResults().get(0);
1746     assertEquals(4, m.getStart());
1747     assertEquals(6, m.getEnd());
1748     sr = MappingUtils.buildSearchResults(pep1, 3, pep1CdsMappings);
1749     m = sr.getResults().get(0);
1750     assertEquals(7, m.getStart());
1751     assertEquals(9, m.getEnd());
1752     sr = MappingUtils.buildSearchResults(pep1, 4, pep1CdsMappings);
1753     m = sr.getResults().get(0);
1754     assertEquals(10, m.getStart());
1755     assertEquals(12, m.getEnd());
1756
1757     /*
1758      * Get mapping of pep2 to cds2 and verify it
1759      * maps GPG in pep2 to 1-3,4-6,7-9 in second CDS sequence
1760      */
1761     List<AlignedCodonFrame> pep2Mappings = MappingUtils
1762             .findMappingsForSequence(pep2, mappings);
1763     assertEquals(2, pep2Mappings.size());
1764     List<AlignedCodonFrame> pep2CdsMappings = MappingUtils
1765             .findMappingsForSequence(cds.getSequenceAt(1), pep2Mappings);
1766     assertEquals(1, pep2CdsMappings.size());
1767     sr = MappingUtils.buildSearchResults(pep2, 1, pep2CdsMappings);
1768     assertEquals(1, sr.getResults().size());
1769     m = sr.getResults().get(0);
1770     assertEquals(cds.getSequenceAt(1).getDatasetSequence(),
1771             m.getSequence());
1772     assertEquals(1, m.getStart());
1773     assertEquals(3, m.getEnd());
1774     sr = MappingUtils.buildSearchResults(pep2, 2, pep2CdsMappings);
1775     m = sr.getResults().get(0);
1776     assertEquals(4, m.getStart());
1777     assertEquals(6, m.getEnd());
1778     sr = MappingUtils.buildSearchResults(pep2, 3, pep2CdsMappings);
1779     m = sr.getResults().get(0);
1780     assertEquals(7, m.getStart());
1781     assertEquals(9, m.getEnd());
1782   }
1783
1784   /**
1785    * Test the method that realigns protein to match mapped codon alignment.
1786    */
1787   @Test(groups = { "Functional" })
1788   public void testAlignProteinAsDna_incompleteStartCodon()
1789   {
1790     // seq1: incomplete start codon (not mapped), then [3, 11]
1791     SequenceI dna1 = new Sequence("Seq1", "ccAAA-TTT-GGG-");
1792     // seq2 codons are [4, 5], [8, 11]
1793     SequenceI dna2 = new Sequence("Seq2", "ccaAA-ttT-GGG-");
1794     // seq3 incomplete start codon at 'tt'
1795     SequenceI dna3 = new Sequence("Seq3", "ccaaa-ttt-GGG-");
1796     AlignmentI dna = new Alignment(new SequenceI[] { dna1, dna2, dna3 });
1797     dna.setDataset(null);
1798
1799     // prot1 has 'X' for incomplete start codon (not mapped)
1800     SequenceI prot1 = new Sequence("Seq1", "XKFG"); // X for incomplete start
1801     SequenceI prot2 = new Sequence("Seq2", "NG");
1802     SequenceI prot3 = new Sequence("Seq3", "XG"); // X for incomplete start
1803     AlignmentI protein = new Alignment(
1804             new SequenceI[]
1805             { prot1, prot2, prot3 });
1806     protein.setDataset(null);
1807
1808     // map dna1 [3, 11] to prot1 [2, 4] KFG
1809     MapList map = new MapList(new int[] { 3, 11 }, new int[] { 2, 4 }, 3,
1810             1);
1811     AlignedCodonFrame acf = new AlignedCodonFrame();
1812     acf.addMap(dna1.getDatasetSequence(), prot1.getDatasetSequence(), map);
1813
1814     // map dna2 [4, 5] [8, 11] to prot2 [1, 2] NG
1815     map = new MapList(new int[] { 4, 5, 8, 11 }, new int[] { 1, 2 }, 3, 1);
1816     acf.addMap(dna2.getDatasetSequence(), prot2.getDatasetSequence(), map);
1817
1818     // map dna3 [9, 11] to prot3 [2, 2] G
1819     map = new MapList(new int[] { 9, 11 }, new int[] { 2, 2 }, 3, 1);
1820     acf.addMap(dna3.getDatasetSequence(), prot3.getDatasetSequence(), map);
1821
1822     ArrayList<AlignedCodonFrame> acfs = new ArrayList<>();
1823     acfs.add(acf);
1824     protein.setCodonFrames(acfs);
1825
1826     /*
1827      * verify X is included in the aligned proteins, and placed just
1828      * before the first mapped residue
1829      * CCT is between CCC and TTT
1830      */
1831     AlignmentUtils.alignProteinAsDna(protein, dna);
1832     assertEquals("XK-FG", prot1.getSequenceAsString());
1833     assertEquals("--N-G", prot2.getSequenceAsString());
1834     assertEquals("---XG", prot3.getSequenceAsString());
1835   }
1836
1837   /**
1838    * Tests for the method that maps the subset of a dna sequence that has CDS
1839    * (or subtype) feature - case where the start codon is incomplete.
1840    */
1841   @Test(groups = "Functional")
1842   public void testFindCdsPositions_fivePrimeIncomplete()
1843   {
1844     SequenceI dnaSeq = new Sequence("dna", "aaagGGCCCaaaTTTttt");
1845     dnaSeq.createDatasetSequence();
1846     SequenceI ds = dnaSeq.getDatasetSequence();
1847
1848     // CDS for dna 5-6 (incomplete codon), 7-9
1849     SequenceFeature sf = new SequenceFeature("CDS", "", 5, 9, 0f, null);
1850     sf.setPhase("2"); // skip 2 bases to start of next codon
1851     ds.addSequenceFeature(sf);
1852     // CDS for dna 13-15
1853     sf = new SequenceFeature("CDS_predicted", "", 13, 15, 0f, null);
1854     ds.addSequenceFeature(sf);
1855
1856     List<int[]> ranges = AlignmentUtils.findCdsPositions(dnaSeq);
1857
1858     /*
1859      * check the mapping starts with the first complete codon
1860      */
1861     assertEquals(6, MappingUtils.getLength(ranges));
1862     assertEquals(2, ranges.size());
1863     assertEquals(7, ranges.get(0)[0]);
1864     assertEquals(9, ranges.get(0)[1]);
1865     assertEquals(13, ranges.get(1)[0]);
1866     assertEquals(15, ranges.get(1)[1]);
1867   }
1868
1869   /**
1870    * Tests for the method that maps the subset of a dna sequence that has CDS
1871    * (or subtype) feature.
1872    */
1873   @Test(groups = "Functional")
1874   public void testFindCdsPositions()
1875   {
1876     SequenceI dnaSeq = new Sequence("dna", "aaaGGGcccAAATTTttt");
1877     dnaSeq.createDatasetSequence();
1878     SequenceI ds = dnaSeq.getDatasetSequence();
1879
1880     // CDS for dna 10-12
1881     SequenceFeature sf = new SequenceFeature("CDS_predicted", "", 10, 12,
1882             0f, null);
1883     sf.setStrand("+");
1884     ds.addSequenceFeature(sf);
1885     // CDS for dna 4-6
1886     sf = new SequenceFeature("CDS", "", 4, 6, 0f, null);
1887     sf.setStrand("+");
1888     ds.addSequenceFeature(sf);
1889     // exon feature should be ignored here
1890     sf = new SequenceFeature("exon", "", 7, 9, 0f, null);
1891     ds.addSequenceFeature(sf);
1892
1893     List<int[]> ranges = AlignmentUtils.findCdsPositions(dnaSeq);
1894     /*
1895      * verify ranges { [4-6], [12-10] }
1896      * note CDS ranges are ordered ascending even if the CDS
1897      * features are not
1898      */
1899     assertEquals(6, MappingUtils.getLength(ranges));
1900     assertEquals(2, ranges.size());
1901     assertEquals(4, ranges.get(0)[0]);
1902     assertEquals(6, ranges.get(0)[1]);
1903     assertEquals(10, ranges.get(1)[0]);
1904     assertEquals(12, ranges.get(1)[1]);
1905   }
1906
1907   /**
1908    * Tests for the method that maps the subset of a dna sequence that has CDS
1909    * (or subtype) feature, with CDS strand = '-' (reverse)
1910    */
1911   // test turned off as currently findCdsPositions is not strand-dependent
1912   // left in case it comes around again...
1913   @Test(groups = "Functional", enabled = false)
1914   public void testFindCdsPositions_reverseStrand()
1915   {
1916     SequenceI dnaSeq = new Sequence("dna", "aaaGGGcccAAATTTttt");
1917     dnaSeq.createDatasetSequence();
1918     SequenceI ds = dnaSeq.getDatasetSequence();
1919
1920     // CDS for dna 4-6
1921     SequenceFeature sf = new SequenceFeature("CDS", "", 4, 6, 0f, null);
1922     sf.setStrand("-");
1923     ds.addSequenceFeature(sf);
1924     // exon feature should be ignored here
1925     sf = new SequenceFeature("exon", "", 7, 9, 0f, null);
1926     ds.addSequenceFeature(sf);
1927     // CDS for dna 10-12
1928     sf = new SequenceFeature("CDS_predicted", "", 10, 12, 0f, null);
1929     sf.setStrand("-");
1930     ds.addSequenceFeature(sf);
1931
1932     List<int[]> ranges = AlignmentUtils.findCdsPositions(dnaSeq);
1933     /*
1934      * verify ranges { [12-10], [6-4] }
1935      */
1936     assertEquals(6, MappingUtils.getLength(ranges));
1937     assertEquals(2, ranges.size());
1938     assertEquals(12, ranges.get(0)[0]);
1939     assertEquals(10, ranges.get(0)[1]);
1940     assertEquals(6, ranges.get(1)[0]);
1941     assertEquals(4, ranges.get(1)[1]);
1942   }
1943
1944   /**
1945    * Tests for the method that maps the subset of a dna sequence that has CDS
1946    * (or subtype) feature - reverse strand case where the start codon is
1947    * incomplete.
1948    */
1949   @Test(groups = "Functional", enabled = false)
1950   // test turned off as currently findCdsPositions is not strand-dependent
1951   // left in case it comes around again...
1952   public void testFindCdsPositions_reverseStrandThreePrimeIncomplete()
1953   {
1954     SequenceI dnaSeq = new Sequence("dna", "aaagGGCCCaaaTTTttt");
1955     dnaSeq.createDatasetSequence();
1956     SequenceI ds = dnaSeq.getDatasetSequence();
1957
1958     // CDS for dna 5-9
1959     SequenceFeature sf = new SequenceFeature("CDS", "", 5, 9, 0f, null);
1960     sf.setStrand("-");
1961     ds.addSequenceFeature(sf);
1962     // CDS for dna 13-15
1963     sf = new SequenceFeature("CDS_predicted", "", 13, 15, 0f, null);
1964     sf.setStrand("-");
1965     sf.setPhase("2"); // skip 2 bases to start of next codon
1966     ds.addSequenceFeature(sf);
1967
1968     List<int[]> ranges = AlignmentUtils.findCdsPositions(dnaSeq);
1969
1970     /*
1971      * check the mapping starts with the first complete codon
1972      * expect ranges [13, 13], [9, 5]
1973      */
1974     assertEquals(6, MappingUtils.getLength(ranges));
1975     assertEquals(2, ranges.size());
1976     assertEquals(13, ranges.get(0)[0]);
1977     assertEquals(13, ranges.get(0)[1]);
1978     assertEquals(9, ranges.get(1)[0]);
1979     assertEquals(5, ranges.get(1)[1]);
1980   }
1981
1982   @Test(groups = "Functional")
1983   public void testAlignAs_alternateTranscriptsUngapped()
1984   {
1985     SequenceI dna1 = new Sequence("dna1", "cccGGGTTTaaa");
1986     SequenceI dna2 = new Sequence("dna2", "CCCgggtttAAA");
1987     AlignmentI dna = new Alignment(new SequenceI[] { dna1, dna2 });
1988     ((Alignment) dna).createDatasetAlignment();
1989     SequenceI cds1 = new Sequence("cds1", "GGGTTT");
1990     SequenceI cds2 = new Sequence("cds2", "CCCAAA");
1991     AlignmentI cds = new Alignment(new SequenceI[] { cds1, cds2 });
1992     ((Alignment) cds).createDatasetAlignment();
1993
1994     AlignedCodonFrame acf = new AlignedCodonFrame();
1995     MapList map = new MapList(new int[] { 4, 9 }, new int[] { 1, 6 }, 1, 1);
1996     acf.addMap(dna1.getDatasetSequence(), cds1.getDatasetSequence(), map);
1997     map = new MapList(new int[] { 1, 3, 10, 12 }, new int[] { 1, 6 }, 1, 1);
1998     acf.addMap(dna2.getDatasetSequence(), cds2.getDatasetSequence(), map);
1999
2000     /*
2001      * verify CDS alignment is as:
2002      *   cccGGGTTTaaa (cdna)
2003      *   CCCgggtttAAA (cdna)
2004      *
2005      *   ---GGGTTT--- (cds)
2006      *   CCC------AAA (cds)
2007      */
2008     dna.addCodonFrame(acf);
2009     AlignmentUtils.alignAs(cds, dna);
2010     assertEquals("---GGGTTT", cds.getSequenceAt(0).getSequenceAsString());
2011     assertEquals("CCC------AAA",
2012             cds.getSequenceAt(1).getSequenceAsString());
2013   }
2014
2015   @Test(groups = { "Functional" })
2016   public void testAddMappedPositions()
2017   {
2018     SequenceI from = new Sequence("dna", "ggAA-ATcc-TT-g");
2019     SequenceI seq1 = new Sequence("cds", "AAATTT");
2020     from.createDatasetSequence();
2021     seq1.createDatasetSequence();
2022     Mapping mapping = new Mapping(seq1,
2023             new MapList(new int[]
2024             { 3, 6, 9, 10 }, new int[] { 1, 6 }, 1, 1));
2025     Map<Integer, Map<SequenceI, Character>> map = new TreeMap<>();
2026     AlignmentUtils.addMappedPositions(seq1, from, mapping, map);
2027
2028     /*
2029      * verify map has seq1 residues in columns 3,4,6,7,11,12
2030      */
2031     assertEquals(6, map.size());
2032     assertEquals('A', map.get(3).get(seq1).charValue());
2033     assertEquals('A', map.get(4).get(seq1).charValue());
2034     assertEquals('A', map.get(6).get(seq1).charValue());
2035     assertEquals('T', map.get(7).get(seq1).charValue());
2036     assertEquals('T', map.get(11).get(seq1).charValue());
2037     assertEquals('T', map.get(12).get(seq1).charValue());
2038
2039     /*
2040      *
2041      */
2042   }
2043
2044   /**
2045    * Test case where the mapping 'from' range includes a stop codon which is
2046    * absent in the 'to' range
2047    */
2048   @Test(groups = { "Functional" })
2049   public void testAddMappedPositions_withStopCodon()
2050   {
2051     SequenceI from = new Sequence("dna", "ggAA-ATcc-TT-g");
2052     SequenceI seq1 = new Sequence("cds", "AAATTT");
2053     from.createDatasetSequence();
2054     seq1.createDatasetSequence();
2055     Mapping mapping = new Mapping(seq1,
2056             new MapList(new int[]
2057             { 3, 6, 9, 10 }, new int[] { 1, 6 }, 1, 1));
2058     Map<Integer, Map<SequenceI, Character>> map = new TreeMap<>();
2059     AlignmentUtils.addMappedPositions(seq1, from, mapping, map);
2060
2061     /*
2062      * verify map has seq1 residues in columns 3,4,6,7,11,12
2063      */
2064     assertEquals(6, map.size());
2065     assertEquals('A', map.get(3).get(seq1).charValue());
2066     assertEquals('A', map.get(4).get(seq1).charValue());
2067     assertEquals('A', map.get(6).get(seq1).charValue());
2068     assertEquals('T', map.get(7).get(seq1).charValue());
2069     assertEquals('T', map.get(11).get(seq1).charValue());
2070     assertEquals('T', map.get(12).get(seq1).charValue());
2071   }
2072
2073   /**
2074    * Test for the case where the products for which we want CDS are specified.
2075    * This is to represent the case where EMBL has CDS mappings to both Uniprot
2076    * and EMBLCDSPROTEIN. makeCdsAlignment() should only return the mappings for
2077    * the protein sequences specified.
2078    */
2079   @Test(groups = { "Functional" })
2080   public void testMakeCdsAlignment_filterProducts()
2081   {
2082     SequenceI dna1 = new Sequence("dna1", "aaaGGGcccTTTaaa");
2083     SequenceI dna2 = new Sequence("dna2", "GGGcccTTTaaaCCC");
2084     SequenceI pep1 = new Sequence("Uniprot|pep1", "GF");
2085     SequenceI pep2 = new Sequence("Uniprot|pep2", "GFP");
2086     SequenceI pep3 = new Sequence("EMBL|pep3", "GF");
2087     SequenceI pep4 = new Sequence("EMBL|pep4", "GFP");
2088     dna1.createDatasetSequence();
2089     dna2.createDatasetSequence();
2090     pep1.createDatasetSequence();
2091     pep2.createDatasetSequence();
2092     pep3.createDatasetSequence();
2093     pep4.createDatasetSequence();
2094     AlignmentI dna = new Alignment(new SequenceI[] { dna1, dna2 });
2095     dna.setDataset(null);
2096     AlignmentI emblPeptides = new Alignment(new SequenceI[] { pep3, pep4 });
2097     emblPeptides.setDataset(null);
2098
2099     AlignedCodonFrame acf = new AlignedCodonFrame();
2100     MapList map = new MapList(new int[] { 4, 6, 10, 12 },
2101             new int[]
2102             { 1, 2 }, 3, 1);
2103     acf.addMap(dna1.getDatasetSequence(), pep1.getDatasetSequence(), map);
2104     acf.addMap(dna1.getDatasetSequence(), pep3.getDatasetSequence(), map);
2105     dna.addCodonFrame(acf);
2106
2107     acf = new AlignedCodonFrame();
2108     map = new MapList(new int[] { 1, 3, 7, 9, 13, 15 }, new int[] { 1, 3 },
2109             3, 1);
2110     acf.addMap(dna2.getDatasetSequence(), pep2.getDatasetSequence(), map);
2111     acf.addMap(dna2.getDatasetSequence(), pep4.getDatasetSequence(), map);
2112     dna.addCodonFrame(acf);
2113
2114     /*
2115      * execute method under test to find CDS for EMBL peptides only
2116      */
2117     AlignmentI cds = AlignmentUtils
2118             .makeCdsAlignment(new SequenceI[]
2119             { dna1, dna2 }, dna.getDataset(),
2120                     emblPeptides.getSequencesArray());
2121
2122     assertEquals(2, cds.getSequences().size());
2123     assertEquals("GGGTTT", cds.getSequenceAt(0).getSequenceAsString());
2124     assertEquals("GGGTTTCCC", cds.getSequenceAt(1).getSequenceAsString());
2125
2126     /*
2127      * verify shared, extended alignment dataset
2128      */
2129     assertSame(dna.getDataset(), cds.getDataset());
2130     assertTrue(dna.getDataset().getSequences()
2131             .contains(cds.getSequenceAt(0).getDatasetSequence()));
2132     assertTrue(dna.getDataset().getSequences()
2133             .contains(cds.getSequenceAt(1).getDatasetSequence()));
2134
2135     /*
2136      * Verify mappings from CDS to peptide, cDNA to CDS, and cDNA to peptide
2137      * the mappings are on the shared alignment dataset
2138      */
2139     List<AlignedCodonFrame> cdsMappings = cds.getDataset().getCodonFrames();
2140     /*
2141      * 6 mappings, 2*(DNA->CDS), 2*(DNA->Pep), 2*(CDS->Pep)
2142      */
2143     assertEquals(6, cdsMappings.size());
2144
2145     /*
2146      * verify that mapping sets for dna and cds alignments are different
2147      * [not current behaviour - all mappings are on the alignment dataset]
2148      */
2149     // select -> subselect type to test.
2150     // Assert.assertNotSame(dna.getCodonFrames(), cds.getCodonFrames());
2151     // assertEquals(4, dna.getCodonFrames().size());
2152     // assertEquals(4, cds.getCodonFrames().size());
2153
2154     /*
2155      * Two mappings involve pep3 (dna to pep3, cds to pep3)
2156      * Mapping from pep3 to GGGTTT in first new exon sequence
2157      */
2158     List<AlignedCodonFrame> pep3Mappings = MappingUtils
2159             .findMappingsForSequence(pep3, cdsMappings);
2160     assertEquals(2, pep3Mappings.size());
2161     List<AlignedCodonFrame> mappings = MappingUtils
2162             .findMappingsForSequence(cds.getSequenceAt(0), pep3Mappings);
2163     assertEquals(1, mappings.size());
2164
2165     // map G to GGG
2166     SearchResultsI sr = MappingUtils.buildSearchResults(pep3, 1, mappings);
2167     assertEquals(1, sr.getResults().size());
2168     SearchResultMatchI m = sr.getResults().get(0);
2169     assertSame(cds.getSequenceAt(0).getDatasetSequence(), m.getSequence());
2170     assertEquals(1, m.getStart());
2171     assertEquals(3, m.getEnd());
2172     // map F to TTT
2173     sr = MappingUtils.buildSearchResults(pep3, 2, mappings);
2174     m = sr.getResults().get(0);
2175     assertSame(cds.getSequenceAt(0).getDatasetSequence(), m.getSequence());
2176     assertEquals(4, m.getStart());
2177     assertEquals(6, m.getEnd());
2178
2179     /*
2180      * Two mappings involve pep4 (dna to pep4, cds to pep4)
2181      * Verify mapping from pep4 to GGGTTTCCC in second new exon sequence
2182      */
2183     List<AlignedCodonFrame> pep4Mappings = MappingUtils
2184             .findMappingsForSequence(pep4, cdsMappings);
2185     assertEquals(2, pep4Mappings.size());
2186     mappings = MappingUtils.findMappingsForSequence(cds.getSequenceAt(1),
2187             pep4Mappings);
2188     assertEquals(1, mappings.size());
2189     // map G to GGG
2190     sr = MappingUtils.buildSearchResults(pep4, 1, mappings);
2191     assertEquals(1, sr.getResults().size());
2192     m = sr.getResults().get(0);
2193     assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
2194     assertEquals(1, m.getStart());
2195     assertEquals(3, m.getEnd());
2196     // map F to TTT
2197     sr = MappingUtils.buildSearchResults(pep4, 2, mappings);
2198     m = sr.getResults().get(0);
2199     assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
2200     assertEquals(4, m.getStart());
2201     assertEquals(6, m.getEnd());
2202     // map P to CCC
2203     sr = MappingUtils.buildSearchResults(pep4, 3, mappings);
2204     m = sr.getResults().get(0);
2205     assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
2206     assertEquals(7, m.getStart());
2207     assertEquals(9, m.getEnd());
2208   }
2209
2210   /**
2211    * Test the method that just copies aligned sequences, provided all sequences
2212    * to be aligned share the aligned sequence's dataset
2213    */
2214   @Test(groups = "Functional")
2215   public void testAlignAsSameSequences()
2216   {
2217     SequenceI dna1 = new Sequence("dna1", "cccGGGTTTaaa");
2218     SequenceI dna2 = new Sequence("dna2", "CCCgggtttAAA");
2219     AlignmentI al1 = new Alignment(new SequenceI[] { dna1, dna2 });
2220     ((Alignment) al1).createDatasetAlignment();
2221
2222     SequenceI dna3 = new Sequence(dna1);
2223     SequenceI dna4 = new Sequence(dna2);
2224     assertSame(dna3.getDatasetSequence(), dna1.getDatasetSequence());
2225     assertSame(dna4.getDatasetSequence(), dna2.getDatasetSequence());
2226     String seq1 = "-cc-GG-GT-TT--aaa";
2227     dna3.setSequence(seq1);
2228     String seq2 = "C--C-Cgg--gtt-tAA-A-";
2229     dna4.setSequence(seq2);
2230     AlignmentI al2 = new Alignment(new SequenceI[] { dna3, dna4 });
2231     ((Alignment) al2).createDatasetAlignment();
2232
2233     /*
2234      * alignment removes gapped columns (two internal, two trailing)
2235      */
2236     assertTrue(AlignmentUtils.alignAsSameSequences(al1, al2));
2237     String aligned1 = "-cc-GG-GTTT-aaa";
2238     assertEquals(aligned1, al1.getSequenceAt(0).getSequenceAsString());
2239     String aligned2 = "C--C-Cgg-gtttAAA";
2240     assertEquals(aligned2, al1.getSequenceAt(1).getSequenceAsString());
2241
2242     /*
2243      * add another sequence to 'aligned' - should still succeed, since
2244      * unaligned sequences still share a dataset with aligned sequences
2245      */
2246     SequenceI dna5 = new Sequence("dna5", "CCCgggtttAAA");
2247     dna5.createDatasetSequence();
2248     al2.addSequence(dna5);
2249     assertTrue(AlignmentUtils.alignAsSameSequences(al1, al2));
2250     assertEquals(aligned1, al1.getSequenceAt(0).getSequenceAsString());
2251     assertEquals(aligned2, al1.getSequenceAt(1).getSequenceAsString());
2252
2253     /*
2254      * add another sequence to 'unaligned' - should fail, since now not
2255      * all unaligned sequences share a dataset with aligned sequences
2256      */
2257     SequenceI dna6 = new Sequence("dna6", "CCCgggtttAAA");
2258     dna6.createDatasetSequence();
2259     al1.addSequence(dna6);
2260     // JAL-2110 JBP Comment: what's the use case for this behaviour ?
2261     assertFalse(AlignmentUtils.alignAsSameSequences(al1, al2));
2262   }
2263
2264   @Test(groups = "Functional")
2265   public void testAlignAsSameSequencesMultipleSubSeq()
2266   {
2267     SequenceI dna1 = new Sequence("dna1", "cccGGGTTTaaa");
2268     SequenceI dna2 = new Sequence("dna2", "CCCgggtttAAA");
2269     SequenceI as1 = dna1.deriveSequence(); // cccGGGTTTaaa/1-12
2270     SequenceI as2 = dna1.deriveSequence().getSubSequence(3, 7); // GGGT/4-7
2271     SequenceI as3 = dna2.deriveSequence(); // CCCgggtttAAA/1-12
2272     as1.insertCharAt(6, 5, '-');
2273     assertEquals("cccGGG-----TTTaaa", as1.getSequenceAsString());
2274     as2.insertCharAt(6, 5, '-');
2275     assertEquals("GGGT-----", as2.getSequenceAsString());
2276     as3.insertCharAt(3, 5, '-');
2277     assertEquals("CCC-----gggtttAAA", as3.getSequenceAsString());
2278     AlignmentI aligned = new Alignment(new SequenceI[] { as1, as2, as3 });
2279
2280     // why do we need to cast this still ?
2281     ((Alignment) aligned).createDatasetAlignment();
2282     SequenceI uas1 = dna1.deriveSequence();
2283     SequenceI uas2 = dna1.deriveSequence().getSubSequence(3, 7);
2284     SequenceI uas3 = dna2.deriveSequence();
2285     AlignmentI tobealigned = new Alignment(
2286             new SequenceI[]
2287             { uas1, uas2, uas3 });
2288     ((Alignment) tobealigned).createDatasetAlignment();
2289
2290     /*
2291      * alignAs lines up dataset sequences and removes empty columns (two)
2292      */
2293     assertTrue(AlignmentUtils.alignAsSameSequences(tobealigned, aligned));
2294     assertEquals("cccGGG---TTTaaa", uas1.getSequenceAsString());
2295     assertEquals("GGGT", uas2.getSequenceAsString());
2296     assertEquals("CCC---gggtttAAA", uas3.getSequenceAsString());
2297   }
2298
2299   @Test(groups = { "Functional" })
2300   public void testTransferGeneLoci()
2301   {
2302     SequenceI from = new Sequence("transcript",
2303             "aaacccgggTTTAAACCCGGGtttaaacccgggttt");
2304     SequenceI to = new Sequence("CDS", "TTTAAACCCGGG");
2305     MapList map = new MapList(new int[] { 1, 12 }, new int[] { 10, 21 }, 1,
2306             1);
2307
2308     /*
2309      * first with nothing to transfer
2310      */
2311     AlignmentUtils.transferGeneLoci(from, map, to);
2312     assertNull(to.getGeneLoci());
2313
2314     /*
2315      * next with gene loci set on 'from' sequence
2316      */
2317     int[] exons = new int[] { 100, 105, 155, 164, 210, 229 };
2318     MapList geneMap = new MapList(new int[] { 1, 36 }, exons, 1, 1);
2319     from.setGeneLoci("human", "GRCh38", "7", geneMap);
2320     AlignmentUtils.transferGeneLoci(from, map, to);
2321
2322     GeneLociI toLoci = to.getGeneLoci();
2323     assertNotNull(toLoci);
2324     // DBRefEntry constructor upper-cases 'source'
2325     assertEquals("HUMAN", toLoci.getSpeciesId());
2326     assertEquals("GRCh38", toLoci.getAssemblyId());
2327     assertEquals("7", toLoci.getChromosomeId());
2328
2329     /*
2330      * transcript 'exons' are 1-6, 7-16, 17-36
2331      * CDS 1:12 is transcript 10-21
2332      * transcript 'CDS' is 10-16, 17-21
2333      * which is 'gene' 158-164, 210-214
2334      */
2335     MapList toMap = toLoci.getMapping();
2336     assertEquals(1, toMap.getFromRanges().size());
2337     assertEquals(2, toMap.getFromRanges().get(0).length);
2338     assertEquals(1, toMap.getFromRanges().get(0)[0]);
2339     assertEquals(12, toMap.getFromRanges().get(0)[1]);
2340     assertEquals(2, toMap.getToRanges().size());
2341     assertEquals(2, toMap.getToRanges().get(0).length);
2342     assertEquals(158, toMap.getToRanges().get(0)[0]);
2343     assertEquals(164, toMap.getToRanges().get(0)[1]);
2344     assertEquals(210, toMap.getToRanges().get(1)[0]);
2345     assertEquals(214, toMap.getToRanges().get(1)[1]);
2346     // or summarised as (but toString might change in future):
2347     assertEquals("[ [1, 12] ] 1:1 to [ [158, 164] [210, 214] ]",
2348             toMap.toString());
2349
2350     /*
2351      * an existing value is not overridden
2352      */
2353     geneMap = new MapList(new int[] { 1, 36 }, new int[] { 36, 1 }, 1, 1);
2354     from.setGeneLoci("inhuman", "GRCh37", "6", geneMap);
2355     AlignmentUtils.transferGeneLoci(from, map, to);
2356     assertEquals("GRCh38", toLoci.getAssemblyId());
2357     assertEquals("7", toLoci.getChromosomeId());
2358     toMap = toLoci.getMapping();
2359     assertEquals("[ [1, 12] ] 1:1 to [ [158, 164] [210, 214] ]",
2360             toMap.toString());
2361   }
2362
2363   /**
2364    * Tests for the method that maps nucleotide to protein based on CDS features
2365    */
2366   @Test(groups = "Functional")
2367   public void testMapCdsToProtein()
2368   {
2369     SequenceI peptide = new Sequence("pep", "KLQ");
2370
2371     /*
2372      * Case 1: CDS 3 times length of peptide
2373      * NB method only checks lengths match, not translation
2374      */
2375     SequenceI dna = new Sequence("dna", "AACGacgtCTCCT");
2376     dna.createDatasetSequence();
2377     dna.addSequenceFeature(new SequenceFeature("CDS", "", 1, 4, null));
2378     dna.addSequenceFeature(new SequenceFeature("CDS", "", 9, 13, null));
2379     MapList ml = AlignmentUtils.mapCdsToProtein(dna, peptide);
2380     assertEquals(3, ml.getFromRatio());
2381     assertEquals(1, ml.getToRatio());
2382     assertEquals("[[1, 3]]",
2383             Arrays.deepToString(ml.getToRanges().toArray()));
2384     assertEquals("[[1, 4], [9, 13]]",
2385             Arrays.deepToString(ml.getFromRanges().toArray()));
2386
2387     /*
2388      * Case 2: CDS 3 times length of peptide + stop codon
2389      * (note code does not currently check trailing codon is a stop codon)
2390      */
2391     dna = new Sequence("dna", "AACGacgtCTCCTCCC");
2392     dna.createDatasetSequence();
2393     dna.addSequenceFeature(new SequenceFeature("CDS", "", 1, 4, null));
2394     dna.addSequenceFeature(new SequenceFeature("CDS", "", 9, 16, null));
2395     ml = AlignmentUtils.mapCdsToProtein(dna, peptide);
2396     assertEquals(3, ml.getFromRatio());
2397     assertEquals(1, ml.getToRatio());
2398     assertEquals("[[1, 3]]",
2399             Arrays.deepToString(ml.getToRanges().toArray()));
2400     assertEquals("[[1, 4], [9, 13]]",
2401             Arrays.deepToString(ml.getFromRanges().toArray()));
2402
2403     /*
2404      * Case 3: CDS longer than 3 * peptide + stop codon - no mapping is made
2405      */
2406     dna = new Sequence("dna", "AACGacgtCTCCTTGATCA");
2407     dna.createDatasetSequence();
2408     dna.addSequenceFeature(new SequenceFeature("CDS", "", 1, 4, null));
2409     dna.addSequenceFeature(new SequenceFeature("CDS", "", 9, 19, null));
2410     ml = AlignmentUtils.mapCdsToProtein(dna, peptide);
2411     assertNull(ml);
2412
2413     /*
2414      * Case 4: CDS shorter than 3 * peptide - no mapping is made
2415      */
2416     dna = new Sequence("dna", "AACGacgtCTCC");
2417     dna.createDatasetSequence();
2418     dna.addSequenceFeature(new SequenceFeature("CDS", "", 1, 4, null));
2419     dna.addSequenceFeature(new SequenceFeature("CDS", "", 9, 12, null));
2420     ml = AlignmentUtils.mapCdsToProtein(dna, peptide);
2421     assertNull(ml);
2422
2423     /*
2424      * Case 5: CDS 3 times length of peptide + part codon - mapping is truncated
2425      */
2426     dna = new Sequence("dna", "AACGacgtCTCCTTG");
2427     dna.createDatasetSequence();
2428     dna.addSequenceFeature(new SequenceFeature("CDS", "", 1, 4, null));
2429     dna.addSequenceFeature(new SequenceFeature("CDS", "", 9, 15, null));
2430     ml = AlignmentUtils.mapCdsToProtein(dna, peptide);
2431     assertEquals(3, ml.getFromRatio());
2432     assertEquals(1, ml.getToRatio());
2433     assertEquals("[[1, 3]]",
2434             Arrays.deepToString(ml.getToRanges().toArray()));
2435     assertEquals("[[1, 4], [9, 13]]",
2436             Arrays.deepToString(ml.getFromRanges().toArray()));
2437
2438     /*
2439      * Case 6: incomplete start codon corresponding to X in peptide
2440      */
2441     dna = new Sequence("dna", "ACGacgtCTCCTTGG");
2442     dna.createDatasetSequence();
2443     SequenceFeature sf = new SequenceFeature("CDS", "", 1, 3, null);
2444     sf.setPhase("2"); // skip 2 positions (AC) to start of next codon (GCT)
2445     dna.addSequenceFeature(sf);
2446     dna.addSequenceFeature(new SequenceFeature("CDS", "", 8, 15, null));
2447     peptide = new Sequence("pep", "XLQ");
2448     ml = AlignmentUtils.mapCdsToProtein(dna, peptide);
2449     assertEquals("[[2, 3]]",
2450             Arrays.deepToString(ml.getToRanges().toArray()));
2451     assertEquals("[[3, 3], [8, 12]]",
2452             Arrays.deepToString(ml.getFromRanges().toArray()));
2453   }
2454
2455   /**
2456    * Tests for the method that locates the CDS sequence that has a mapping to
2457    * the given protein. That is, given a transcript-to-peptide mapping, find the
2458    * cds-to-peptide mapping that relates to both, and return the CDS sequence.
2459    */
2460   @Test
2461   public void testFindCdsForProtein()
2462   {
2463     List<AlignedCodonFrame> mappings = new ArrayList<>();
2464     AlignedCodonFrame acf1 = new AlignedCodonFrame();
2465     mappings.add(acf1);
2466
2467     SequenceI dna1 = new Sequence("dna1", "cgatATcgGCTATCTATGacg");
2468     dna1.createDatasetSequence();
2469
2470     // NB we currently exclude STOP codon from CDS sequences
2471     // the test would need to change if this changes in future
2472     SequenceI cds1 = new Sequence("cds1", "ATGCTATCT");
2473     cds1.createDatasetSequence();
2474
2475     SequenceI pep1 = new Sequence("pep1", "MLS");
2476     pep1.createDatasetSequence();
2477     List<AlignedCodonFrame> seqMappings = new ArrayList<>();
2478     MapList mapList = new MapList(new int[] { 5, 6, 9, 15 },
2479             new int[]
2480             { 1, 3 }, 3, 1);
2481     Mapping dnaToPeptide = new Mapping(pep1.getDatasetSequence(), mapList);
2482
2483     // add dna to peptide mapping
2484     seqMappings.add(acf1);
2485     acf1.addMap(dna1.getDatasetSequence(), pep1.getDatasetSequence(),
2486             mapList);
2487
2488     /*
2489      * first case - no dna-to-CDS mapping exists - search fails
2490      */
2491     SequenceI seq = AlignmentUtils.findCdsForProtein(mappings, dna1,
2492             seqMappings, dnaToPeptide);
2493     assertNull(seq);
2494
2495     /*
2496      * second case - CDS-to-peptide mapping exists but no dna-to-CDS
2497      * - search fails
2498      */
2499     // todo this test fails if the mapping is added to acf1, not acf2
2500     // need to tidy up use of lists of mappings in AlignedCodonFrame
2501     AlignedCodonFrame acf2 = new AlignedCodonFrame();
2502     mappings.add(acf2);
2503     MapList cdsToPeptideMapping = new MapList(new int[] { 1, 9 },
2504             new int[]
2505             { 1, 3 }, 3, 1);
2506     acf2.addMap(cds1.getDatasetSequence(), pep1.getDatasetSequence(),
2507             cdsToPeptideMapping);
2508     assertNull(AlignmentUtils.findCdsForProtein(mappings, dna1, seqMappings,
2509             dnaToPeptide));
2510
2511     /*
2512      * third case - add dna-to-CDS mapping - CDS is now found!
2513      */
2514     MapList dnaToCdsMapping = new MapList(new int[] { 5, 6, 9, 15 },
2515             new int[]
2516             { 1, 9 }, 1, 1);
2517     acf1.addMap(dna1.getDatasetSequence(), cds1.getDatasetSequence(),
2518             dnaToCdsMapping);
2519     seq = AlignmentUtils.findCdsForProtein(mappings, dna1, seqMappings,
2520             dnaToPeptide);
2521     assertSame(seq, cds1.getDatasetSequence());
2522   }
2523
2524   /**
2525    * Tests for the method that locates the CDS sequence that has a mapping to
2526    * the given protein. That is, given a transcript-to-peptide mapping, find the
2527    * cds-to-peptide mapping that relates to both, and return the CDS sequence.
2528    * This test is for the case where transcript and CDS are the same length.
2529    */
2530   @Test
2531   public void testFindCdsForProtein_noUTR()
2532   {
2533     List<AlignedCodonFrame> mappings = new ArrayList<>();
2534     AlignedCodonFrame acf1 = new AlignedCodonFrame();
2535     mappings.add(acf1);
2536
2537     SequenceI dna1 = new Sequence("dna1", "ATGCTATCTTAA");
2538     dna1.createDatasetSequence();
2539
2540     // NB we currently exclude STOP codon from CDS sequences
2541     // the test would need to change if this changes in future
2542     SequenceI cds1 = new Sequence("cds1", "ATGCTATCT");
2543     cds1.createDatasetSequence();
2544
2545     SequenceI pep1 = new Sequence("pep1", "MLS");
2546     pep1.createDatasetSequence();
2547     List<AlignedCodonFrame> seqMappings = new ArrayList<>();
2548     MapList mapList = new MapList(new int[] { 1, 9 }, new int[] { 1, 3 }, 3,
2549             1);
2550     Mapping dnaToPeptide = new Mapping(pep1.getDatasetSequence(), mapList);
2551
2552     // add dna to peptide mapping
2553     seqMappings.add(acf1);
2554     acf1.addMap(dna1.getDatasetSequence(), pep1.getDatasetSequence(),
2555             mapList);
2556
2557     /*
2558      * first case - transcript lacks CDS features - it appears to be
2559      * the CDS sequence and is returned
2560      */
2561     SequenceI seq = AlignmentUtils.findCdsForProtein(mappings, dna1,
2562             seqMappings, dnaToPeptide);
2563     assertSame(seq, dna1.getDatasetSequence());
2564
2565     /*
2566      * second case - transcript has CDS feature - this means it is
2567      * not returned as a match for CDS (CDS sequences don't have CDS features)
2568      */
2569     dna1.addSequenceFeature(
2570             new SequenceFeature(SequenceOntologyI.CDS, "cds", 1, 12, null));
2571     seq = AlignmentUtils.findCdsForProtein(mappings, dna1, seqMappings,
2572             dnaToPeptide);
2573     assertNull(seq);
2574
2575     /*
2576      * third case - CDS-to-peptide mapping exists but no dna-to-CDS
2577      * - search fails
2578      */
2579     // todo this test fails if the mapping is added to acf1, not acf2
2580     // need to tidy up use of lists of mappings in AlignedCodonFrame
2581     AlignedCodonFrame acf2 = new AlignedCodonFrame();
2582     mappings.add(acf2);
2583     MapList cdsToPeptideMapping = new MapList(new int[] { 1, 9 },
2584             new int[]
2585             { 1, 3 }, 3, 1);
2586     acf2.addMap(cds1.getDatasetSequence(), pep1.getDatasetSequence(),
2587             cdsToPeptideMapping);
2588     assertNull(AlignmentUtils.findCdsForProtein(mappings, dna1, seqMappings,
2589             dnaToPeptide));
2590
2591     /*
2592      * fourth case - add dna-to-CDS mapping - CDS is now found!
2593      */
2594     MapList dnaToCdsMapping = new MapList(new int[] { 1, 9 },
2595             new int[]
2596             { 1, 9 }, 1, 1);
2597     acf1.addMap(dna1.getDatasetSequence(), cds1.getDatasetSequence(),
2598             dnaToCdsMapping);
2599     seq = AlignmentUtils.findCdsForProtein(mappings, dna1, seqMappings,
2600             dnaToPeptide);
2601     assertSame(seq, cds1.getDatasetSequence());
2602   }
2603 }