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