from dgrec.example_data import get_example_data_diranalysis
Functions for analyzing mutation rates in DGRec experiments. Computes the mutation rate per base within a specified sequence range from genotype data. The hallmark DGRec signature is elevated A-specific mutation rates within the VR.
mut_rate
def mut_rate(
gen_list:list, # a genotype list with the number of molecules detected
ran:tuple, # the position range in which to compute the mutation rate. If None the rate is computed for the full sequence.
ref_seq:str, # reference sequence
):Computes the mutation rate per base within the specified range. The rate can be computed for specific bases using the base_restriction argument.
data_path=get_example_data_dir()
gen_list=parse_genotypes(os.path.join(data_path,"sacB_genotypes.csv"))
read_ref_file="sacB_ref.fasta"
ref=next(SeqIO.parse(os.path.join(data_path,read_ref_file),"fasta"))
ref_seq=str(ref.seq)
#showing a few example lines
for g,n in gen_list[1:200:20]:
print(n,"\t",g)279 A91G
28 A68C
15 A72G,A79T,A91T
10 A61G,A72G
6 A61G,A68G
6 A68G,A76G,A91G
5 A61T,A79G
4 A86T
4 A72G,A76G,A86G,A91T
3 A61T,A76G,A91GTR_range=(50,119)
before_TR_range=(5,50)
mut_rate_TR=mut_rate(gen_list,TR_range,ref_seq)
for b in mut_rate_TR:
print(f"Mutation rate in VR at {b} positions: {mut_rate_TR[b]:.1e}")
mut_rate_outside_TR=mut_rate(gen_list,before_TR_range,ref_seq)
for b in mut_rate_outside_TR:
print(f"Mutation rate outside VR at {b} positions: {mut_rate_outside_TR[b]:.1e}")Mutation rate in VR at A positions: 1.9e-02
Mutation rate in VR at T positions: 9.9e-04
Mutation rate in VR at G positions: 1.1e-04
Mutation rate in VR at C positions: 2.2e-04
Mutation rate in VR at all positions: 2.9e-03
Mutation rate outside VR at A positions: 2.1e-05
Mutation rate outside VR at T positions: 4.8e-05
Mutation rate outside VR at G positions: 1.8e-04
Mutation rate outside VR at C positions: 5.9e-05
Mutation rate outside VR at all positions: 6.7e-05TR-aware analysis
A DGR copies its template repeat into the variable repeat, making errors at TR adenines and, an order of magnitude lower but still well above background, at TR thymines. Classifying a mutation against that expectation needs two facts kept apart — what kind of position it is, and what the read carries there — because merging them mixes was this molecule converted with did the reverse transcriptase err here.
These helpers key on the TR base, not the reference base. Where the VR is designed to match the TR the two coincide, which is why is_dgrec tests the reference base. At a natural DGR locus the VR is already a retrohoming product with its adenines spent, so the positions that report a new conversion have a non-A reference base and that test is blind to them.
POS_CLASSES, OBSERVATIONS(('A-div', 'A-shared', 'N-div', 'invariant', 'outside'),
('to_TR', 'to_third', 'indel'))TRAlignment
def TRAlignment(
ref_seq:str, # reference amplicon sequence
tr_seq:str, # template repeat, in the same orientation as ref_seq
vr_start:int=None, # 0-based start of the VR; located automatically if None
min_identity:float=0.5, # identity required when locating the VR
):Positional annotation of an amplicon against its template repeat.
Two facts are kept apart, because merging them mixes was this molecule converted with did the reverse transcriptase err here:
position class, a property of the reference A-div TR=A, TR!=VR conversion and RT fidelity, confounded A-shared TR=VR=A RT error only, and only if converted N-div TR!=A, TR!=VR conversion alone, no mutagenesis expected invariant TR=VR!=A nothing expected; mutations are noise outside not covered by the TR
observation, what the read carries to_TR, to_third, or indel (never a conversion signal)
Note this keys on the TR base, not the reference base. Where the VR is designed to match the TR the two coincide, but at a natural DGR locus the VR is already a retrohoming product with its adenines spent, so the positions reporting a new conversion have a non-A reference base.
# a toy amplicon: an 8 bp flank, a 10 bp VR, another 8 bp flank
demo_ref = "ACGTACGT" + "GCTACGATCG" + "ACGTACGT"
demo_tr = "ACTACGATCA"
aln = TRAlignment(demo_ref, demo_tr) # the VR is located automatically
alnTRAlignment(ref=26bp, VR=[8,18), A-div=2, A-shared=2, N-div=0, invariant=6)aln.positions{'A-div': (8, 17),
'A-shared': (11, 14),
'N-div': (),
'invariant': (9, 10, 12, 13, 15, 16)}# The test uses TR adenines against TR G/C. TR-T is excluded from both:
# the reverse transcriptase errs there well above the G/C background.
aln.signal_positions, aln.null_positions, aln.thymine_positions((8, 17, 11, 14), (9, 12, 13, 16), (10, 15))classify_genotype
def classify_genotype(
genotype:str, # genotype string, e.g. "G90A,A101G"
aln:TRAlignment, # positional annotation of the amplicon
):Counts per (position class, observation), plus molecule-level summaries.
dgr_consistent means every mutation is one a conversion carrying reverse-transcriptase errors could produce. fidelity is the fraction of mutations at divergent positions that carry the TR base; it is measured within observed events, so it needs no conversion-tract model.
classify_genotype('G8A,A11G,G17A,A14C', aln){'n_mutations': 4,
'counts': {('A-div', 'to_TR'): 2, ('A-shared', 'to_third'): 2},
'conversion_evidence': 2,
'div_positions_mutated': 2,
'rt_errors': 2,
'unexplained': 0,
'dgr_consistent': True,
'fidelity': 1.0,
'k_TR_A': 4,
'k_TR_GC': 0,
'k_TR_T': 0}annotate_genotype
def annotate_genotype(
genotype:str, # genotype string
aln:TRAlignment, # positional annotation of the amplicon
sep:str='|', # separator between a mutation and its tag
):Render a genotype with each mutation tagged by class, for inspection.
A derived view: the canonical genotype string is unchanged, stays comparable between analyses, and remains the key for UMI grouping. Annotated strings must not be used as keys.
annotate_genotype('G8A,A11G,G17A,A14C', aln)'G8A|Adiv>TR,A11G|Ash>3rd,G17A|Adiv>TR,A14C|Ash>3rd'calibrate_null
def calibrate_null(
gen_list:list, # (genotype, count) pairs from get_genotypes
aln:TRAlignment, # positional annotation of the amplicon
covered:tuple=None, # (start, end) of ref_seq the reads actually span
):Null probability that a noise mutation lands on a TR adenine.
Estimated from sequence outside the VR, where neither conversion nor mutagenesis acts, split by base: adenines and G/C do not carry the same sequencing error rate, and assuming they do misstates the null badly.
Pools per-position rates after trimming the noisiest, rather than pooling all. Off-target reads - failed clusters, index hopping, adapter readthrough - disagree with the reference almost everywhere, so the few surviving read-level filtering inflate whichever positions they reach. Pooling lets one such position set the rate for its whole base class: in the data this was built on that raised the adenine null threefold and reported zero DGR molecules in a sample that is ~7% converted.
Returns (rate at A, rate at G/C, p0).
# a toy library: one candidate molecule, a quiet background at flanking adenines,
# a noisier one at flanking G/C (as real data shows), and the unmutated majority
_flankA = [q for q in aln.outside if aln.ref_seq[q] == "A"]
_flankGC = [q for q in aln.outside if aln.ref_seq[q] in ("G", "C")]
_gen_list = ([("G8A,A11G,G17A,A14C", 3)]
+ [(f"A{q}G", 1) for q in _flankA]
+ [(f"{aln.ref_seq[q]}{q}A", 50) for q in _flankGC]
+ [("", 100000)])
calibrate_null(_gen_list, aln) # (rate at A, rate at G/C, p0)(1.1619375807131642e-05, 0.000498684639232609, 0.022769516728624536)dgr_pvalue
def dgr_pvalue(
genotype:str, # genotype string
aln:TRAlignment, # positional annotation of the amplicon
p0:float, # null probability from calibrate_null
):One-sided binomial tail: are this molecule’s mutations concentrated on TR adenines?
Conditions on the molecule’s own mutation count, so a poor-quality molecule gets more trials and its adenine count is correctly less surprising. TR-T positions are excluded from both classes.
dgr_pvalue('G8A,A11G,G17A,A14C', aln, p0=0.15)0.00050625find_dgr_molecules
def find_dgr_molecules(
gen_list:list, # (genotype, count) pairs from get_genotypes
aln:TRAlignment, # positional annotation of the amplicon
fdr:float=0.05, # Benjamini-Hochberg false discovery rate
p0:float=None, # null probability; calibrated from gen_list if None
covered:tuple=None, # (start, end) of ref_seq the reads span
):Search for molecules whose mutations are concentrated on TR adenines.
Every molecule is a test and the correction runs over all of them, so a handful of strongly mutated molecules can be called against a background of millions. Returns (results, p0), results being (genotype, count, p, info) sorted by p, for the molecules passing the FDR threshold.
hits, p0 = find_dgr_molecules(_gen_list, aln)
[(g, n, f'{p:.2e}') for g, n, p, info in hits][('G8A,A11G,G17A,A14C', 3, '2.69e-07')]