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