test/jalview/analysis/AlignmentUtilsTests.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.2)
   3  * Copyright (C) 2014 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.junit.Assert.assertEquals;
  24 import static org.junit.Assert.assertSame;
  25 import static org.junit.Assert.assertTrue;
  26 import jalview.analysis.AlignmentUtils.MappingResult;
  27 import jalview.datamodel.AlignedCodonFrame;
  28 import jalview.datamodel.Alignment;
  29 import jalview.datamodel.AlignmentI;
  30 import jalview.datamodel.Mapping;
  31 import jalview.datamodel.Sequence;
  32 import jalview.datamodel.SequenceI;
  33 import jalview.io.AppletFormatAdapter;
  34 import jalview.io.FormatAdapter;
  35 import jalview.util.MapList;
  36
  37 import java.io.IOException;
  38 import java.util.Arrays;
  39 import java.util.List;
  40 import java.util.Map;
  41
  42 import org.junit.Test;
  43
  44 public class AlignmentUtilsTests
  45 {
  46   // @formatter:off
  47   private static final String TEST_DATA =
  48           "# STOCKHOLM 1.0\n" +
  49           "#=GS D.melanogaster.1 AC AY119185.1/838-902\n" +
  50           "#=GS D.melanogaster.2 AC AC092237.1/57223-57161\n" +
  51           "#=GS D.melanogaster.3 AC AY060611.1/560-627\n" +
  52           "D.melanogaster.1          G.AGCC.CU...AUGAUCGA\n" +
  53           "#=GR D.melanogaster.1 SS  ................((((\n" +
  54           "D.melanogaster.2          C.AUUCAACU.UAUGAGGAU\n" +
  55           "#=GR D.melanogaster.2 SS  ................((((\n" +
  56           "D.melanogaster.3          G.UGGCGCU..UAUGACGCA\n" +
  57           "#=GR D.melanogaster.3 SS  (.(((...(....(((((((\n" +
  58           "//";
  59
  60   private static final String AA_SEQS_1 =
  61           ">Seq1Name\n" +
  62           "K-QY--L\n" +
  63           ">Seq2Name\n" +
  64           "-R-FP-W-\n";
  65
  66   private static final String CDNA_SEQS_1 =
  67           ">Seq1Name\n" +
  68           "AC-GG--CUC-CAA-CT\n" +
  69           ">Seq2Name\n" +
  70           "-CG-TTA--ACG---AAGT\n";
  71
  72   private static final String CDNA_SEQS_2 =
  73           ">Seq1Name\n" +
  74           "GCTCGUCGTACT\n" +
  75           ">Seq2Name\n" +
  76           "GGGTCAGGCAGT\n";
  77   // @formatter:on
  78
  79   public static Sequence ts=new Sequence("short","ASDASDASDASDASDASDASDASDASDASDASDASDASD");
  80
  81   @Test
  82   public void testExpandFlanks()
  83   {
  84     AlignmentI al = new Alignment(new Sequence[] {});
  85     for (int i=4;i<14;i+=3)
  86     {
  87       SequenceI s1=ts.deriveSequence().getSubSequence(i, i+7);
  88       al.addSequence(s1);
  89     }
  90     System.out.println(new AppletFormatAdapter().formatSequences("Clustal", al, true));
  91     for (int flnk=-1;flnk<25; flnk++)
  92     {
  93       AlignmentI exp;
  94       System.out.println("\nFlank size: "+flnk);
  95       System.out.println(new AppletFormatAdapter().formatSequences("Clustal", exp=AlignmentUtils.expandContext(al, flnk), true));
  96       if (flnk==-1) {
  97         for (SequenceI sq:exp.getSequences())
  98       {
  99           String ung = sq.getSequenceAsString().replaceAll("-+", "");
 100           assertTrue("Flanking sequence not the same as original dataset sequence.\n"+ung+"\n"+sq.getDatasetSequence().getSequenceAsString(),ung.equalsIgnoreCase(sq.getDatasetSequence().getSequenceAsString()));
 101       }
 102       }
 103     }
 104     }
 105
 106   /**
 107    * Test method that returns a map of lists of sequences by sequence name.
 108    *
 109    * @throws IOException
 110    */
 111   @Test
 112   public void testGetSequencesByName() throws IOException
 113   {
 114     final String data = ">Seq1Name\nKQYL\n" + ">Seq2Name\nRFPW\n"
 115             + ">Seq1Name\nABCD\n";
 116     AlignmentI al = loadAlignment(data, "FASTA");
 117     Map<String, List<SequenceI>> map = AlignmentUtils
 118             .getSequencesByName(al);
 119     assertEquals(2, map.keySet().size());
 120     assertEquals(2, map.get("Seq1Name").size());
 121     assertEquals("KQYL", map.get("Seq1Name").get(0).getSequenceAsString());
 122     assertEquals("ABCD", map.get("Seq1Name").get(1).getSequenceAsString());
 123     assertEquals(1, map.get("Seq2Name").size());
 124     assertEquals("RFPW", map.get("Seq2Name").get(0).getSequenceAsString());
 125   }
 126   /**
 127    * Helper method to load an alignment and ensure dataset sequences are set up.
 128    *
 129    * @param data
 130    * @param format TODO
 131    * @return
 132    * @throws IOException
 133    */
 134   protected AlignmentI loadAlignment(final String data, String format) throws IOException
 135   {
 136     Alignment a = new FormatAdapter().readFile(data,
 137             AppletFormatAdapter.PASTE, format);
 138     a.setDataset(null);
 139     return a;
 140   }
 141   /**
 142    * Test mapping of protein to cDNA.
 143    *
 144    * @throws IOException
 145    */
 146   @Test
 147   public void testMapProteinToCdna() throws IOException
 148   {
 149     // protein: Human + Mouse, 3 residues
 150     AlignmentI protein = loadAlignment(
 151             ">Human\nKQY\n>Mouse\nAFP\n>Worm\nRST\n",
 152             "FASTA");
 153     // cDNA: Mouse, Human, Mouse, 9 bases
 154     // @formatter:off
 155     String dnaData =
 156             ">Mouse\nGAAATCCAG\n" +
 157             ">Human\nTTCGATTAC\n" +
 158             ">Mouse\nGTCGTTTGC\n" +
 159             ">Mouse\nGTCGTTTGCgac\n" + // not mapped - wrong length
 160             ">Fly\nGTCGTTTGC\n"; // not mapped - no name match
 161     // @formatter:on
 162     AlignmentI cdna1 = loadAlignment(
 163             dnaData,
 164             "FASTA");
 165     MappingResult mapped = AlignmentUtils.mapProteinToCdna(protein, cdna1);
 166     assertEquals(mapped, MappingResult.Mapped);
 167
 168     /*
 169      * Check two mappings (one for Mouse, one for Human)
 170      */
 171     assertEquals(2, protein.getCodonFrames().size());
 172     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(0)).size());
 173     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(1)).size());
 174
 175     /*
 176      * Inspect mapping for Human protein
 177      */
 178     AlignedCodonFrame humanMapping = protein.getCodonFrame(
 179             protein.getSequenceAt(0)).get(0);
 180     assertEquals(1, humanMapping.getdnaSeqs().length);
 181     assertEquals(cdna1.getSequenceAt(1).getDatasetSequence(),
 182             humanMapping.getdnaSeqs()[0]);
 183     Mapping[] protMappings = humanMapping.getProtMappings();
 184     assertEquals(1, protMappings.length);
 185     MapList mapList = protMappings[0].getMap();
 186     assertEquals(3, mapList.getFromRatio());
 187     assertEquals(1, mapList.getToRatio());
 188     assertTrue(Arrays.equals(new int[]
 189     { 1, 9 }, mapList.getFromRanges()));
 190     assertTrue(Arrays.equals(new int[]
 191     { 1, 3 }, mapList.getToRanges()));
 192
 193     /*
 194      * Inspect mappings for Mouse protein
 195      */
 196     AlignedCodonFrame mouseMapping1 = protein.getCodonFrame(
 197             protein.getSequenceAt(1)).get(0);
 198     assertEquals(2, mouseMapping1.getdnaSeqs().length);
 199     assertEquals(cdna1.getSequenceAt(0).getDatasetSequence(),
 200             mouseMapping1.getdnaSeqs()[0]);
 201     assertEquals(cdna1.getSequenceAt(2).getDatasetSequence(),
 202             mouseMapping1.getdnaSeqs()[1]);
 203     protMappings = mouseMapping1.getProtMappings();
 204     assertEquals(2, protMappings.length);
 205     for (int i = 0; i < 2; i++)
 206     {
 207       mapList = protMappings[i].getMap();
 208       assertEquals(3, mapList.getFromRatio());
 209       assertEquals(1, mapList.getToRatio());
 210       assertTrue(Arrays.equals(new int[]
 211       { 1, 9 }, mapList.getFromRanges()));
 212       assertTrue(Arrays.equals(new int[]
 213       { 1, 3 }, mapList.getToRanges()));
 214     }
 215   }
 216
 217   /**
 218    * Test mapping of protein to cDNA which may include start and/or stop codons.
 219    *
 220    * @throws IOException
 221    */
 222   @Test
 223   public void testMapProteinToCdna_stopStartCodons() throws IOException
 224   {
 225     // protein: Human + Mouse, 3 residues
 226     AlignmentI protein = loadAlignment(
 227             ">Human\nKQY\n>Mouse\nAFP\n>Worm\nRST\n", "FASTA");
 228     // @formatter:off
 229     String dnaData =
 230             ">Mouse\natgGAAATCCAG\n" + // Mouse with start codon
 231             ">Human\nTTCGATtactaa\n" + // Human with stop codon TAA
 232             ">Mouse\nGTCGTTTGctaG\n" + // Mouse with stop codon TAG
 233             ">Human\nGTCGTTTgctGa\n" + // Human with stop codon TGA
 234             ">Mouse\nATGGTCGTTTGCtag\n"; // Mouse with start and stop codons
 235     // @formatter:on
 236     AlignmentI cdna1 = loadAlignment(
 237             dnaData,
 238             "FASTA");
 239     MappingResult mapped = AlignmentUtils.mapProteinToCdna(protein, cdna1);
 240     assertEquals(mapped, MappingResult.Mapped);
 241
 242     /*
 243      * Check two mappings (one for Mouse, one for Human)
 244      */
 245     assertEquals(2, protein.getCodonFrames().size());
 246     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(0)).size());
 247     assertEquals(1, protein.getCodonFrame(protein.getSequenceAt(1)).size());
 248
 249     /*
 250      * Inspect mapping for Human protein - should map to 2nd and 4th cDNA seqs
 251      */
 252     AlignedCodonFrame humanMapping = protein.getCodonFrame(
 253             protein.getSequenceAt(0)).get(0);
 254     assertEquals(2, humanMapping.getdnaSeqs().length);
 255     assertEquals(cdna1.getSequenceAt(1).getDatasetSequence(),
 256             humanMapping.getdnaSeqs()[0]);
 257     assertEquals(cdna1.getSequenceAt(3).getDatasetSequence(),
 258             humanMapping.getdnaSeqs()[1]);
 259     Mapping[] protMappings = humanMapping.getProtMappings();
 260     // two mappings, both to cDNA with stop codon
 261     assertEquals(2, protMappings.length);
 262     MapList mapList = protMappings[0].getMap();
 263     assertEquals(3, mapList.getFromRatio());
 264     assertEquals(1, mapList.getToRatio());
 265     assertTrue(Arrays.equals(new int[]
 266     { 1, 9 }, mapList.getFromRanges()));
 267     assertTrue(Arrays.equals(new int[]
 268     { 1, 3 }, mapList.getToRanges()));
 269     mapList = protMappings[1].getMap();
 270     assertEquals(3, mapList.getFromRatio());
 271     assertEquals(1, mapList.getToRatio());
 272     assertTrue(Arrays.equals(new int[]
 273     { 1, 9 }, mapList.getFromRanges()));
 274     assertTrue(Arrays.equals(new int[]
 275     { 1, 3 }, mapList.getToRanges()));
 276
 277     /*
 278      * Inspect mapping for Mouse protein - should map to 1st/3rd/5th cDNA seqs
 279      */
 280     AlignedCodonFrame mouseMapping = protein.getCodonFrame(
 281             protein.getSequenceAt(1)).get(0);
 282     assertEquals(3, mouseMapping.getdnaSeqs().length);
 283     assertEquals(cdna1.getSequenceAt(0).getDatasetSequence(),
 284             mouseMapping.getdnaSeqs()[0]);
 285     assertEquals(cdna1.getSequenceAt(2).getDatasetSequence(),
 286             mouseMapping.getdnaSeqs()[1]);
 287     assertEquals(cdna1.getSequenceAt(4).getDatasetSequence(),
 288             mouseMapping.getdnaSeqs()[2]);
 289
 290     // three mappings
 291     protMappings = mouseMapping.getProtMappings();
 292     assertEquals(3, protMappings.length);
 293
 294     // first mapping to cDNA with start codon
 295     mapList = protMappings[0].getMap();
 296     assertEquals(3, mapList.getFromRatio());
 297     assertEquals(1, mapList.getToRatio());
 298     assertTrue(Arrays.equals(new int[]
 299     { 4, 12 }, mapList.getFromRanges()));
 300     assertTrue(Arrays.equals(new int[]
 301     { 1, 3 }, mapList.getToRanges()));
 302
 303     // second mapping to cDNA with stop codon
 304     mapList = protMappings[1].getMap();
 305     assertEquals(3, mapList.getFromRatio());
 306     assertEquals(1, mapList.getToRatio());
 307     assertTrue(Arrays.equals(new int[]
 308     { 1, 9 }, mapList.getFromRanges()));
 309     assertTrue(Arrays.equals(new int[]
 310     { 1, 3 }, mapList.getToRanges()));
 311
 312     // third mapping to cDNA with start and stop codon
 313     mapList = protMappings[2].getMap();
 314     assertEquals(3, mapList.getFromRatio());
 315     assertEquals(1, mapList.getToRatio());
 316     assertTrue(Arrays.equals(new int[]
 317     { 4, 12 }, mapList.getFromRanges()));
 318     assertTrue(Arrays.equals(new int[]
 319     { 1, 3 }, mapList.getToRanges()));
 320   }
 321
 322   /**
 323    * Test for the alignSequenceAs method that takes two sequences and a mapping.
 324    */
 325   @Test
 326   public void testAlignSequenceAs_withMapping_noIntrons()
 327   {
 328     MapList map = new MapList(new int[]
 329     { 1, 6 }, new int[]
 330     { 1, 2 }, 3, 1);
 331
 332     /*
 333      * No existing gaps in dna:
 334      */
 335     checkAlignSequenceAs("GGGAAA", "-A-L-", false, false, map,
 336             "---GGG---AAA");
 337
 338     /*
 339      * Now introduce gaps in dna but ignore them when realigning.
 340      */
 341     checkAlignSequenceAs("-G-G-G-A-A-A-", "-A-L-", false, false, map,
 342             "---GGG---AAA");
 343
 344     /*
 345      * Now include gaps in dna when realigning. First retaining 'mapped' gaps
 346      * only, i.e. those within the exon region.
 347      */
 348     checkAlignSequenceAs("-G-G--G-A--A-A-", "-A-L-", true, false, map,
 349             "---G-G--G---A--A-A");
 350
 351     /*
 352      * Include all gaps in dna when realigning (within and without the exon
 353      * region). The leading gap, and the gaps between codons, are subsumed by
 354      * the protein alignment gap.
 355      */
 356     checkAlignSequenceAs("-G-GG--AA-A-", "-A-L-", true, true, map,
 357             "---G-GG---AA-A-");
 358
 359     /*
 360      * Include only unmapped gaps in dna when realigning (outside the exon
 361      * region). The leading gap, and the gaps between codons, are subsumed by
 362      * the protein alignment gap.
 363      */
 364     checkAlignSequenceAs("-G-GG--AA-A-", "-A-L-", false, true, map,
 365             "---GGG---AAA-");
 366   }
 367
 368   /**
 369    * Test for the alignSequenceAs method that takes two sequences and a mapping.
 370    */
 371   @Test
 372   public void testAlignSequenceAs_withMapping_withIntrons()
 373   {
 374     /*
 375      * Exons at codon 2 (AAA) and 4 (TTT)
 376      */
 377     MapList map = new MapList(new int[]
 378     { 4, 6, 10, 12 }, new int[]
 379     { 1, 2 }, 3, 1);
 380
 381     /*
 382      * Simple case: no gaps in dna
 383      */
 384     checkAlignSequenceAs("GGGAAACCCTTTGGG", "--A-L-", false, false, map,
 385             "GGG---AAACCCTTTGGG");
 386
 387     /*
 388      * Add gaps to dna - but ignore when realigning.
 389      */
 390     checkAlignSequenceAs("-G-G-G--A--A---AC-CC-T-TT-GG-G-", "--A-L-",
 391             false, false, map, "GGG---AAACCCTTTGGG");
 392
 393     /*
 394      * Add gaps to dna - include within exons only when realigning.
 395      */
 396     checkAlignSequenceAs("-G-G-G--A--A---A-C-CC-T-TT-GG-G-", "--A-L-",
 397             true, false, map, "GGG---A--A---ACCCT-TTGGG");
 398
 399     /*
 400      * Include gaps outside exons only when realigning.
 401      */
 402     checkAlignSequenceAs("-G-G-G--A--A---A-C-CC-T-TT-GG-G-", "--A-L-",
 403             false, true, map, "-G-G-GAAAC-CCTTT-GG-G-");
 404
 405     /*
 406      * Include gaps following first intron if we are 'preserving mapped gaps'
 407      */
 408     checkAlignSequenceAs("-G-G-G--A--A---A-C-CC-T-TT-GG-G-", "--A-L-",
 409             true, true, map, "-G-G-G--A--A---A-C-CC-T-TT-GG-G-");
 410
 411     /*
 412      * Include all gaps in dna when realigning.
 413      */
 414     checkAlignSequenceAs("-G-G-G--A--A---A-C-CC-T-TT-GG-G-", "--A-L-",
 415             true, true, map, "-G-G-G--A--A---A-C-CC-T-TT-GG-G-");
 416   }
 417
 418   /**
 419    * Test for the case where not all of the protein sequence is mapped to cDNA.
 420    */
 421   @Test
 422   public void testAlignSequenceAs_withMapping_withUnmappedProtein()
 423   {
 424
 425     /*
 426      * Exons at codon 2 (AAA) and 4 (TTT) mapped to A and P
 427      */
 428     final MapList map = new MapList(new int[]
 429     { 4, 6, 10, 12 }, new int[]
 430     { 1, 1, 3, 3 }, 3, 1);
 431
 432
 433     /*
 434      * Expect alignment does nothing (aborts realignment). Change this test
 435      * first if different behaviour wanted.
 436      */
 437     checkAlignSequenceAs("GGGAAACCCTTTGGG", "-A-L-P-", false,
 438             false, map, "GGGAAACCCTTTGGG");
 439   }
 440
 441   /**
 442    * Helper method that performs and verifies the method under test.
 443    *
 444    * @param dnaSeq
 445    * @param proteinSeq
 446    * @param preserveMappedGaps
 447    * @param preserveUnmappedGaps
 448    * @param map
 449    * @param expected
 450    */
 451   protected void checkAlignSequenceAs(final String dnaSeq,
 452           final String proteinSeq, final boolean preserveMappedGaps,
 453           final boolean preserveUnmappedGaps, MapList map,
 454           final String expected)
 455   {
 456     SequenceI dna = new Sequence("Seq1", dnaSeq);
 457     dna.createDatasetSequence();
 458     SequenceI protein = new Sequence("Seq1", proteinSeq);
 459     protein.createDatasetSequence();
 460     AlignedCodonFrame acf = new AlignedCodonFrame();
 461     acf.addMap(dna.getDatasetSequence(), protein.getDatasetSequence(), map);
 462
 463     AlignmentUtils.alignSequenceAs(dna, protein, acf, "---", '-',
 464             preserveMappedGaps, preserveUnmappedGaps);
 465     assertEquals(expected, dna.getSequenceAsString());
 466   }
 467
 468   /**
 469    * Test for the alignSequenceAs method where we preserve gaps in introns only.
 470    */
 471   @Test
 472   public void testAlignSequenceAs_keepIntronGapsOnly()
 473   {
 474
 475     /*
 476      * Intron GGGAAA followed by exon CCCTTT
 477      */
 478     MapList map = new MapList(new int[]
 479     { 7, 12 }, new int[]
 480     { 1, 2 }, 3, 1);
 481
 482     checkAlignSequenceAs("GG-G-AA-A-C-CC-T-TT", "AL",
 483             false, true, map, "GG-G-AA-ACCCTTT");
 484   }
 485
 486   /**
 487    * Test for the method that generates an aligned translated sequence from one
 488    * mapping.
 489    */
 490   @Test
 491   public void testGetAlignedTranslation_dnaLikeProtein()
 492   {
 493     // dna alignment will be replaced
 494     SequenceI dna = new Sequence("Seq1", "T-G-CC-A--T-TAC-CAG-");
 495     dna.createDatasetSequence();
 496     // protein alignment will be 'applied' to dna
 497     SequenceI protein = new Sequence("Seq1", "-CH-Y--Q-");
 498     protein.createDatasetSequence();
 499     MapList map = new MapList(new int[]
 500     { 1, 12 }, new int[]
 501     { 1, 4 }, 3, 1);
 502     AlignedCodonFrame acf = new AlignedCodonFrame();
 503     acf.addMap(dna.getDatasetSequence(), protein.getDatasetSequence(), map);
 504
 505     final SequenceI aligned = AlignmentUtils
 506                 .getAlignedTranslation(protein, '-', acf);
 507     assertEquals("---TGCCAT---TAC------CAG---", aligned.getSequenceAsString());
 508     assertSame(aligned.getDatasetSequence(), dna.getDatasetSequence());
 509   }
 510 }