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