Capstone - Clinical Variant Report

The finished program reads a VCF as raw text and writes a ranked report: it separates the metadata from the records, pulls the depth and the allele frequency out of the INFO field, discards calls that fail a quality and depth threshold, splits multi-allelic records into one row per alternate allele, left-normalises every surviving allele against the reference sequence, names the gene each one falls in, orders the result so the clinically interesting rows come first, and closes with a single line accounting for the whole run.

Nothing here is a call into a variant library. You will do the work those libraries do, which is the only way to know what they were deciding on your behalf.

The data

Three inputs, and the program has to be told about all three because a VCF alone cannot answer the questions the report asks.

The VCF. The sandbox mounts variants.vcf, which is the shape you are parsing:

##fileformat=VCFv4.2
##reference=GRCh38
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total read depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele frequency">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
chr17	7676154	.	G	A	99.0	PASS	DP=120;AF=0.51
chr17	43093220	rs80357906	C	T	45.2	PASS	DP=88;AF=0.49
chr7	55191822	.	T	G	88.7	PASS	DP=143;AF=0.32
chr12	25245350	.	C	T	12.4	LowQual	DP=19;AF=0.21

Lines beginning ## are metadata and carry no calls. The single line beginning #CHROM is the column header. Everything after it is a record of at least eight tab-separated fields, of which QUAL may be . and INFO is a semicolon-separated list of key=value pairs mixed with bare flags that have no = at all. ALT may hold several alleles separated by commas, in which case AF holds one value per alternate allele in corresponding order (or fewer values if the upstream caller omitted annotations for some alleles).

The reference sequence. Normalisation is not a string operation; it needs the bases either side of the variant. The program is given the sequence of each chromosome it must handle, keyed by name, and positions in the VCF are 1-based offsets into it.

The gene intervals. A table of gene name, chromosome, start and end, 1-based and inclusive. Intervals overlap: a gene sits inside another gene's span often enough that picking one is a decision the program has to make rather than a case it can ignore.

Requirements

1. Read the metadata lines and the header without treating either as a record. Count every remaining line as one record, whether or not it survives.
2. Treat a record as malformed, count it, and drop it when it has fewer than eight tab-separated fields, when POS is not a run of digits, when its chromosome has no reference sequence, when REF or any ALT allele is empty or holds a character other than A, C, G or T after uppercasing, when an alternate allele is identical to the reference allele, when DP, AF or QUAL is present but not numeric, or when REF does not match the reference sequence at POS. A slice that runs past the end of the chromosome is not equal to REF, so a position beyond the sequence is caught by the same test.
3. Parse INFO into a mapping, ignoring entries without an =. A missing DP is 0. A QUAL of . is 0.0.
4. Split a multi-allelic record into one allele per alternate, taking the matching element of the AF list by position and 0.0 when the list is shorter than the allele list. Every count after this point is per allele, not per record.
5. Filter an allele, count it, and drop it when FILTER is neither PASS nor ., when QUAL is below 30, or when DP is below 20.
6. Left-normalise each surviving allele. While both alleles end in the same base: truncate that base from both if both are longer than one base, otherwise step the position back by one and prepend the reference base at the new position to both. An allele that would need to step back past position 1 has no anchoring base, so count it with the filtered alleles and drop it. Then, while both alleles are longer than one base and begin with the same base, drop that base from both and advance the position by one.
7. Annotate the normalised allele with the gene whose interval overlaps [pos, pos + len(ref) - 1]. When several overlap, take the shortest interval, and break a tie on interval length by the first gene name alphabetically. When none overlaps, the gene is intergenic.
8. Classify each allele as snv when both alleles are one base, mnv when they are equal in length and longer than one base, insertion when the alternate allele is longer, and deletion otherwise. Classify after normalising, not before: trimming changes the answer.
9. Rank annotated alleles above intergenic ones, then by quality descending, then chromosome, position, reference allele and alternate allele ascending. Emit the four counters, the ranked rows, and a summary line holding the number of alleles kept, the number of distinct genes excluding intergenic, and the mean depth of the kept alleles.

Milestones

Run these in order. They build the two pieces the rest of the program leans on and leave the pipeline itself to you.

Parse the file into records with a usable INFO mapping.

Trim a few alleles, and confirm that trimming alone is not normalisation.

The last case trims to (4, 'GG', 'G') and stops there, because both alleles are down to their last shared base and there is nothing left to trim. Against the reference AAAGGGTTTCCC the correct answer is (3, 'AG', 'A'), one base further left. Getting from the first to the second is what requires the reference sequence, and it is requirement 6.

Do the filter accounting over the records you parsed, keeping the reasons apart.

One record here fails three tests at once. The report counts it once, and the order of the tests decides which reason it is recorded under, establishing the deterministic evaluation order specified in the requirements.

What to watch for

• Counting malformed records and filtered alleles on the same footing. A record is malformed once; its alleles are filtered individually. A multi-allelic record can contribute several filtered alleles and one record.
• Classifying the variant before normalising it. CACA>CA looks like a four-base event and is a one-base deletion.
• Trimming without shifting. A deletion in a homopolymer has as many representations as the run is long, and only the leftmost is comparable with anyone else's call.
• Stepping past position 1. The loop that prepends reference bases will run off the front of the chromosome if you do not stop it.
• Trusting the AF list to match the ALT list. It often does not, and indexing past its end is the failure that reaches the report as a wrong allele frequency rather than as an error.
• Letting a bare INFO flag through split("=", 1). An entry with no = unpacks to one element and raises.
• Assuming one gene per position. Overlapping intervals are the normal case, so the tie-break has to be written down and applied.

Exercise

Write the whole report generator in one function. hgvs is not installed here, so the left-normalisation step is one you implement against a plain reference string; every other stage of the brief -- header handling, INFO parsing, the quality and depth filter, multi-allelic splitting, interval annotation and ranking -- is present.

Write report(text), returning a list of output lines. The input arrives in three sections separated by lines containing exactly --: the reference sequences (a chromosome name and its sequence, tab separated, 1-based), the VCF text, and a gene interval table (gene, chrom, start, end, tab separated, 1-based and inclusive). Ignore blank lines. In the VCF section, discard every line beginning with #; count each remaining line as one record.

A record is malformed, and is counted and dropped, when it has fewer than eight tab-separated fields, when POS is not a run of digits, when its chromosome is absent from the reference section, when REF or any ALT allele is empty or contains a character other than A, C, G or T after uppercasing, when an alternate allele equals the reference allele, or when REF does not match the reference sequence at POS. Read DP and AF out of the semicolon-separated INFO field, ignoring any entry without an =; a missing DP is 0, a QUAL of . is 0.0, and a non-numeric DP, AF or QUAL makes the record malformed.

Split a record carrying several comma-separated alternate alleles into one allele each, taking the matching element of the comma-separated AF list by position and 0.0 when the list is shorter. Every later count is per allele. An allele is filtered, and is counted and dropped, when FILTER is neither PASS nor ., when QUAL is below 30, when DP is below 20, or when left-normalisation runs off the start of the chromosome.

Normalise each surviving allele: while both alleles end in the same base, truncate that base from both if both are longer than one base, otherwise step the position back by one and prepend the reference base at the new position to both -- and if the position is already 1 there is no base to anchor to, so the allele is filtered. Then, while both alleles are longer than one base and begin with the same base, drop that base from both and advance the position by one. Annotate the normalised allele with the gene whose interval overlaps [pos, pos + len(ref) - 1] on that chromosome, choosing the shortest interval and then the first gene name alphabetically when several overlap, and intergenic when none does. Classify it as snv when both alleles are one base, mnv when they are equal in length and longer than that, insertion when the alternate allele is longer, and deletion otherwise.

Print records=<n>, then malformed=<n>, then filtered=<n>, then kept=<n>. Follow those four lines with one line per kept allele, ordered by annotated alleles before intergenic ones, then by quality descending, then chromosome ascending, then position ascending, then reference allele, then alternate allele; each line holds <chrom>:<pos>, the reference allele, the alternate allele, the gene, the depth, the allele frequency to exactly four decimal places, and the class, separated by tabs. If nothing is kept, print no variants reported in place of the table. Close with summary: <k> variants in <g> genes, mean depth <d>, where k is the number kept, g the number of distinct genes excluding intergenic, and d the mean depth of the kept alleles to exactly one decimal place, or 0.0 when nothing is kept.