export_test.go

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