test/jalview/io/vcf/VCFLoaderTest.java

   1 package jalview.io.vcf;
   2
   3 import static org.testng.Assert.assertEquals;
   4
   5 import jalview.datamodel.AlignmentI;
   6 import jalview.datamodel.DBRefEntry;
   7 import jalview.datamodel.GeneLoci;
   8 import jalview.datamodel.Mapping;
   9 import jalview.datamodel.Sequence;
  10 import jalview.datamodel.SequenceFeature;
  11 import jalview.datamodel.SequenceI;
  12 import jalview.gui.AlignFrame;
  13 import jalview.io.DataSourceType;
  14 import jalview.io.FileLoader;
  15 import jalview.io.gff.Gff3Helper;
  16 import jalview.io.gff.SequenceOntologyI;
  17 import jalview.util.MapList;
  18
  19 import java.io.File;
  20 import java.io.IOException;
  21 import java.io.PrintWriter;
  22 import java.util.Arrays;
  23 import java.util.List;
  24
  25 import org.testng.annotations.Test;
  26
  27 public class VCFLoaderTest
  28 {
  29   // columns 9717- of gene P30419 from Ensembl (modified)
  30   private static final String FASTA = ">ENSG00000136448/1-25 chromosome:GRCh38:17:45051610:45051634:1\n"
  31           + "CAAGCTGGCGGACGAGAGTGTGACA\n"
  32           // and a 'made up' mini-transcript with two exons
  33           + ">ENST00000592782/1-18\n--AGCTGGCG----AGAGTGTGAC-\n";
  34
  35   private static final String[] VCF = { "##fileformat=VCFv4.2",
  36       "##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency, for each ALT allele, in the same order as listed\">",
  37       "##reference=GRCh38",
  38       "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO",
  39       // SNP A/T in position 2 of gene sequence (precedes transcript)
  40       "17\t45051611\t.\tA\tT\t1666.64\tRF\tAC=15;AF=5.08130e-03",
  41       // SNP G/C in position 4 of gene sequence, position 2 of transcript
  42       // this is a mixed variant, the insertion G/GA is not transferred
  43       "17\t45051613\t.\tG\tGA,C\t1666.64\tRF\tAC=15;AF=3.08130e-03" };
  44
  45   @Test(groups = "Functional")
  46   public void testLoadVCF() throws IOException
  47   {
  48     AlignmentI al = buildAlignment();
  49     VCFLoader loader = new VCFLoader(al);
  50
  51     File f = makeVcf();
  52
  53     loader.loadVCF(f.getPath(), null);
  54
  55     /*
  56      * verify variant feature(s) added to gene
  57      */
  58     List<SequenceFeature> geneFeatures = al.getSequenceAt(0)
  59             .getSequenceFeatures();
  60     assertEquals(geneFeatures.size(), 2);
  61     SequenceFeature sf = geneFeatures.get(0);
  62     assertEquals(sf.getFeatureGroup(), "VCF");
  63     assertEquals(sf.getBegin(), 2);
  64     assertEquals(sf.getEnd(), 2);
  65     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
  66     assertEquals(sf.getScore(), 5.08130e-03, 0.000001f);
  67     assertEquals("A,T", sf.getValue(Gff3Helper.ALLELES));
  68
  69     sf = geneFeatures.get(1);
  70     assertEquals(sf.getFeatureGroup(), "VCF");
  71     assertEquals(sf.getBegin(), 4);
  72     assertEquals(sf.getEnd(), 4);
  73     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
  74     assertEquals(sf.getScore(), 3.08130e-03, 0.000001f);
  75     assertEquals("G,C", sf.getValue(Gff3Helper.ALLELES));
  76
  77     /*
  78      * verify variant feature(s) added to transcript
  79      */
  80     List<SequenceFeature> transcriptFeatures = al.getSequenceAt(1)
  81             .getSequenceFeatures();
  82     assertEquals(transcriptFeatures.size(), 1);
  83     sf = transcriptFeatures.get(0);
  84     assertEquals(sf.getFeatureGroup(), "VCF");
  85     assertEquals(sf.getBegin(), 2);
  86     assertEquals(sf.getEnd(), 2);
  87     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
  88     assertEquals(sf.getScore(), 3.08130e-03, 0.000001f);
  89     assertEquals("G,C", sf.getValue(Gff3Helper.ALLELES));
  90
  91     /*
  92      * verify variant feature(s) computed and added to protein
  93      * first codon AGC varies to ACC giving S/T
  94      */
  95     SequenceI peptide = al.getSequenceAt(1)
  96             .getDBRefs()[0].getMap().getTo();
  97     List<SequenceFeature> proteinFeatures = peptide.getSequenceFeatures();
  98     assertEquals(proteinFeatures.size(), 1);
  99     sf = proteinFeatures.get(0);
 100     assertEquals(sf.getFeatureGroup(), "VCF");
 101     assertEquals(sf.getBegin(), 1);
 102     assertEquals(sf.getEnd(), 1);
 103     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
 104     assertEquals(sf.getDescription(), "p.Ser1Thr");
 105   }
 106
 107   private File makeVcf() throws IOException
 108   {
 109     File f = File.createTempFile("Test", ".vcf");
 110     f.deleteOnExit();
 111     PrintWriter pw = new PrintWriter(f);
 112     for (String vcfLine : VCF)
 113     {
 114       pw.println(vcfLine);
 115     }
 116     pw.close();
 117     return f;
 118   }
 119
 120   /**
 121    * Make a simple alignment with one 'gene' and one 'transcript'
 122    *
 123    * @return
 124    */
 125   private AlignmentI buildAlignment()
 126   {
 127     AlignFrame af = new FileLoader().LoadFileWaitTillLoaded(FASTA,
 128             DataSourceType.PASTE);
 129
 130     /*
 131      * map gene sequence to chromosome (normally done when the sequence is fetched
 132      * from Ensembl and transcripts computed)
 133      */
 134     AlignmentI alignment = af.getViewport().getAlignment();
 135     int[][] to = new int[][] { new int[] { 45051610, 45051634 } };
 136     List<int[]> toRanges = Arrays.asList(to);
 137     SequenceI gene = alignment.getSequenceAt(0);
 138     List<int[]> fromRanges = Arrays.asList(new int[][] { new int[] {
 139         gene.getStart(), gene.getEnd() } });
 140     ((Sequence) gene).setGeneLoci(new GeneLoci("human", "GRCh38", "17",
 141             new MapList(fromRanges, toRanges, 1, 1)));
 142
 143     /*
 144      * map 'transcript' to chromosome via 'gene'
 145      * transcript/1-18 is gene/3-10,15-24
 146      * which is chromosome 45051612-45051619,45051624-45051633
 147      */
 148     to = new int[][] { new int[] { 45051612, 45051619 },
 149         new int[] { 45051624, 45051633 } };
 150     toRanges = Arrays.asList(to);
 151     SequenceI transcript = alignment.getSequenceAt(1);
 152     fromRanges = Arrays.asList(new int[][] { new int[] {
 153         transcript.getStart(), transcript.getEnd() } });
 154     ((Sequence) transcript).setGeneLoci(new GeneLoci("human", "GRCh38",
 155             "17", new MapList(fromRanges, toRanges, 1, 1)));
 156
 157     /*
 158      * add a protein product as a DBRef on the transcript
 159      */
 160     SequenceI peptide = new Sequence("ENSP001", "SWRECD");
 161     MapList mapList = new MapList(new int[] { 1, 18 }, new int[] { 1, 6 },
 162             3, 1);
 163     Mapping map = new Mapping(peptide, mapList);
 164     DBRefEntry product = new DBRefEntry("", "", "ENSP001", map);
 165     transcript.addDBRef(product);
 166
 167     return alignment;
 168   }
 169
 170   /**
 171    * Test with 'gene' and 'transcript' mapped to the reverse strand of the
 172    * chromosome. The VCF variant positions (in forward coordinates) should get
 173    * correctly located on sequence positions.
 174    *
 175    * @throws IOException
 176    */
 177   @Test(groups = "Functional")
 178   public void testLoadVCF_reverseStrand() throws IOException
 179   {
 180     AlignmentI al = buildAlignment();
 181
 182     /*
 183      * invert forward to reverse strand mappings
 184      */
 185     List<int[]> to = al.getSequenceAt(0).getGeneLoci().mapping
 186             .getToRanges();
 187     int temp = to.get(0)[0];
 188     to.get(0)[0] = to.get(0)[1];
 189     to.get(0)[1] = temp;
 190     to = al.getSequenceAt(1).getGeneLoci().mapping.getToRanges();
 191     to.get(0)[0] = to.get(0)[1];
 192     to.get(0)[1] = temp;
 193     to.get(1)[0] = to.get(1)[1];
 194     to.get(1)[1] = temp;
 195     int[] tmp2 = to.get(0);
 196     to.set(0, to.get(1));
 197     to.set(1, tmp2);
 198
 199     VCFLoader loader = new VCFLoader(al);
 200
 201     File f = makeVcf();
 202
 203     loader.loadVCF(f.getPath(), null);
 204
 205     /*
 206      * verify variant feature(s) added to gene
 207      */
 208     List<SequenceFeature> geneFeatures = al.getSequenceAt(0)
 209             .getSequenceFeatures();
 210     assertEquals(geneFeatures.size(), 2);
 211     SequenceFeature sf = geneFeatures.get(0);
 212     assertEquals(sf.getFeatureGroup(), "VCF");
 213     assertEquals(sf.getBegin(), 2);
 214     assertEquals(sf.getEnd(), 2);
 215     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
 216     assertEquals(sf.getScore(), 5.08130e-03, 0.000001f);
 217     assertEquals("A,T", sf.getValue(Gff3Helper.ALLELES));
 218
 219     sf = geneFeatures.get(1);
 220     assertEquals(sf.getFeatureGroup(), "VCF");
 221     assertEquals(sf.getBegin(), 4);
 222     assertEquals(sf.getEnd(), 4);
 223     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
 224     assertEquals(sf.getScore(), 3.08130e-03, 0.000001f);
 225     assertEquals("G,C", sf.getValue(Gff3Helper.ALLELES));
 226
 227     /*
 228      * verify variant feature(s) added to transcript
 229      */
 230     List<SequenceFeature> transcriptFeatures = al.getSequenceAt(1)
 231             .getSequenceFeatures();
 232     assertEquals(transcriptFeatures.size(), 1);
 233     sf = transcriptFeatures.get(0);
 234     assertEquals(sf.getFeatureGroup(), "VCF");
 235     assertEquals(sf.getBegin(), 2);
 236     assertEquals(sf.getEnd(), 2);
 237     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
 238     assertEquals(sf.getScore(), 3.08130e-03, 0.000001f);
 239     assertEquals("G,C", sf.getValue(Gff3Helper.ALLELES));
 240
 241     /*
 242      * verify variant feature(s) computed and added to protein
 243      * first codon AGC varies to ACC giving S/T
 244      */
 245     SequenceI peptide = al.getSequenceAt(1).getDBRefs()[0].getMap().getTo();
 246     List<SequenceFeature> proteinFeatures = peptide.getSequenceFeatures();
 247     assertEquals(proteinFeatures.size(), 1);
 248     sf = proteinFeatures.get(0);
 249     assertEquals(sf.getFeatureGroup(), "VCF");
 250     assertEquals(sf.getBegin(), 1);
 251     assertEquals(sf.getEnd(), 1);
 252     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
 253     assertEquals(sf.getDescription(), "p.Ser1Thr");
 254   }
 255 }