Skip to content

feat: Rel num diff #40

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
WardDeb merged 4 commits into from
Oct 24, 2026
Merged

feat: Rel num diff #40

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9e9a2a2
figure order, include num diff
WardDeb Oct 23, 2026
6d88f1f
Merge branch 'main' of github.com:deeptools/deeptools25 into refine_v…
WardDeb Oct 24, 2026
0bd5488
include relative shift
WardDeb Oct 24, 2026
ea95d64
num diffs finalized
WardDeb Oct 24, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 20 additions & 0 deletions benchmark/benchmark.smk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,26 @@ computeMatrix_samples = {
] * 5
}

VERIFY_bamCompare_samples = {
**bamCompare_samples,
'bcom_human_rna': ['human_rna_SRR28012902', 'human_wgs_SRR15494527'],
'bcom_human_chiprna': ['human_chip_SRR28592124', 'human_rna_SRR28012902'],
'bcom_triticum_rna': ['triticum_rna_SRR27822150', 'triticum_wgs_SRR27887047'],
'bcom_triticum_chiprna': ['triticum_chip_SRR1686799', 'triticum_rna_SRR27822150']
}

VERIFY_multibamSummary_samples = {
**multibamSummary_samples,
'mbs_human_6': multibamSummary_samples['mbs_human_3'] * 2,
'mbs_triticum_6': multibamSummary_samples['mbs_triticum_3'] * 2
}

VERIFY_computeMatrix_samples = {
**computeMatrix_samples,
'cm_human_6': computeMatrix_samples['cm_human_3'] * 2,
'cm_triticum_6': computeMatrix_samples['cm_triticum_3'] * 2
}

ALSIEVE_THREADS = 20
BAMCOVERAGE_THREADS = 20
BAMCOMPARE_THREADS = 20
Expand Down
54 changes: 31 additions & 23 deletions benchmark/results/numerical_diff.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,23 +1,31 @@
tool sample n_reads_dt3 n_reads_dt4 n_identical n_only_dt3 n_only_dt4 n_flag_mismatch n_qlen_mismatch n_cigar_mismatch n_mapq_mismatch chroms_only_dt3 chroms_only_dt4 n_values pearson_r n_nan_mismatch n_nonidentical frac_identical max_abs_diff max_abs_diff_loc max_rel_diff n_diff_gt_1e-06 n_diff_gt_0.001 n_diff_gt_0.1 shape_dt3 shape_dt4 shape_mismatch n_bins_dt3 n_bins_dt4 n_common_bins
alignmentSieve alsieve_human_chip 29727000.0 29727000.0 29727000.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_human_rna 250070850.0 250070850.0 250070850.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_human_wgs 934945514.0 934945514.0 934945514.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_triticum_chip 46427552.0 46427552.0 46427552.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_triticum_rna 22650946.0 22650946.0 22650946.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_triticum_wgs 1051310232.0 1051310232.0 1051310232.0 0.0 0.0 0.0 0.0 0.0 0.0
bamCoverage bcov_human_chip 3099750718.0 0.999999999001162 0.0 1835797727.0 0.4077595606834859 0.05078125 1:125179260 0.0014522272584559 1835797727.0 1352121305.0 0.0
bamCoverage bcov_human_rna 3099750718.0 0.9999999999991716 0.0 301140058.0 0.9028502336490144 4.75 14:49586730 0.0002870736248038 301140058.0 42751849.0 24690.0
bamCoverage bcov_human_wgs 3099750718.0 0.9999999894310986 0.0 2908998687.0 0.0615378617036261 0.05078125 3:93470430 0.0431545081092797 2908998687.0 2310496679.0 0.0
bamCoverage bcov_triticum_chip 14547261565.0 0.9999999922069556 0.0 2889008442.0 0.8014053415420251 0.010009765625 1B:80735790 0.0004098165701451 2889008442.0 1416822765.0 0.0
bamCoverage bcov_triticum_rna 14547261565.0 0.9999999999885468 0.0 154092080.0 0.9894074854355588 0.0625 3A:430336380 0.0010941214464004 154092080.0 97398130.0 0.0
bamCoverage bcov_triticum_wgs 14547261565.0 0.9999999675942488 0.0 13761865696.0 0.0539892587680986 0.0498046875 3B:101914160 0.0387326507523322 13761865696.0 10339414376.0 0.0
bamCompare bcom_human_chip 3099750718.0 0.9999956228802924 0.0 2909691427.0 0.0613143792164773 0.0050001144409179 1:2659940 1.0 2909691407.0 2313910977.0 0.0
bamCompare bcom_triticum_chip 14547261565.0 0.9999911613864244 0.0 13807666841.0 0.0508408211879154 0.0050001144409179 1A:722100 1.0 13807666711.0 12187591376.0 0.0
multiBamSummary mbs_human_3 930288.0 0.9999999999999984 0.0 6.0 0.9999935503843972 3.0 21:8210000-8220000 sample1 1.3432190674782218e-07 6.0 6.0 4.0 (310096, 6) (310096, 6) False 310096.0 310096.0 310096.0
multiBamSummary mbs_human_15 4651440.0 0.9999999999999982 0.0 30.0 0.9999935503843972 3.0 21:8210000-8220000 sample1 1.3432190674782218e-07 30.0 30.0 20.0 (310096, 18) (310096, 18) False 310096.0 310096.0 310096.0
multiBamSummary mbs_triticum_3 4364211.0 0.9999999999999996 0.0 4.0 0.999999083454031 0.0600000000558793 Un:93410000-93420000 sample2 7.382661280557464e-08 4.0 4.0 0.0 (1454737, 6) (1454737, 6) False 1454737.0 1454737.0 1454737.0
multiBamSummary mbs_triticum_15 21821055.0 0.9999999999999996 0.0 20.0 0.999999083454031 0.0600000000558793 Un:93410000-93420000 sample2 7.382661280557464e-08 20.0 20.0 0.0 (1454737, 18) (1454737, 18) False 1454737.0 1454737.0 1454737.0
computeMatrix cm_human_3 120185364.0 0.9999999999999996 0.0 1173160.0 0.9902387448774544 0.0500000000465661 row94221 col459 2.821192800318795e-05 847660.0 984.0 0.0 (200310, 600) (200310, 600) False
computeMatrix cm_human_15 600926820.0 0.9999999999999998 0.0 5865800.0 0.9902387448774544 0.0500000000465661 row94221 col459 2.821192800318795e-05 4238300.0 4920.0 0.0 (200310, 3000) (200310, 3000) False
computeMatrix cm_triticum_3 87957837.0 1.0 0.0 590633.0 0.9932850440603718 0.0007810000024619 row1809 col59 6.230529601420289e-08 547708.0 0.0 0.0 (146597, 600) (146597, 600) False
computeMatrix cm_triticum_15 439789185.0 0.9999999999999998 0.0 2953165.0 0.9932850440603718 0.0007810000024619 row1809 col59 6.230529601420289e-08 2738540.0 0.0 0.0 (146597, 3000) (146597, 3000) False
tool sample canonical_only n_reads_dt3 n_reads_dt4 n_identical n_only_dt3 n_only_dt4 n_flag_mismatch n_qlen_mismatch n_cigar_mismatch n_mapq_mismatch chroms_only_dt3 chroms_only_dt4 n_values pearson_r n_nan_mismatch n_nonidentical frac_identical max_abs_diff max_abs_diff_loc max_rel_diff n_diff_gt_1e-06 n_diff_gt_0.001 n_diff_gt_0.1 shape_dt3 shape_dt4 shape_mismatch n_noncanonical_bins_excluded n_bins_dt3 n_bins_dt4 n_common_bins n_noncanonical_rows_excluded n_region_mismatch
alignmentSieve alsieve_human_chip True 29026022.0 29026022.0 29026022.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_human_rna True 163156764.0 163156764.0 163156764.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_human_wgs True 925339008.0 925339008.0 925339008.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_triticum_chip True 42233778.0 42233778.0 42233778.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_triticum_rna True 22211940.0 22211940.0 22211940.0 0.0 0.0 0.0 0.0 0.0 0.0
alignmentSieve alsieve_triticum_wgs True 1010572518.0 1010572518.0 1010572518.0 0.0 0.0 0.0 0.0 0.0 0.0
bamCoverage bcov_human_chip True 3088286401.0 0.9999999983742904 0.0 1826772360.0 0.4084835009445745 0.05078125 1:125179260 0.0014522272584559 1826772360.0 1344897310.0 0.0
bamCoverage bcov_human_rna True 3088286401.0 0.9999999999990576 0.0 300545879.0 0.9026819925436054 4.75 14:49586730 0.0002870736248038 300545879.0 42661279.0 13110.0
bamCoverage bcov_human_wgs True 3088286401.0 0.9999999889088896 0.0 2898179779.0 0.0615573160372829 0.05078125 3:93470430 0.0431545081092797 2898179779.0 2301816059.0 0.0
bamCoverage bcov_triticum_chip True 14066280851.0 0.999999991245448 0.0 2757018522.0 0.803998046732872 0.010009765625 1B:80735790 0.0004098165701451 2757018522.0 1334648905.0 0.0
bamCoverage bcov_triticum_rna True 14066280851.0 0.9999999999885788 0.0 150925920.0 0.989270374905868 0.0625 3A:430336380 0.0010941214464004 150925920.0 95480100.0 0.0
bamCoverage bcov_triticum_wgs True 14066280851.0 0.9999999634638266 0.0 13387921476.0 0.0482259228424103 0.0498046875 3B:101914160 0.0387326507523322 13387921476.0 10049695846.0 0.0
bamCompare bcom_human_chip True 3088286401.0 0.9999956112004996 0.0 2898793239.0 0.0613586751340942 0.0050001144409179 1:2659940 1.0 2898793219.0 2305249729.0 0.0
bamCompare bcom_triticum_chip True 14066280851.0 0.9999906145946544 0.0 13425959101.0 0.0455217521093699 0.0050001144409179 1A:722100 1.0 13425959041.0 11871436416.0 0.0
bamCompare bcom_human_rna True 3088286401.0 0.9999968923489088 0.0 2898384589.0 0.0614909977062065 0.0050001144409179 1:945650 1.0 2898384589.0 2328556520.0 0.0
bamCompare bcom_human_chiprna True 3088286401.0 0.9999981993701154 0.0 1083576779.0 0.6491333256367889 0.0050010681152343 MT:3180 1.0 1083576779.0 1039113469.0 0.0
bamCompare bcom_triticum_rna True 14066280851.0 0.9999763789353484 0.0 13390033346.0 0.0480757857861144 0.0050001144409179 1A:33210 1.0 13390033336.0 11918196715.0 0.0
bamCompare bcom_triticum_chiprna True 14066280851.0 0.9999993970710996 0.0 2851195582.0 0.7973028114395069 0.0050001144409179 1A:65480 1.0 2851195582.0 412341680.0 0.0
multiBamSummary mbs_human_3 True 926517.0 0.9999999999999968 0.0 5.0 0.9999946034449448 3.0 21:8210000-8220000 sample1 1.3432190674782218e-07 5.0 5.0 3.0 (310096, 6) (310096, 6) False 1257.0 310096.0 310096.0 308839.0
multiBamSummary mbs_human_15 True 4632585.0 0.9999999999999972 0.0 25.0 0.9999946034449448 3.0 21:8210000-8220000 sample1 1.3432190674782218e-07 25.0 25.0 15.0 (310096, 18) (310096, 18) False 1257.0 310096.0 310096.0 308839.0
multiBamSummary mbs_triticum_3 True 4219914.0 1.0 0.0 3.0 0.9999992890850382 0.0600000000558793 1A:554460000-554470000 sample2 7.209311553508035e-08 3.0 3.0 0.0 (1454737, 6) (1454737, 6) False 48099.0 1454737.0 1454737.0 1406638.0
multiBamSummary mbs_triticum_15 True 21099570.0 1.0 0.0 15.0 0.9999992890850382 0.0600000000558793 1A:554460000-554470000 sample2 7.209311553508035e-08 15.0 15.0 0.0 (1454737, 18) (1454737, 18) False 48099.0 1454737.0 1454737.0 1406638.0
multiBamSummary mbs_human_6 True 1853034.0 0.9999999999999968 0.0 10.0 0.9999946034449448 3.0 21:8210000-8220000 sample1 1.3432190674782218e-07 10.0 10.0 6.0 (310096, 9) (310096, 9) False 1257.0 310096.0 310096.0 308839.0
multiBamSummary mbs_triticum_6 True 8439828.0 1.0 0.0 6.0 0.9999992890850382 0.0600000000558793 1A:554460000-554470000 sample2 7.209311553508035e-08 6.0 6.0 0.0 (1454737, 9) (1454737, 9) False 48099.0 1454737.0 1454737.0 1406638.0
computeMatrix cm_human_3 True 120145557.0 1.0 0.0 1171253.0 0.9902513831618426 0.0500000000465661 row94221 col459 2.821192800318795e-05 846173.0 970.0 0.0 (200310, 600) (200310, 600) False 67.0 0.0
computeMatrix cm_human_15 True 600727785.0 1.0 0.0 5856265.0 0.9902513831618426 0.0500000000465661 row94221 col459 2.821192800318795e-05 4230865.0 4850.0 0.0 (200310, 3000) (200310, 3000) False 67.0 0.0
computeMatrix cm_triticum_3 True 83199237.0 1.0000000000000002 0.0 578304.0 0.993049167025414 0.0007810000024619 row1809 col59 6.230529601420289e-08 537345.0 0.0 0.0 (146597, 600) (146597, 600) False 7931.0 0.0
computeMatrix cm_triticum_15 True 415996185.0 1.0 0.0 2891520.0 0.993049167025414 0.0007810000024619 row1809 col59 6.230529601420289e-08 2686725.0 0.0 0.0 (146597, 3000) (146597, 3000) False 7931.0 0.0
computeMatrix cm_human_6 True 240291114.0 1.0000000000000002 0.0 2342506.0 0.9902513831618426 0.0500000000465661 row94221 col459 2.821192800318795e-05 1692346.0 1940.0 0.0 (200310, 1200) (200310, 1200) False 67.0 0.0
computeMatrix cm_triticum_6 True 166398474.0 1.0000000000000002 0.0 1156608.0 0.993049167025414 0.0007810000024619 row1809 col59 6.230529601420289e-08 1074690.0 0.0 0.0 (146597, 1200) (146597, 1200) False 7931.0 0.0
Binary file added benchmark/results/numerical_diff_hist.pdf
View file
Open in desktop
Binary file not shown.
Binary file added benchmark/results/numerical_diff_hist.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added benchmark/results/numerical_diff_hist.tiff
View file
Open in desktop
Binary file not shown.
20 changes: 13 additions & 7 deletions benchmark/rules/benchmark_summary.smk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,9 @@ rule combine_numerical_diff:
input:
expand("results/numerical_diff/alignmentsieve_{s}.tsv", s=alignmentSieve_samples.keys()),
expand("results/numerical_diff/bamcoverage_{s}.tsv", s=bamCoverage_samples.keys()),
expand("results/numerical_diff/bamcompare_{s}.tsv", s=bamCompare_samples.keys()),
expand("results/numerical_diff/multibamsummary_{s}.tsv", s=multibamSummary_samples.keys()),
expand("results/numerical_diff/computematrix_{s}.tsv", s=computeMatrix_samples.keys())
expand("results/numerical_diff/bamcompare_{s}.tsv", s=VERIFY_bamCompare_samples.keys()),
expand("results/numerical_diff/multibamsummary_{s}.tsv", s=VERIFY_multibamSummary_samples.keys()),
expand("results/numerical_diff/computematrix_{s}.tsv", s=VERIFY_computeMatrix_samples.keys())
output:
tsv = 'results/numerical_diff.tsv'
run:
Expand All @@ -40,10 +40,16 @@ rule plot_numerical_diff:
tsv = 'results/numerical_diff.tsv',
hists = lambda wc: (
expand("results/numerical_diff/bamcoverage_{s}.hist.tsv", s=bamCoverage_samples.keys())
+ expand("results/numerical_diff/bamcompare_{s}.hist.tsv", s=bamCompare_samples.keys())
+ expand("results/numerical_diff/multibamsummary_{s}.hist.tsv", s=multibamSummary_samples.keys())
+ expand("results/numerical_diff/computematrix_{s}.hist.tsv", s=computeMatrix_samples.keys())
)
+ expand("results/numerical_diff/bamcompare_{s}.hist.tsv", s=VERIFY_bamCompare_samples.keys())
+ expand("results/numerical_diff/multibamsummary_{s}.hist.tsv", s=VERIFY_multibamSummary_samples.keys())
+ expand("results/numerical_diff/computematrix_{s}.hist.tsv", s=VERIFY_computeMatrix_samples.keys())
),
relhists = lambda wc: (
expand("results/numerical_diff/bamcoverage_{s}.relhist.tsv", s=bamCoverage_samples.keys())
+ expand("results/numerical_diff/bamcompare_{s}.relhist.tsv", s=VERIFY_bamCompare_samples.keys())
+ expand("results/numerical_diff/multibamsummary_{s}.relhist.tsv", s=VERIFY_multibamSummary_samples.keys())
+ expand("results/numerical_diff/computematrix_{s}.relhist.tsv", s=VERIFY_computeMatrix_samples.keys())
),
output:
png = 'results/numerical_diff_hist.png',
pdf = 'results/numerical_diff_hist.pdf',
Expand Down
53 changes: 41 additions & 12 deletions benchmark/rules/scripts/compare_dt3_dt4.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
import re

import numpy as np
import pandas as pd

Expand All @@ -8,7 +10,16 @@
row = {"tool": tool, "sample": sample}

CHUNK_SIZE = 50_000_000
CANONICAL = re.compile(r"^(chr)?(\d{1,2}|X|Y|M|MT)$|^[1-7][ABD]$")


def is_canonical(name):
return CANONICAL.match(str(name)) is not None


row["canonical_only"] = True
HIST_EDGES = np.logspace(-9, 4, 13 * 4 + 1)
REL_EDGES = np.logspace(-12, 0, 12 * 4 + 1) # |dt3 - dt4| / max(|dt3|, |dt4|) is at most 1
THRESHOLDS = (1e-6, 1e-3, 1e-1)


Expand All @@ -23,6 +34,7 @@ def __init__(self):
self.max_abs_loc = None
self.max_rel = 0.0
self.counts = np.zeros(len(HIST_EDGES) + 1, dtype=np.int64) # under, bins..., over
self.rel_counts = np.zeros(len(REL_EDGES) + 1, dtype=np.int64)
self.n_zero = 0
self.n_gt = {t: 0 for t in THRESHOLDS}

Expand Down Expand Up @@ -61,8 +73,10 @@ def update(self, x, y, loc=lambda flat_idx: flat_idx):
if n_nz:
dn = d[nz]
denom = np.maximum(np.abs(x[nz]), np.abs(y[nz]))
self.max_rel = max(self.max_rel, float((dn / denom).max()))
rel = dn / denom
self.max_rel = max(self.max_rel, float(rel.max()))
self.counts += np.bincount(np.searchsorted(HIST_EDGES, dn), minlength=len(self.counts))
self.rel_counts += np.bincount(np.searchsorted(REL_EDGES, rel), minlength=len(self.rel_counts))
for t in THRESHOLDS:
self.n_gt[t] += int((dn > t).sum())

Expand Down Expand Up @@ -90,24 +104,27 @@ def finish(self):
out[f"n_diff_gt_{t:g}"] = self.n_gt[t]
return out

def histogram(self):
lo = np.concatenate([[0.0], [0.0], HIST_EDGES[:-1], [HIST_EDGES[-1]]])
hi = np.concatenate([[0.0], [HIST_EDGES[0]], HIST_EDGES[1:], [np.inf]])
counts = np.concatenate([[self.n_zero], self.counts])
def histogram(self, relative=False):
edges, counts = (REL_EDGES, self.rel_counts) if relative else (HIST_EDGES, self.counts)
lo = np.concatenate([[0.0], [0.0], edges[:-1], [edges[-1]]])
hi = np.concatenate([[0.0], [edges[0]], edges[1:], [np.inf]])
counts = np.concatenate([[self.n_zero], counts])
return pd.DataFrame({"lower": lo, "upper": hi, "count": counts})


def write_hist(stats):
if "hist" in snakemake.output.keys():
stats.histogram().to_csv(snakemake.output.hist, sep="\t", index=False)
if "relhist" in snakemake.output.keys():
stats.histogram(relative=True).to_csv(snakemake.output.relhist, sep="\t", index=False)


if mode == "bigwig":
import pyBigWig

def paired_chunks(bw3, bw4):
for chrom, length in bw3.chroms().items():
if chrom not in bw4.chroms():
if chrom not in bw4.chroms() or not is_canonical(chrom):
continue
for start in range(0, length, CHUNK_SIZE):
end = min(start + CHUNK_SIZE, length)
Expand Down Expand Up @@ -142,6 +159,9 @@ def paired_chunks(bw3, bw4):
d.columns = ["chrom", "start", "end"] + [f"s{i}" for i in range(d.shape[1] - 3)]
d["chrom"] = d["chrom"].astype(str)
merged = d3.merge(d4, on=["chrom", "start", "end"], suffixes=("_3", "_4"))
canon = merged["chrom"].map(is_canonical)
row["n_noncanonical_bins_excluded"] = int((~canon).sum())
merged = merged[canon]
ncol = d3.shape[1] - 3
v3 = merged[[f"s{i}_3" for i in range(ncol)]].to_numpy(dtype=np.float64)
v4 = merged[[f"s{i}_4" for i in range(ncol)]].to_numpy(dtype=np.float64)
Expand Down Expand Up @@ -179,12 +199,20 @@ def loc(i):

n_rows = min(m3.shape[0], m4.shape[0])
ncol = min(m3.shape[1], m4.shape[1])
# rows are regions; keep canonical chromosomes only, and check both versions list the same regions
r3, r4 = hm3.matrix.regions[:n_rows], hm4.matrix.regions[:n_rows]
keep = np.array([is_canonical(r[0]) for r in r3], dtype=bool)
row["n_noncanonical_rows_excluded"] = int((~keep).sum())
row["n_region_mismatch"] = sum(
(a[0], a[1], a[2]) != (b[0], b[1], b[2]) for a, b, k in zip(r3, r4, keep) if k
)
rows = np.flatnonzero(keep)
stats = DiffStats()
for a in range(0, n_rows, ROW_CHUNK):
b = min(a + ROW_CHUNK, n_rows)
c3 = np.ma.filled(m3[a:b, :ncol].astype(np.float64), np.nan)
c4 = np.ma.filled(m4[a:b, :ncol].astype(np.float64), np.nan)
stats.update(c3, c4, loc=lambda i, a=a: f"row{a + i // ncol} col{i % ncol}")
for a in range(0, len(rows), ROW_CHUNK):
idx = rows[a:a + ROW_CHUNK]
c3 = np.ma.filled(m3[idx, :ncol].astype(np.float64), np.nan)
c4 = np.ma.filled(m4[idx, :ncol].astype(np.float64), np.nan)
stats.update(c3, c4, loc=lambda i, idx=idx: f"row{idx[i // ncol]} col{i % ncol}")

row.update(stats.finish())
write_hist(stats)
Expand Down Expand Up @@ -268,7 +296,8 @@ def compare_contig(args):
with pysam.AlignmentFile(sorted3, "rb") as f3, pysam.AlignmentFile(sorted4, "rb") as f4:
if f3.references != f4.references:
raise RuntimeError("dt3 and dt4 BAMs have different reference sequences")
order = sorted(zip(f3.references, f3.lengths), key=lambda t: -t[1])
order = sorted(((n, l) for n, l in zip(f3.references, f3.lengths) if is_canonical(n)),
key=lambda t: -t[1])
jobs = [(sorted3, sorted4, name) for name, _ in order]

total = dict.fromkeys(COUNTERS, 0)
Expand Down
Loading
Loading