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