test/jalview/analysis/AlignmentUtilsTests.java

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