test/jalview/analysis/AlignmentUtilsTests.java

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