export_test.go

   1 // Copyright (C) The Lightning Authors. All rights reserved.
   2 //
   3 // SPDX-License-Identifier: AGPL-3.0
   4
   5 package lightning
   6
   7 import (
   8         "bytes"
   9         "io/ioutil"
  10         "os"
  11         "os/exec"
  12
  13         "github.com/kshedden/gonpy"
  14         "gopkg.in/check.v1"
  15 )
  16
  17 type exportSuite struct{}
  18
  19 var _ = check.Suite(&exportSuite{})
  20
  21 func (s *exportSuite) TestFastaToHGVS(c *check.C) {
  22         tmpdir := c.MkDir()
  23
  24         err := ioutil.WriteFile(tmpdir+"/chr1-12-100.bed", []byte("chr1\t12\t100\ttest.1\n"), 0644)
  25         c.Check(err, check.IsNil)
  26
  27         var buffer bytes.Buffer
  28         exited := (&importer{}).RunCommand("import", []string{"-local=true", "-tag-library", "testdata/tags", "-output-tiles", "-save-incomplete-tiles", "testdata/pipeline1", "testdata/ref.fasta"}, &bytes.Buffer{}, &buffer, os.Stderr)
  29         c.Assert(exited, check.Equals, 0)
  30         ioutil.WriteFile(tmpdir+"/library.gob", buffer.Bytes(), 0644)
  31
  32         exited = (&exporter{}).RunCommand("export", []string{
  33                 "-local=true",
  34                 "-input-dir=" + tmpdir,
  35                 "-output-dir=" + tmpdir,
  36                 "-output-format=hgvs-onehot",
  37                 "-output-labels=" + tmpdir + "/labels.csv",
  38                 "-ref=testdata/ref.fasta",
  39         }, &buffer, os.Stderr, os.Stderr)
  40         c.Check(exited, check.Equals, 0)
  41         output, err := ioutil.ReadFile(tmpdir + "/out.chr1.tsv")
  42         if !c.Check(err, check.IsNil) {
  43                 out, _ := exec.Command("find", tmpdir, "-ls").CombinedOutput()
  44                 c.Logf("%s", out)
  45         }
  46         c.Check(sortLines(string(output)), check.Equals, sortLines(`chr1.1_3delinsGGC   1       0
  47 chr1.41_42delinsAA      1       0
  48 chr1.161A>T     1       0
  49 chr1.178A>T     1       0
  50 chr1.222_224del 1       0
  51 chr1.302_305delinsAAAA  1       0
  52 `))
  53         output, err = ioutil.ReadFile(tmpdir + "/out.chr2.tsv")
  54         c.Check(err, check.IsNil)
  55         c.Check(sortLines(string(output)), check.Equals, sortLines(`chr2.1_3delinsAAA   0       1
  56 chr2.125_127delinsAAA   0       1
  57 chr2.241_254del 1       0
  58 chr2.258_269delinsAA    1       0
  59 chr2.315C>A     1       0
  60 chr2.470_472del 1       0
  61 chr2.471_472delinsAA    1       0
  62 `))
  63         labels, err := ioutil.ReadFile(tmpdir + "/labels.csv")
  64         c.Check(err, check.IsNil)
  65         c.Check(string(labels), check.Equals, `0,"input1","out.tsv"
  66 1,"input2","out.tsv"
  67 `)
  68
  69         exited = (&exporter{}).RunCommand("export", []string{
  70                 "-local=true",
  71                 "-input-dir=" + tmpdir,
  72                 "-output-dir=" + tmpdir,
  73                 "-output-format=pvcf",
  74                 "-ref=testdata/ref.fasta",
  75         }, &buffer, os.Stderr, os.Stderr)
  76         c.Check(exited, check.Equals, 0)
  77         output, err = ioutil.ReadFile(tmpdir + "/out.chr1.vcf")
  78         c.Check(err, check.IsNil)
  79         c.Log(string(output))
  80         c.Check(sortLines(string(output)), check.Equals, sortLines(`##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
  81 #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  testdata/pipeline1/input1.1.fasta       testdata/pipeline1/input2.1.fasta
  82 chr1    1       .       NNN     GGC     .       .       .       GT      1/1     0/0
  83 chr1    41      .       TT      AA      .       .       .       GT      1/0     0/0
  84 chr1    161     .       A       T       .       .       .       GT      0/1     0/0
  85 chr1    178     .       A       T       .       .       .       GT      0/1     0/0
  86 chr1    221     .       TCCA    T       .       .       .       GT      1/1     0/0
  87 chr1    302     .       TTTT    AAAA    .       .       .       GT      0/1     0/0
  88 `))
  89         output, err = ioutil.ReadFile(tmpdir + "/out.chr2.vcf")
  90         c.Check(err, check.IsNil)
  91         c.Log(string(output))
  92         c.Check(sortLines(string(output)), check.Equals, sortLines(`##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
  93 #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  testdata/pipeline1/input1.1.fasta       testdata/pipeline1/input2.1.fasta
  94 chr2    1       .       TTT     AAA     .       .       .       GT      0/0     0/1
  95 chr2    125     .       CTT     AAA     .       .       .       GT      0/0     1/1
  96 chr2    240     .       ATTTTTCTTGCTCTC A       .       .       .       GT      1/0     0/0
  97 chr2    258     .       CCTTGTATTTTT    AA      .       .       .       GT      1/0     0/0
  98 chr2    315     .       C       A       .       .       .       GT      1/0     0/0
  99 chr2    469     .       GTGG    G       .       .       .       GT      1/0     0/0
 100 chr2    471     .       GG      AA      .       .       .       GT      0/1     0/0
 101 `))
 102
 103         exited = (&exporter{}).RunCommand("export", []string{
 104                 "-local=true",
 105                 "-input-dir=" + tmpdir,
 106                 "-output-dir=" + tmpdir,
 107                 "-output-format=vcf",
 108                 "-ref=testdata/ref.fasta",
 109         }, &buffer, os.Stderr, os.Stderr)
 110         c.Check(exited, check.Equals, 0)
 111         output, err = ioutil.ReadFile(tmpdir + "/out.chr1.vcf")
 112         c.Check(err, check.IsNil)
 113         c.Log(string(output))
 114         c.Check(sortLines(string(output)), check.Equals, sortLines(`#CHROM      POS     ID      REF     ALT     QUAL    FILTER  INFO
 115 chr1    1       .       NNN     GGC     .       .       AC=2
 116 chr1    41      .       TT      AA      .       .       AC=1
 117 chr1    161     .       A       T       .       .       AC=1
 118 chr1    178     .       A       T       .       .       AC=1
 119 chr1    221     .       TCCA    T       .       .       AC=2
 120 chr1    302     .       TTTT    AAAA    .       .       AC=1
 121 `))
 122         output, err = ioutil.ReadFile(tmpdir + "/out.chr2.vcf")
 123         c.Check(err, check.IsNil)
 124         c.Log(string(output))
 125         c.Check(sortLines(string(output)), check.Equals, sortLines(`#CHROM      POS     ID      REF     ALT     QUAL    FILTER  INFO
 126 chr2    1       .       TTT     AAA     .       .       AC=1
 127 chr2    125     .       CTT     AAA     .       .       AC=2
 128 chr2    240     .       ATTTTTCTTGCTCTC A       .       .       AC=1
 129 chr2    258     .       CCTTGTATTTTT    AA      .       .       AC=1
 130 chr2    315     .       C       A       .       .       AC=1
 131 chr2    469     .       GTGG    G       .       .       AC=1
 132 chr2    471     .       GG      AA      .       .       AC=1
 133 `))
 134
 135         outdir := c.MkDir()
 136         exited = (&exporter{}).RunCommand("export", []string{
 137                 "-local=true",
 138                 "-input-dir=" + tmpdir,
 139                 "-output-dir=" + outdir,
 140                 "-output-format=hgvs-numpy",
 141                 "-ref=testdata/ref.fasta",
 142         }, &buffer, os.Stderr, os.Stderr)
 143         c.Check(exited, check.Equals, 0)
 144
 145         f, err := os.Open(outdir + "/matrix.chr1.npy")
 146         c.Assert(err, check.IsNil)
 147         defer f.Close()
 148         npy, err := gonpy.NewReader(f)
 149         c.Assert(err, check.IsNil)
 150         variants, err := npy.GetInt8()
 151         c.Assert(err, check.IsNil)
 152         c.Check(variants, check.HasLen, 6*2*2) // 6 variants * 2 alleles * 2 genomes
 153         c.Check(variants, check.DeepEquals, []int8{
 154                 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, // input1.1.fasta
 155                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // input2.1.fasta
 156         })
 157
 158         f, err = os.Open(outdir + "/matrix.chr2.npy")
 159         c.Assert(err, check.IsNil)
 160         defer f.Close()
 161         npy, err = gonpy.NewReader(f)
 162         c.Assert(err, check.IsNil)
 163         variants, err = npy.GetInt8()
 164         c.Assert(err, check.IsNil)
 165         c.Check(variants, check.HasLen, 7*2*2) // 6 variants * 2 alleles * 2 genomes
 166         c.Check(variants, check.DeepEquals, []int8{
 167                 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, // input1.1.fasta
 168                 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // input2.1.fasta
 169         })
 170
 171         annotations, err := ioutil.ReadFile(outdir + "/annotations.chr1.csv")
 172         c.Check(err, check.IsNil)
 173         c.Logf("%s", string(annotations))
 174         c.Check(string(annotations), check.Equals, `0,"chr1.1_3delinsGGC"
 175 1,"chr1.41_42delinsAA"
 176 2,"chr1.161A>T"
 177 3,"chr1.178A>T"
 178 4,"chr1.222_224del"
 179 5,"chr1.302_305delinsAAAA"
 180 `)
 181         annotations, err = ioutil.ReadFile(outdir + "/annotations.chr2.csv")
 182         c.Check(err, check.IsNil)
 183         c.Check(string(annotations), check.Equals, `0,"chr2.1_3delinsAAA"
 184 1,"chr2.125_127delinsAAA"
 185 2,"chr2.241_254del"
 186 3,"chr2.258_269delinsAA"
 187 4,"chr2.315C>A"
 188 5,"chr2.470_472del"
 189 6,"chr2.471_472delinsAA"
 190 `)
 191 }