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