test/jalview/analysis/AlignmentUtilsTests.java

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