Skip to content

Commit ae98b71

Browse files
committed
VcfCalculatePRS: improved debug output
1 parent 2300d34 commit ae98b71

File tree

1 file changed

+14
-6
lines changed

1 file changed

+14
-6
lines changed

src/VcfCalculatePRS/main.cpp

+14-6
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#include <QFile>
77
#include <QTextStream>
88

9+
//TODO: Handle missing variants (check in BAM for depth, check target region if region is contained > use AF if missing?)
910

1011
class ConcreteTool
1112
: public ToolBase
@@ -103,6 +104,7 @@ class ConcreteTool
103104
}
104105

105106
//iterate over all variants in PRS
107+
int c_found = 0;
106108
for(int i = 0; i < prs_variant_list.count(); ++i)
107109
{
108110
const VcfLine& prs_variant = prs_variant_list[i];
@@ -113,7 +115,6 @@ class ConcreteTool
113115
THROW(FileParseException, "Multi-allelic variants in PRS VCF files are not supported: " + prs_variant.variantToString());
114116
}
115117

116-
int allele_count = 0;
117118
//get all matching variants at this position
118119
QByteArrayList matching_lines = sample_vcf.getMatchingLines(prs_variant.chr(), prs_variant.start(), prs_variant.end(), true);
119120
QByteArrayList matching_variants;
@@ -142,6 +143,7 @@ class ConcreteTool
142143
if(genotype_idx < 0) THROW(FileParseException, "Genotype information is missing for sample variant: " + matching_variants[0]);
143144
QByteArray genotype = format_value_items[genotype_idx].trimmed().replace("|", "/").replace(".", "0");
144145

146+
int allele_count = 0;
145147
if(genotype == "0/1") allele_count = 1;
146148
else if(genotype == "1/1") allele_count = 2;
147149
else THROW(FileParseException, "Invalid genotype '" + genotype + "' in sample variant: " + matching_variants[0]);
@@ -150,16 +152,22 @@ class ConcreteTool
150152
double weight = Helper::toDouble(prs_variant.info("WEIGHT"), "PRS weight");
151153
prs += weight * allele_count;
152154

155+
++c_found;
153156
}
154-
155157
}
158+
156159
// compute percentile
157160
int percentile = -1;
158161
if (percentiles.size() == 100)
159162
{
160-
int i = 0;
161-
for (i = 0; i < percentiles.size(); ++i) if (prs < percentiles[i]) break;
162-
percentile = i + 1;
163+
for (int i = 0; i < percentiles.size(); ++i)
164+
{
165+
if (prs < percentiles[i])
166+
{
167+
percentile = i + 1;
168+
break;
169+
}
170+
}
163171
}
164172

165173
QByteArray percentile_string = (percentile == -1) ? "." : QByteArray::number(percentile);
@@ -171,7 +179,7 @@ class ConcreteTool
171179

172180

173181
//print final PRS
174-
out << column_entries["pgs_id"] << ":\t" << prs << endl;
182+
out << column_entries["pgs_id"] << ": variants_found=" << c_found << " prs=" << prs << " percentile=" << percentile_string << endl;
175183
}
176184

177185
output_tsv->flush();

0 commit comments

Comments
 (0)