Geneticist & bioinformatician: MQL recipes

Recipes

GC content, CpG ratio, TATA boxes and divergence of a promoter
Hardy-Weinberg equilibrium at a SNP
Recessive disease: binomial and Bayesian risks in a family
Cumulative GC skew to locate the replication origin
Linkage analysis: LOD score and map distance
Three-point testcross: interference and map functions
Microsatellite diversity and loss by drift
Heritability by parent-offspring regression (QR)
Somatic variant call: beta-binomial noise model
Ts/Tv ratio with 2-bit nucleotide codes

executed on MentDB

GC content, CpG ratio, TATA boxes and divergence of a promoter

For a promoter fragment and the same region in a related species, compute the GC content, the CpG observed/expected ratio (above 0.6 suggests a CpG island), the TATA box motifs with a regex, the codon count, and the divergence: edit distance, p-distance and the Jukes-Cantor corrected distance.

Pipeline

01 · countG, C and CpG counts, TATA motif regexstring regexp_count · string count→02 · ratiosGC % and CpG o/ebigdec div→03 · codonsCodons and remaining basesarith iquo · arith irem→04 · divergenceEdit distance and Jukes-Cantor dstring levenshtein_distance · math calc

Script

admin

#1. A promoter fragment and the same region in a related species;
-> "[s]" "GCGCGGTATAAAAGGCGCCGCTCGAGCGTTACGCCGATCGCGGCTATAAATACGCGTTCGCGGCAGCGCGCTTAACGATCG";
-> "[h]" "GCGCGGTATAAAAGGCGTCGCTCGAGCATTACGCCGATCGCGACTATAAATACGCGTTCACGGCAGCGTGCTTAACGATCG";
-> "[n]" (string length [s]);
#2. GC content and CpG observed / expected ratio;
-> "[gc]" (string regexp_count [s] "[GC]");
-> "[oe]" (bigdec div (math calc (concat (string count [s] "CG") "*[n]")) (math calc (concat (string regexp_count [s] "C") "*" (string regexp_count [s] "G"))) 4);
#3. TATA box motifs and complete codons;
-> "[tata]" (string regexp_count [s] "TATA[AT]A[AT]");
#4. Divergence: edit distance, p-distance and Jukes-Cantor correction;
-> "[dist]" (string levenshtein_distance [s] [h]);
-> "[p]" (math calc "[dist]/[n]");
concat [n] " bp = " (arith iquo [n] 3) " codons + " (arith irem [n] 3) " bp, GC = " (bigdec div (math calc "100*[gc]") [n] 4) " %, CpG o/e = " [oe] ", TATA boxes: " [tata] ", " [dist] " substitutions, p = " (math round [p] 4) ", Jukes-Cantor d = " (math round (math calc "-3/4*ln(1-4/3*[p])") 4)

Result

mentdb

81 bp = 27 codons + 0 bp, GC = 61.73 %, CpG o/e = 2.337, TATA boxes: 2, 5 substitutions, p = 0.0617, Jukes-Cantor d = 0.0644

Go further

Search the reverse complement with string reverse and string replace
Use the Kimura 2-parameter distance when transitions and transversions are counted apart
Scan a whole gene in windows with a for loop and string substring

executed on MentDB

Hardy-Weinberg equilibrium at a SNP

From genotype counts in 1000 people, get the allele frequencies as exact fractions, the expected heterozygosity 2pq, the chi-square test against the Hardy-Weinberg proportions (1 degree of freedom, not 2) and the exact probability of the observed number of heterozygotes given the allele counts (Levene formula with lgamma).

Pipeline

01 · frequenciesp and q as exact fractions, 2pqarith to_fraction · arith frac_sub · arith frac_mul→02 · chi2Chi-square statisticstat chi2_test→03 · p-valueChi2 law with 1 dfstat dist_cdf→04 · exactLevene probability with lgammalgamma(x) · math calc

Script

admin

#1. Genotype counts at a SNP in 1000 people;
-> "[AA]" 298; -> "[AB]" 489; -> "[BB]" 213; -> "[n]" 1000;
#2. Allele frequency of A as an exact fraction, expected heterozygosity 2pq;
-> "[p]" (arith to_fraction (concat (math calc "2*[AA]+[AB]") "/" (math calc "2*[n]")));
-> "[q]" (arith frac_sub 1 [p]);
-> "[He]" (arith frac_mul 2 (arith frac_mul [p] [q]));
#3. Chi-square against the Hardy-Weinberg proportions, 1 degree of freedom;
-> "[pd]" (math calc "(2*[AA]+[AB])/(2*[n])");
-> "[x2]" (json load select "rx_hw" (stat chi2_test (concat "[" [AA] ", " [AB] ", " [BB] "]") (concat "[" (math calc "[pd]^2") ", " (math calc "2*[pd]*(1-[pd])") ", " (math calc "(1-[pd])^2") "]")) "/chi2");
-> "[pv]" (- 1 (stat dist_cdf "chi2" "[1]" [x2]));
#4. Exact probability of 489 heterozygotes given the allele counts, with lgamma;
-> "[nA]" (math calc "2*[AA]+[AB]"); -> "[nB]" (math calc "2*[BB]+[AB]");
-> "[Pex]" (math calc "exp(lgamma([n]+1)-lgamma([AA]+1)-lgamma([AB]+1)-lgamma([BB]+1)+[AB]*ln(2)+lgamma([nA]+1)+lgamma([nB]+1)-lgamma(2*[n]+1))");
concat "p = " [p] ", q = " [q] ", 2pq = " [He] " = " (bigdec div (string regexp_extract [He] "^([0-9]+)/" 1) (string regexp_extract [He] "/([0-9]+)$" 1) 4) " (observed " (/ [AB] [n]) "), chi2 = " (math round [x2] 3) ", p-value = " (math round [pv] 3) ", exact P(nAB = 489) = " (math round [Pex] 4)

Result

mentdb

p = 217/400, q = 183/400, 2pq = 39711/80000 = 0.4964 (observed 0.489), chi2 = 0.221, p-value = 0.638, exact P(nAB = 489) = 0.0451

Go further

Sum the exact probabilities not larger than the observed one to get the exact HWE p-value
Compute the inbreeding coefficient F = 1 - Ho / He
Test several SNPs with a json parse_array loop

executed on MentDB

Recessive disease: binomial and Bayesian risks in a family

Both parents carry a recessive allele (risk 1/4 per child). The script gives exact binomial probabilities for 5 children with fractions, checks the tail P(X >= 2) with the regularized incomplete beta function and the binomial law, then the Bayesian carrier risk of an unaffected sibling (2/3) and the risk for her child with a partner from a population where 1 person in 25 is a carrier.

Pipeline

01 · binomialC(5,2) (1/4)^2 (3/4)^3 exactlyarith ncr · arith frac_mul · arith frac_sub→02 · tailP(X >= 2) = I_p(2, 4)beta_inc(a, b, x) · stat dist_cdf→03 · BayesCarrier risk 1/2 / (3/4)arith frac_div→04 · riskRisk for the next generationarith frac_mul

Script

admin

#1. Both parents carriers of a recessive allele: risk 1/4 per child, 5 children;
-> "[pa]" "1/4"; -> "[q]" (arith frac_sub 1 [pa]); -> "[n]" 5;
#2. Exact binomial probabilities with fractions;
-> "[q3]" (arith frac_mul [q] (arith frac_mul [q] [q]));
-> "[p2]" (arith frac_mul (arith ncr [n] 2) (arith frac_mul (arith frac_mul [pa] [pa]) [q3]));
-> "[ge1]" (arith frac_sub 1 (arith frac_mul [q3] (arith frac_mul [q] [q])));
#3. P of X at least 2 with the regularized incomplete beta function;
-> "[ge2]" (math eval "beta_inc(2, 4, 0.25)");
-> "[chk]" (- 1 (stat dist_cdf "binomial" "[5, 0.25]" 1));
#4. Bayes: an unaffected sibling is a carrier with 2/3, partner carrier frequency 1/25;
-> "[car]" (arith frac_div "1/2" (arith frac_sub 1 [pa]));
-> "[risk]" (arith frac_mul [car] (arith frac_mul "1/25" [pa]));
concat "P(2 affected of 5) = " [p2] ", P(at least one) = " [ge1] ", P(at least 2) = " (math round [ge2] 6) " (binomial cdf " (math round [chk] 6) "), sibling carrier risk " [car] ", risk for her child " [risk] " (1 in " (arith frac_div 1 [risk]) ")"

Result

mentdb

P(2 affected of 5) = 135/512, P(at least one) = 781/1024, P(at least 2) = 0.367188 (binomial cdf 0.367188), sibling carrier risk 2/3, risk for her child 1/150 (1 in 150)

Go further

Update the carrier risk after a negative carrier test of known sensitivity
Use stat dist_pdf "binomial" for every number of affected children
Compute the carrier frequency from the disease incidence with Hardy-Weinberg: q = sqrt(incidence)

executed on MentDB

Cumulative GC skew to locate the replication origin

In bacterial genomes, the leading strand is rich in G: the GC skew (G - C) / (G + C) changes sign at the origin and at the terminus of replication. From the G and C counts of 12 blocks of 500 bp, the script computes the skew of each block and its cumulative sum, then draws it: the maximum marks the terminus and the minimum the origin.

Pipeline

01 · skew(G - C) / (G + C) per blockmath calc→02 · cumulateRunning sum in a loopjson load select→03 · plotCumulative skew along the sequenceplot signal

Script

admin

#1. A 6 kb genome window read in 12 blocks of 500 bp: G and C counts per block;
-> "[G]" "[131, 128, 135, 140, 138, 142, 118, 112, 109, 115, 111, 120]";
-> "[C]" "[119, 121, 116, 112, 113, 110, 133, 139, 141, 134, 138, 128]";
#2. GC skew per block (G - C) / (G + C) and its cumulative sum;
-> "[sk]" "["; -> "[cum]" "["; -> "[s]" 0; -> "[sep]" "";
for (-> "[i]" 0) (< [i] 12) (++ "[i]") {
  -> "[g]" (json load select "rx_g" [G] (concat "/[" [i] "]")); -> "[c]" (json load select "rx_c" [C] (concat "/[" [i] "]"));
  -> "[k]" (math calc "([g]-[c])/([g]+[c])"); -> "[s]" (math calc "[s]+[k]");
  -> "[sk]" (concat [sk] [sep] [k]); -> "[cum]" (concat [cum] [sep] [s]); -> "[sep]" ",";
};
-> "[cum]" (concat [cum] "]");
#3. The maximum of the cumulative skew marks the replication terminus side, the sign switch the origin;
plot signal [cum]

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Count G and C of each window directly from the sequence with string substring and string regexp_count
Find the extremum with dsp peaks on the cumulative skew (and on its opposite for the minimum)
Smooth the skew per block with dsp savgol before the sum

executed on MentDB

Linkage analysis: LOD score and map distance

In 20 informative meioses a marker and a disease gene show 3 recombinants. Maximize the LOD score over the recombination fraction, get the 1-LOD support interval and convert theta into centimorgans with the Haldane and Kosambi map functions.

Pipeline

01 · LODlog10 L(theta)/L(0.5)num maximize→02 · intervalRoots of LOD = Zmax - 1num solve→03 · mapHaldane and Kosambimath calc

Script

admin

#1. Linkage of a marker and a dominant disease gene: 20 informative meioses, 3 recombinants;
-> "[n]" 20; -> "[r]" 3;
-> "[lod]" (concat "log(t^" [r] "*(1-t)^(" [n] "-" [r] ")/0.5^" [n] ")");
#2. Maximum likelihood recombination fraction and maximum LOD score;
json load "rx_tu1_lod" (num maximize [lod] "t" 0.001 0.5);
-> "[th]" (json select "rx_tu1_lod" "/x"); -> "[Z]" (json select "rx_tu1_lod" "/value");
#3. 1-LOD support interval: roots of LOD = Zmax - 1 on each side;
-> "[lo]" (num solve (concat [lod] "-" [Z] "+1") "t" 0.001 [th]);
-> "[hi]" (num solve (concat [lod] "-" [Z] "+1") "t" [th] 0.5);
#4. Map distance in cM, Haldane and Kosambi map functions;
-> "[dH]" (math calc "-50*ln(1-2*[th])");
-> "[dK]" (math calc "25*ln((1+2*[th])/(1-2*[th]))");
concat "theta = " (math round [th] 4) ", Zmax = " (math round [Z] 3) " (linkage if >= 3), 1-LOD interval " (math round [lo] 3) " to " (math round [hi] 3) "; distance Haldane = " (math round [dH] 1) " cM, Kosambi = " (math round [dK] 1) " cM"

Result

mentdb

theta = 0.15, Zmax = 2.349 (linkage if >= 3), 1-LOD interval 0.034 to 0.366; distance Haldane = 17.8 cM, Kosambi = 15.5 cM

Go further

Add a second family: LOD scores add, check if the total reaches 3
Draw the LOD curve with plot function between 0 and 0.5

executed on MentDB

Three-point testcross: interference and map functions

From the progeny classes of a Drosophila three-point testcross, the script computes both recombination fractions, the coefficient of coincidence and the interference, the Haldane and Kosambi map distances, and checks the additivity of the Kosambi distances on the outer interval.

Pipeline

01 · fractionsSCO + DCO over nmath round→02 · interference1 - CoCmath calc→03 · Haldane-50.ln(1-2r)math log1p→04 · Kosambi50.atanh(2r)math atanh→05 · backr = tanh(2d)/2math tanh

Script

admin

#1. Drosophila three-point testcross, gene order v - cv - ct, 1448 progeny;
-> "[par]" 1084; -> "[sco1]" 202; -> "[sco2]" 144; -> "[dco]" 18; -> "[n]" 1448;
#2. Recombination fractions of both intervals, double crossovers count in each;
-> "[r1]" (/ (+ [sco1] [dco]) [n]); -> "[r2]" (/ (+ [sco2] [dco]) [n]);
#3. Coefficient of coincidence and interference;
-> "[coc]" (/ (/ [dco] [n]) (* [r1] [r2]));
#4. Map distances in cM: Haldane -50.ln 1-2r with log1p, Kosambi 50.atanh 2r;
-> "[h1]" (* -50 (math log1p (* -2 [r1])));
-> "[k1]" (* 50 (math atanh (* 2 [r1])));
-> "[k2]" (* 50 (math atanh (* 2 [r2])));
#5. Kosambi is additive: back to r for the outer interval with r = tanh 2d / 2;
-> "[r13]" (/ (math tanh (/ (* 2 (+ [k1] [k2])) 100)) 2);
concat "r v-cv = " (math round [r1] 4) ", r cv-ct = " (math round [r2] 4) ", CoC = " (math round [coc] 3) ", interference = " (math round (- 1 [coc]) 3) ", v-cv: Haldane " (math round [h1] 2) " cM, Kosambi " (math round [k1] 2) " cM; cv-ct Kosambi " (math round [k2] 2) " cM; predicted r v-ct = " (math round [r13] 4) " vs observed " (math round (/ (+ [sco1] [sco2]) [n]) 4)

Result

mentdb

r v-cv = 0.1519, r cv-ct = 0.1119, CoC = 0.731, interference = 0.269, v-cv: Haldane 18.11 cM, Kosambi 15.69 cM; cv-ct Kosambi 11.38 cM; predicted r v-ct = 0.247 vs observed 0.239

Go further

Compare with the LOD score of each interval (math log10).
Order the three genes automatically: the DCO class is the rarest one (stat min).

executed on MentDB

Microsatellite diversity and loss by drift

Allele sizes of a microsatellite typed in 15 hatchery fish give the allele frequencies and the unbiased expected heterozygosity. The effective size of the stock over 5 generations is the harmonic mean of the census sizes, which predicts the heterozygosity left after 10 generations.

Pipeline

01 · tableAllele frequenciesstat frequency→02 · loopHe = 1 - sum p^2json count→03 · NeHarmonic mean of the censusstat harm_mean→04 · driftHe.(1 - 1/2Ne)^tmath power→05 · modeMost common allelestat mode

Script

admin

#1. Microsatellite D5S818 typed in 15 diploid fish: 30 allele sizes in bp;
-> "[al]" "[112, 116, 116, 120, 112, 116, 124, 116, 120, 112, 128, 116, 120, 116, 124, 112, 116, 120, 116, 112, 124, 116, 120, 112, 116, 128, 116, 120, 124, 116]";
#2. Allele frequencies and expected heterozygosity He = 1 - sum p2;
json load "rx_p4d_ms" (stat frequency [al]);
-> "[k]" (json count "rx_p4d_ms" "/");
-> "[sp2]" 0;
for (-> "[i]" 0) (< [i] [k]) (++ "[i]") {
	-> "[f]" (json select "rx_p4d_ms" (concat "/[" [i] "]/frequency"));
	-> "[sp2]" (+ [sp2] (* [f] [f]));
};
-> "[he]" (* (/ 30 29) (- 1 [sp2]));
#3. Hatchery stock: census of breeders over 5 generations, Ne = harmonic mean;
-> "[ne]" (stat harm_mean "[60, 12, 45, 50, 8]");
#4. Expected He after 10 generations of drift: He.1 - 1/2Ne ^ t;
-> "[he10]" (* [he] (math power (- 1 (/ 1 (* 2 [ne]))) 10));
concat [k] " alleles, most common " (stat mode [al]) " bp, unbiased He = " (math round [he] 4) ", Ne = " (math round [ne] 1) " vs arithmetic mean " (stat mean "[60, 12, 45, 50, 8]") ", He after 10 generations = " (math round [he10] 4)

Result

mentdb

5 alleles, most common [116] bp, unbiased He = 0.7632, Ne = 18.7 vs arithmetic mean 35, He after 10 generations = 0.5821

Go further

Test the Hardy-Weinberg equilibrium of the genotypes with stat chi2_test.
Add the variance in family size to Ne (Ne = 4N/(Vk + 2)).

executed on MentDB

Heritability by parent-offspring regression (QR)

Heights of 8 sons are regressed on the heights of the father and the mother. The least squares problem is solved with a QR decomposition (R.b = Qt.y) and the narrow-sense heritability is the sum of the two parental slopes.

Pipeline

01 · decomposeA = Q.R, Householdermatrix qr→02 · projectQt.ymatrix transpose→03 · solveR.b = Qt.ymatrix solve→04 · h2b father + b motherjson load select→05 · check|R11| = sqrt(n)matrix get

Script

admin

#1. Height in cm of 8 adult sons with the heights of the father and the mother, design matrix 1, father, mother;
-> "[x]" "[[1, 178, 165], [1, 170, 160], [1, 185, 170], [1, 172, 158], [1, 180, 168], [1, 168, 155], [1, 175, 172], [1, 183, 162]]";
-> "[y]" "[179, 172, 183, 175, 178, 171, 180, 177]";
#2. Least squares by QR decomposition: R.b = Qt.y;
json load "rx_p4d_qr" (matrix qr [x]);
-> "[q]" (json select "rx_p4d_qr" "/Q");
-> "[r]" (json select "rx_p4d_qr" "/R");
-> "[qty]" (matrix mul (matrix transpose [q]) [y]);
-> "[b]" (matrix solve [r] [qty]);
#3. Narrow-sense heritability from the parent-offspring slopes: h2 = b father + b mother;
-> "[bf]" (json load select "rx_p4d_b" [b] "/[1]");
-> "[bm]" (json select "rx_p4d_b" "/[2]");
concat "b = " [b] ", h2 = bf + bm = " (math round (+ [bf] [bm]) 3) ", |R11| = " (math round (math abs (matrix get [r] 1 1)) 4) " = sqrt n = " (math round (math sqrt 8) 4)

Result

mentdb

b = [59.94669063572463,0.2860924987743512,0.4059159993462544], h2 = bf + bm = 0.692, |R11| = 2.8284 = sqrt n = 2.8284

Go further

Same fit in one call with matrix solve on the design matrix (QR inside).
Get the standard error of h2 from (RtR)^-1 with matrix inverse.

executed on MentDB

Somatic variant call: beta-binomial noise model

A candidate somatic variant has 18 alternate reads on a depth of 120. The script computes the tail probability of the panel noise model (beta-binomial) with the Beta function, the Phred quality, and the variant allele fraction with its Beta credible interval.

Pipeline

01 · noiseB(alpha, beta)math beta→02 · tailSum of beta-binomial termsarith ncr→03 · Phred-10.log10(p)math log10→04 · VAFBeta posterior quantilesstat dist_inv

Script

admin

#1. Tumour sample: 18 alternate reads on a depth of 120 at a candidate site;
-> "[n]" 120; -> "[k]" 18;
#2. Noise model of the panel: sequencing errors follow a beta-binomial with alpha 1 and beta 199, mean 0.5 %;
-> "[a]" 1; -> "[b]" 199;
-> "[bab]" (math beta [a] [b]);
#3. Tail probability of the noise P X >= 18, beta-binomial terms C n,j . B j+a, n-j+b / B a,b;
-> "[p]" 0;
for (-> "[j]" [k]) (<= [j] [n]) (++ "[j]") {
	-> "[p]" (+ [p] (/ (* (arith ncr [n] [j]) (math beta (+ [j] [a]) (+ (- [n] [j]) [b]))) [bab]));
};
#4. Variant allele fraction: Beta posterior with a uniform prior and its 95 % credible interval;
-> "[lo]" (stat dist_inv "beta" (concat "[" (+ [k] 1) ", " (+ (- [n] [k]) 1) "]") 0.025);
-> "[hi]" (stat dist_inv "beta" (concat "[" (+ [k] 1) ", " (+ (- [n] [k]) 1) "]") 0.975);
concat "P noise gives >= 18 alt reads = " (math round_sig [p] 3) ", Phred QUAL = " (math round (* -10 (math log10 [p])) 0) ", VAF = " (math round (/ [k] [n]) 3) " [" (math round [lo] 3) "; " (math round [hi] 3) "], " (if (< [p] 0.000001) {"PASS: somatic call"} {"FILTER: noise"})

Result

mentdb

P noise gives >= 18 alt reads = 0.00000000969, Phred QUAL = 80, VAF = 0.15 [0.097; 0.225], PASS: somatic call

Go further

Fit alpha and beta of the noise on normal samples with the method of moments (stat mean, stat variance).
Estimate the tumour purity from the VAF of clonal heterozygous variants (VAF = purity / 2).

executed on MentDB

Ts/Tv ratio with 2-bit nucleotide codes

SNVs of an exome are coded on 2 bits (A=0, C=1, G=2, T=3). Bitwise operations classify them: a transition gives ref XOR alt = 2, bit 0 tells purines from pyrimidines, and a packed 4-bit code finds the C>T/G>A deamination signature. The Ts/Tv ratio checks the quality of the call set.

Pipeline

01 · loopRead each ref>altatom get→02 · transitionref XOR alt = 2math bit_xor→03 · ringPurine bitmath bit_and→04 · packref.4 OR altmath bit_or→05 · ratioTs/Tvmath round

Script

admin

#1. SNVs called in an exome, written ref>alt;
-> "[snv]" "C>T,G>A,A>G,T>C,C>A,G>T,C>T,A>T,G>A,T>C,C>T,G>C,A>G,C>T,G>A";
-> "[n]" (atom size [snv] ",");
-> "[ts]" 0; -> "[cpg]" 0; -> "[pur]" 0;
#2. 2-bit code A=0, C=1, G=2, T=3: purines have bit 0 off, a transition keeps the ring class, so ref XOR alt = 2;
for (-> "[i]" 1) (<= [i] [n]) (++ "[i]") {
	-> "[s]" (atom get [snv] [i] ",");
	-> "[r]" (case (equal (string substring [s] 0 1) "A") {0} (equal (string substring [s] 0 1) "C") {1} (equal (string substring [s] 0 1) "G") {2} {3});
	-> "[a]" (case (equal (string substring [s] 2 3) "A") {0} (equal (string substring [s] 2 3) "C") {1} (equal (string substring [s] 2 3) "G") {2} {3});
	if (== (math bit_xor [r] [a]) 2) {++ "[ts]";};
	if (== (math bit_and [r] 1) 0) {++ "[pur]";};
	#3. Pack ref and alt in one 4-bit code: ref.4 OR alt, C>T = 7 and G>A = 8 are the deamination signature;
	-> "[code]" (math bit_or (* [r] 4) [a]);
	if (or (== [code] 7) (== [code] 8)) {++ "[cpg]";};
};
concat [n] " SNVs, transitions = " [ts] ", transversions = " (- [n] [ts]) ", Ts/Tv = " (math round (/ [ts] (- [n] [ts])) 2) ", purine ref = " [pur] ", C>T/G>A = " [cpg] " = " (math round (* 100 (/ [cpg] [n])) 1) " %"

Result

mentdb

15 SNVs, transitions = 11, transversions = 4, Ts/Tv = 2.75, purine ref = 8, C>T/G>A = 7 = 46.7 %

Go further

A whole-genome Ts/Tv near 2.0-2.1 and an exome near 2.8-3.0 are expected: a low value means false positives.
Count the 96 trinucleotide contexts of the mutational signatures with a 6-bit code.