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