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