1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
/*
* Copyright (C) 2014-2021 Brian L. Browning
*
* This file is part of Beagle
*
* Beagle is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Beagle is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package vcf;
import java.util.Arrays;
import java.util.stream.IntStream;
/**
* <p>Class {@code LowMafGTRc} stores genotypes for a list of samples
* at a marker. Instances of class {@code LowMafGTRec} store lists of
* haplotypes carrying each non-major or missing allele. All genotypes are
* considered to be unphased if any sample has an unphased or missing
* genotype.</p>
*
* <p>Instances of class {@code LowMafGTRec} are immutable.</p>
*
* @author Brian L. Browning {@code <browning@uw.edu>}
*/
public final class LowMafGTRec implements GTRec {
private final Marker marker;
private final Samples samples;
private final int nHaps;
private final int majorAllele;
private final int[][] hapIndices;
private final int[] missingSamples;
private final boolean isPhased;
/**
* Constructs a new {@code LowMafGTRec} representing
* the specified VCF record's GT format field data.
* @param listRep the VCF record genotype data
* @throws NullPointerException if {@code listRep == null}
*/
public LowMafGTRec(VcfRecGTParser.HapListRep listRep) {
this.marker = listRep.marker();
this.samples = listRep.samples();
this.nHaps = samples.size()<<1;
this.majorAllele = listRep.majorAllele();
this.hapIndices = listRep.hapLists(true);
this.missingSamples = listRep.missingSamples();
this.isPhased = listRep.isPhased();
}
@Override
public boolean isPhased(int sample) {
if (sample < 0 || sample >= this.samples().size()) {
throw new IndexOutOfBoundsException(String.valueOf(sample));
}
return isPhased;
}
@Override
public boolean isPhased() {
return isPhased;
}
@Override
public Samples samples() {
return samples;
}
@Override
public int size() {
return nHaps;
}
@Override
public Marker marker() {
return marker;
}
@Override
public int allele1(int sample) {
return get(sample<<1);
}
@Override
public int allele2(int sample) {
return get((sample<<1) | 0b1);
}
@Override
public int get(int hap) {
if (hap < 0 || hap >= nHaps) {
throw new IndexOutOfBoundsException(String.valueOf(hap));
}
for (int j=0; j<hapIndices.length; ++j) {
if (j != majorAllele) {
if (Arrays.binarySearch(hapIndices[j], hap) >= 0) {
return j;
}
}
}
return Arrays.binarySearch(missingSamples, (hap>>1))>=0 ? -1 : majorAllele;
}
@Override
public int[] alleles() {
int[] ia = IntStream.range(0, nHaps)
.map(h -> majorAllele)
.toArray();
for (int al=0; al<hapIndices.length; ++al) {
if (al != majorAllele) {
for (int h : hapIndices[al]) {
ia[h] = al;
}
}
}
for (int s : missingSamples) {
int h1 = s<<1;
int h2 = h1 | 0b1;
ia[h1] = -1;
ia[h2] = -1;
}
return ia;
}
/**
* Returns the major allele. If there are two major alleles with the same
* allele count, the smaller allele is returned.
* @return the major allele
*/
public int majorAllele() {
return majorAllele;
}
/**
* Returns the number of copies of the specified allele.
* @param allele an allele
* @return the number of copies of the specified allele
* @throws IndexOutOfBoundsException if
* {@code allele < 0 || allele >= this.marker().nAlleles()}
*/
public int alleleCount(int allele) {
if (allele==majorAllele) {
int n = nHaps - (missingSamples.length<<1);
for (int al=0; al<hapIndices.length; ++al) {
if (al!=majorAllele) {
n -= hapIndices.length;
}
}
return n;
}
else {
return hapIndices[allele].length;
}
}
/**
* Returns the data represented by {@code this} as a VCF
* record with a GT format field. The returned VCF record
* will have missing QUAL and INFO fields, will have "PASS"
* in the filter field, and will have a GT format field.
* @return the data represented by {@code this} as a VCF
* record with a GT format field
*/
@Override
public String toString() {
return GTRec.toVcfRec(this);
}
}
|
