Repository navigation
Read each sample against the cohort: deletions, relatedness, series and minority variants - #8
Merged
Merged
Conversation
…ohort Each sample's all-positions VCF is reduced by a new per-sample step, DEPTH_PROFILE, to its depth per window, the stretches no read covers and a per-gene table, published under stats/. The report compares the stretches across the samples of each reference: one a sample lacks while the others read it is a deletion, private or shared; one nobody reads is a repeat or a part of the reference none of these genomes has, and is listed apart rather than counted as anyone's. Samples under --report_depth_min are not assessed, since stretches without reads turn up there by chance. A Deletions panel lists the regions with the genes they remove and zooms the landscape onto each, and <samplesheet>_deletions.tsv keeps them with every carrier's own coordinates. The landscape gains a Deletions track and a SNPs / kb track that divides each bin's SNPs by the positions deep enough to call, so a bin half of which was not read no longer looks half as variable.
…roup A new step, SNP_DISTANCES, counts the SNPs between every two consensus sequences over the positions both called, so a no-call, a masked position or a mixed site is never a difference: the distances a tree built on the same sequences sees. It writes the square matrix and one row per pair with how many of the reference's variable positions each comparison rests on. Samples of different references are not compared. The report gains a Relatedness page: a heatmap in the order that keeps every single-linkage cluster together, the clusters at a threshold that moves live (--snp_cluster_threshold, 12 by default, the usual TB cut), and the samples further than the threshold from every other member of their samplesheet group, in groups that are otherwise tight. Those are flagged GROUP_MISMATCH (WARN). On the cohort that motivated it, with distances approximated from its SNP matrix, 20 samples are flagged, 18 of them the samples typed as another lineage than their line; a rule that also asked for a close sample of another group flagged none.
For every samplesheet group followed over time (a patient, a passage line) the report lists the SNPs each later sample carries fixed that the group's first time point did not: new where that first time point was read at the site without the allele, risen where it held the allele as a minority, unknown where it was not read there, and lost for those fixed at the start and read without the allele later. The report now reads the SNP matrix for the depth, so a site the start never read is not counted as new; the matrix is built before the report for that reason. A chart follows each series through the median of each time point, coloured by treatment, and leaves out the samples the QC fails or places outside their series, which would otherwise flatten every real series: one mislabelled sample carries thousands of SNPs. On the Resistance page each series' grade 1-2 mutations are laid out per time point with those acquired since the start set apart.
A Minority variants panel counts each sample's calls below fixation, their spread of allele fractions and how many rest on three alternate reads or fewer, where a sequencing error and a real minority look the same. When the samplesheet names each sample's DNA extract (dna_id, extract, biosample...), libraries of the same DNA are compared: a real minority is in the DNA and another library calls it too, while an error is not reproduced. The share reproduced in each band of allele fraction says where the noise ends, against fixed calls as the ceiling. Pairs that share a FASTQ file are left out, since a merged sample shares the reads of the runs it was merged from and its agreement with them proves nothing, and so are libraries mapped against different references. A call counts as missed only where the other library was read, which the SNP matrix's depth says; without the matrix the comparison is not made rather than made biased.
The stub run already checks that each optional cohort output stops when its flag is false; the distances step joins the list.
The coverage card now counts the stretches some samples have no reads for while others do, private and shared, and those nobody reads. The relatedness card says how many samples belong to a cluster in plainer words.
The depth profile counted a position as callable above --consensus_min_dp (7), but the consensus writes a base only from --allpos_min_cov reads (30) and N in between, so SNPs per callable kb read low. Callable now means --allpos_min_cov reads or more, which the depth profile is given.
Stretches without reads were merged into a region whenever they overlapped at all, and the region was classed as a whole. A sample's 6 kb private deletion that overlapped a gap nine of ten samples share was swallowed into a "nobody reads it" region; ten neighbouring 700 bp deletions, none lacked by more than two samples, chained into one; and a carrier's two separate stretches were reported, and drawn, as one span across the region. Two samples' stretches are now one deletion when each covers at least half of the other, measured against the stretch that opened the region, and every carrier keeps its own. Profiles are keyed on the sample and its reference, so a sample mapped against two references is no longer a mix of both, and they are read into little memory: the windows reduced to the landscape's bins as they are read, gene breadths kept as one array per sample. A thousand samples used to take about 3 GB of the report's fixed 4 GB, which now also grows on retry. Naming a region's genes looks only at the genes each carrier barely reads. The depth tables and the SNP distances are produced whether or not a report is made. The class nearly every sample lacks is labelled as such rather than "no sample reads it", a sample alone on its reference no longer lights up the deletion track, clicking a stretch offers to reset the zoom and opens a nearly-unread one on the Missing track, and closing the full view redraws these panels.
The distances were keyed on the sample alone. A sample mapped against two references was listed under only one of them, with the distances of whichever reference came last, and in the square matrix its second row overwrote the first. Pairs are now keyed on the reference too, the matrix names such a sample sample@reference, and a sample far from its group is looked for within each reference. The page rounded the share of variable positions a pair compared, the flag used the exact share: a pair at 49.5% was drawn and clustered but left out of the flag. Both now round down. A reference without variable positions blanked every pair, although its samples compared everything there was; they now count as fully compared. A truncated consensus stopped the distances of every reference; it is now left out with a warning that names it. The histogram's bins meet at the threshold, so no bar straddles it. The controls are built once, so dragging the slider no longer loses it under the pointer. The heatmap is drawn when it is seen, its order and colours worked out once per reference, and a large cluster's members wrap, the first forty named.
A series was read against its earliest sample whatever the QC made of it: a failed or swapped first time point made every later sample gain what that one lacked. The start is now the earliest time point with a sample the QC keeps and places in the series, and a sample whose time cannot be read has no place in it rather than being taken for the last. Times were read by their first number. A day-first date sorted by its day, so 03/06/2020 came before 15/01/2020; an ISO date by its year, so a year's samples were one time point; and "baseline", with no number, came last. Dates are now read as dates and a word for the start as time 0. A DNA-extract column such as replicate_of, which matches the group pattern, is no longer taken for the group. A site read at exactly --consensus_min_dp reads, which the consensus leaves uncalled, no longer counts as read. The chart placed times that are not numbers by row, one step per sample, so a time point's samples could not share a median; they are now categories, and every tick carries the time as the samplesheet writes it. Groups beyond the twelfth get hues of their own, and a sample name is no longer run as HTML by the tooltip. The summary's gain leaves out the samples the chart leaves out and takes the median at each series' last time point.
A call counted as missed by another library of the same DNA wherever that library had --consensus_min_dp reads. At 8x a 5% minority is 0.4 alternate reads: its absence says nothing, and counting it made the low bands look like noise however real they were. A miss now counts only above --consensus_min_dp reads and where at least five alternate reads were expected at the call's fraction. The floor was the first band from the bottom to reach 80% reproduced, even when a band above it fell short again. It is now the lowest band from which every band up to fixation reaches it, and a failing top band, or too few comparisons to judge, is said as such. A DNA-extract column under which no two independent libraries share a DNA no longer reads as a floor never reached: the page says nothing was compared.
The SNP matrix writes NA where a sample was called at a fraction its files do not keep. The report read it as an empty cell, so with the site's depth beside it, a series' first time point that carried the allele counted as read without it, which made the later call new, and a library that called it counted as having missed it. An NA cell is now a call of unknown fraction: unknown in a series, reproduced among the libraries of one DNA.
The relatedness entry credited GROUP_MISMATCH with finding exactly the samples typed as another lineage than their line, at distances that were only approximated; on those distances 18 of the 20 samples it flags are. The resistance example now says which lines acquire which atpE mutation. The entries also say what the fixes changed: a missed minority needs a library that could have called it, a series starts at a sample it can be read against and reads dates as dates, deletions are joined by reciprocal overlap, and callable means the depth from which the consensus calls a base.
A call counted as reproduced wherever the other library of its DNA made it, but as missed only where that library was deep enough to have made it. Whether a comparison counted depended on its outcome: a 5% call the other library made at 40 reads was reproduced, one it did not make at 40 reads was left out as too thin to tell, and the thin bands filled with the calls made by luck. The other library's depth, from its call where it made one and from the SNP matrix where it did not, now decides for both outcomes.
5 of 6 tasks
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What changed and why
Four analyses that read each sample against the rest of the cohort. Each one is driven by generic samplesheet columns and hidden when its data are absent, so the report works the same for any organism.
DEPTH_PROFILE, reduces each all-positions VCF to its depth per window, the stretches no read covers and a per-gene table. The report compares those stretches across the samples of each reference. One a sample lacks while the others read it is a deletion, private or shared. One that nearly every sample lacks (90% or more) is a repeat, or a part of the reference none of these genomes has, and is listed apart. Two samples' stretches are one deletion when each covers at least half of the other. The landscape gains a Deletions track and a SNPs / kb track, whose callable positions are those the consensus calls a base at (--allpos_min_covreads).SNP_DISTANCES, counts the SNPs between every two consensus sequences over the positions both called. The Relatedness page draws the heatmap and the clusters at a threshold that moves live (--snp_cluster_threshold, 12). It flagsGROUP_MISMATCH(WARN) on a sample far from the rest of a samplesheet group that is otherwise tight. On the cohort that motivated it, with distances approximated from its SNP matrix, 18 of the 20 samples it flags are the ones typed as another lineage than their line.--consensus_min_dpreads, with at least five alternate reads expected at the call's fraction.Found and fixed in review before this PR:
NAcell of the matrix was read as a site without the allele.Stacked on #7: retarget once that is merged.
Outputs and defaults
DEPTH_PROFILE(per sample) andSNP_DISTANCES(--make_snp_distances true). Both run whether or not a report is made.<samplesheet>_deletions.tsv,<samplesheet>_snp_distances.tsvand<samplesheet>_snp_distances_pairs.tsv. Each sample'sstats/folder gainsdepth_windows.tsv,zero_depth.tsvandgene_depth.tsv.--snp_cluster_threshold(12),--make_snp_distances(true),--depth_window(1000) and--deletion_min_len(200).qc_flags.tsvcan carryGROUP_MISMATCH, which turns a PASS sample into WARN.QC_REPORTreads the SNP matrix, the depth tables and the distances, and its memory grows with each retry (4 GB times the attempt).How it was verified
Each commit was also tested on its own. CI passes on the branch head, every job green, the stub run included on Nextflow 24.04.2 and the latest stable release: run 35875982736.
The cohort's existing report was re-rendered from its own payload with this code. All eight pages load without console errors. The deletions, distances and replicate depths in that preview were stand-ins, so its numbers for those panels are not results.
Checklist
tests/run_tests.shpasses (unit, js and lint locally; the pipeline leg in CI)ruff check .is clean, and I did not restyle code the change does not touchmodules/has astub:block (DEPTH_PROFILE,SNP_DISTANCES), and the stub run reaches them in CItests/data/was not edited by hand (fixtures come frommake_test_data.py)CHANGELOG.mdentry and, where relevant, a docs update