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