OSCR

Chromosome-scale genome assembly of the European common cuttlefish <i>Sepia officinalis</i>.

Code ↔ Paper

40 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 40 matches · 7 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ analysis/genome_alignments/breakpoints/breakpoint_analysis.ipynb, lines 81–128 · score 0.95 · DToL_31, DToL_41, DToL_44, DToL_5, MPIBR_2, MPIBR_40
  2. [2] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ analysis/genome_alignments/wga/assembly_comparison.R, lines 176–234 · score 0.94 · DToL_31, DToL_41, DToL_44, DToL_5, MPIBR_2, MPIBR_40
  3. [3] § Materials and methods › Nuclear genome annotation ↔ analysis/gene_families/cafe5/large_families_heatmap.R, lines 50–126 · score 0.92 · Pecten maximus, Nautilus pompilius, Octopus vulgaris, Octopus bimaculoides, Doryteuthis pealeii, Euprymna scolopes
  4. [4] § Materials and methods › Nuclear genome annotation ↔ annotation/task_braker_softmasked.sh, the whole file · a weak match · score 0.91 · addUTR, prot_seq, metazoa_odb10, BUSCO lineage, softmasked, BRAKER3
  5. [5] § Results › Gene modeling and annotation ↔ analysis/gene_families/cafe5/large_families_heatmap.R, lines 50–126 · score 0.89 · Pecten maximus, Nautilus pompilius, Octopus vulgaris, Octopus bimaculoides, Doryteuthis pealeii, Euprymna scolopes
  6. [6] § Materials and methods › Bulk RNA-seq analysis ↔ alignment/task_star_rna.sh, lines 26–89 · score 0.85 · RemoveNoncanonical, SortedByCoordinate, outFilterIntronMotifs, outSAMmultNmax, outSAMtype, STAR
  7. [7] § Materials and methods › Bulk RNA-seq analysis ↔ analysis/rna_seq/bulkRNA_Deseq2_braker.Rmd, lines 200–273 · score 0.85 · log2 fold change, DESeq2, target tissue, Wald, apeglm, binary
  8. [8] § Materials and methods › Gene family expansion analysis ↔ analysis/gene_families/cafe5/gene_repeat_overlap.py, lines 271–346 · score 0.84 · avoid double, overlap fraction, CDS length, repeat intervals, gene families, intersect
  9. [9] § Materials and methods › Bulk RNA-seq analysis ↔ alignment/task_star_scell.sh, lines 23–91 · score 0.84 · RemoveNoncanonical, SortedByCoordinate, outFilterIntronMotifs, outSAMmultNmax, outSAMtype, STAR
  10. [10] § Materials and methods › Coverage analysis ↔ analysis/genome_alignments/breakpoints/breakpoint_analysis.ipynb, lines 1143–1258 · score 0.81 · breakpoint rate, intra scaffold, trans pair, DToL, contact rates, empirical
  11. [11] § Materials and methods › Nuclear genome assembly ↔ assembly/run_scaffolding.sh, lines 44–101 · score 0.80 · JBAT, telo, motif, tool, YAHS, mapping
  12. [12] § Materials and methods › Coverage analysis ↔ analysis/genome_alignments/breakpoints/count_hic_trans_background.sh, lines 1–49 · score 0.78 · intra scaffold pairs, trans pair, scaffold lengths, product, background, Mb2
  13. [13] § Materials and methods › Bulk RNA-seq analysis ↔ analysis/rna_seq/parse_interproscan_go.R, lines 328–409 · score 0.76 · DESeq2, clusterProfiler, DE genes, GO terms, enricher, background
  14. [14] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ analysis/genome_alignments/breakpoints/extract_hic_trans_pairs.sh, lines 1–68 · score 0.76 · breakpoint scaffold pairs, scaffold_40, scaffold length, pairtools, BP1, BP3
  15. [15] § Materials and methods › Bulk RNA-seq analysis ↔ analysis/gene_families/gene_family_expression_analysis.Rmd, lines 414–472 · score 0.73 · clusterProfiler, GO enrichment, GO annotations, GO terms, enricher, FDR
  16. [16] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ analysis/genome_alignments/breakpoints/count_hic_trans_background.sh, lines 1–49 · score 0.73 · intra chromosomal contacts, intra scaffold pairs, scaffold length, background, pairtools, HiC
  17. [17] § Materials and methods › Bulk RNA-seq analysis ↔ alignment/task_star_count_braker.sh, lines 44–93 · score 0.72 · FeatureCounts, countReadPairs, gene_id, STAR, exon
  18. [18] § Materials and methods › Bulk RNA-seq analysis ↔ analysis/rna_seq/parse_interproscan_go.R, lines 159–223 · score 0.72 · reduce noise, unique GO terms, generic, union, collapsed, enrichment
  19. [19] § Materials and methods › Coverage analysis ↔ analysis/genome_alignments/breakpoints/breakpoint_analysis.ipynb, lines 131–240 · score 0.71 · mapping quality threshold, count_coverage, pysam, querying, depth, alignments
  20. [20] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ analysis/genome_alignments/breakpoints/run_repeatmasker_junctions.sh, lines 67–144 · score 0.71 · breakpoint junctions, junction window, scaffold_40, scaffold ends, BP2, RepeatMasker
  21. [21] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ analysis/genome_alignments/breakpoints/breakpoint_analysis.ipynb, lines 1143–1258 · score 0.69 · intra scaffold, control scaffolds, background pairs, DToL, Contact rates, empirical
  22. [22] § Materials and methods › Gene family expansion analysis ↔ analysis/gene_families/cafe5/run_cafe5.sh, lines 45–100 · score 0.66 · Gamma models, log likelihood, Gene family, lambda, lowest, CAFE5
  23. [23] § Materials and methods › Bulk RNA-seq analysis ↔ analysis/rna_seq/parse_interproscan_go.R, lines 328–409 · score 0.64 · DESeq2, clusterProfiler, GO term, enricher, tissue, RNA
  24. [24] § Materials and methods › Bulk RNA-seq analysis ↔ analysis/gene_families/gene_family_expression_analysis.Rmd, lines 414–472 · score 0.64 · clusterProfiler, GO enrichment, GO term, enricher, gene family, VST
  25. [25] § Materials and methods › Coverage analysis ↔ analysis/genome_alignments/breakpoints/extract_hic_trans_pairs.sh, lines 1–68 · score 0.64 · uninvolved scaffold, trans pairs, UU, pairtools, breakpoints, Hi
  26. [26] § Materials and methods › Gene family expansion analysis ↔ analysis/gene_families/cafe5/large_families_heatmap.R, lines 1–48 · score 0.63 · OrthoFinder, species tree, ape, ultrametric, rooted, Orthogroups
  27. [27] § Results › Gene modeling and annotation ↔ analysis/rna_seq/emapper_interproscan_LUT.R, the whole file · a weak match · score 0.63 · officinalis proteome, eggNOG, InterProScan, mapper, transcriptomic, filters
  28. [28] § Materials and methods › Bulk RNA-seq analysis ↔ alignment/task_countreads.sh, the whole file · a weak match · score 0.61 · featureCounts, countReadPairs, Subread, GTF, exon, id
  29. [29] § Materials and methods › Short-read RNA library preparation and sequencing ↔ analysis/rna_seq/bulkRNA_Deseq2_braker.Rmd, lines 48–65 · score 0.60 · NextSeq2000, library prep, mRNA, RNA seq, stranded, Illumina
  30. [30] § Materials and methods › Gene family expansion analysis ↔ analysis/gene_families/orthofinder/orthofinder_to_cafe.R, the whole file · a weak match · score 0.60 · OrthoFinder, species tree, ultrametric, rooted, proteomes, Orthogroups
  31. [31] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ analysis/genome_alignments/breakpoints/breakpoint_analysis.ipynb, lines 739–748 · score 0.59 · Hi Fi, DToL, kb, boundary, windows, breakpoint
  32. [32] § Materials and methods › Gene family expansion analysis ↔ analysis/gene_families/gene_family_expression_analysis.Rmd, lines 722–755 · score 0.57 · birth death model, extreme copy, Gene families, likelihood, orthogroups, species
  33. [33] § Materials and methods › Nuclear genome annotation ↔ analysis/rna_seq/emapper_interproscan_LUT.R, the whole file · a weak match · score 0.56 · eggNOG, InterProScan, lookup, mapper, GO, officinalis
  34. [34] § Materials and methods › Coverage analysis ↔ analysis/genome_alignments/breakpoints/breakpoint_analysis.ipynb, lines 1488–1527 · score 0.56 · LTR, RepeatMasker, SINE, parsed, classes, unknown
  35. [35] § Results › Genome size and heterozygosity ↔ assembly/run_busco.sh, the whole file · a weak match · score 0.55 · mollusca_odb12, metazoa_odb12, BUSCO, assembly, genome
  36. [36] § Materials and methods › Nuclear genome annotation ↔ annotation/task_repeatmasker_custom_softmask.sh, lines 47–77 · score 0.55 · RepeatModeler, xsmall, gff, softmasked, genome
  37. [37] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ analysis/genome_alignments/breakpoints/extract_spanning_reads.py, lines 1–37 · score 0.53 · Hi Fi, terminal, breakpoint, boundary, windows, alignments
  38. [38] § Materials and methods › Nuclear genome annotation ↔ analysis/genome_alignments/breakpoints/run_repeatmasker_junctions.sh, lines 146–198 · score 0.51 · repeat library, RepeatMasker, xsmall, gff, genome
  39. [39] § Results › Tissue-specific expression of expanded gene families ↔ analysis/rna_seq/bulkRNA_Deseq2_braker.Rmd, lines 83–117 · score 0.50 · Bulk RNA seq, DESeq2, posterior, vertical, retina, skin
  40. [40] § Results › Comparison with another chromosome-scale Sepia officinalis assembly ↔ assembly/run_busco.sh, the whole file · a weak match · score 0.50 · mollusca_odb12, metazoa_odb12, BUSCO, assemblies

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,598 lines · 67 KB · MIT · 6 matches

  1. # %% [markdown]
  2. # # Sepia officinalis Genome Assembly — Breakpoint Analysis
  3. # Jupyter notebook consolidating all visualisation and statistical analysis scripts for the MPIBR vs DToL assembly comparison.
  4. #
  5. # **Sections**
  6. # 1. HiFi spanning reads + trans HiC contacts at breakpoints
  7. # 2. Trans HiC contact rate statistics
  8. # 3. Repeat content at junction windows
  9. #
  10. # ---
  12. # %% [markdown]
  13. # ## Shared imports and plot settings
  15. # %%
  16. import os, sys, warnings
  17. import numpy as np
  18. import pandas as pd
  19. import matplotlib
  20. import matplotlib.pyplot as plt
  21. import matplotlib.patches as mpatches
  22. from matplotlib.patches import ConnectionPatch
  23. from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec
  24. from matplotlib.ticker import MaxNLocator
  25. from scipy import stats
  27. warnings.filterwarnings("ignore")
  29. matplotlib.rcParams.update({
  30. "font.family": "sans-serif",
  31. "font.sans-serif": ["Arial", "Helvetica", "DejaVu Sans"],
  32. "font.size": 7,
  33. "axes.linewidth": 0.5,
  34. "xtick.major.width": 0.5,
  35. "ytick.major.width": 0.5,
  36. "xtick.major.size": 2,
  37. "ytick.major.size": 2,
  38. "pdf.fonttype": 42,
  39. })
  41. try:
  42. import pysam
  43. except ImportError:
  44. sys.exit("ERROR: pysam is required — pip install pysam")
  46. print("Imports OK")
  48. # %% [markdown]
  49. # ---
  50. # ## 1. HiFi spanning reads + trans HiC contacts
  51. #
  52. # Visualises HiFi read alignments at scaffold junction ends alongside trans HiC contact dots and arcs.
  53. # Each breakpoint row shows: genome-wide contact density | terminal read track | HiFi coverage.
  54. # Control rows show terminal read track and coverage only.
  56. # %% [markdown]
  57. # ### 1.1 Parameters
  59. # %%
  60. # ── Input paths ─────────────────────────────────────────────────────────────
  61. BASE = "/gpfs/scic/data/projects/CuttlefishOmics/genome/assembly/soff250801"
  62. WD = "/gpfs/scic/data/projects/CuttlefishOmics/sandbox/resubmission"
  64. SPANNING_DIR = f"{WD}/spanning_reads"
  65. TRANS_HIC_DIR = f"{WD}/hic_trans_pairs"
  67. MPIBR_BAM = f"{BASE}/bams_hifi_scaffolds/soff250801_mpibr.hic_hifi_scaffolds.bam"
  68. DTOL_BAM = f"{BASE}/bams_hifi_scaffolds/soff250801_sanger.hic_hifi_scaffolds.bam"
  70. MPIBR_FAI = (f"{BASE}/yahs/soff250801_mpibr.hic/"
  71. "soff250801_mpibr.hic_scaffolds_final.fa.fai")
  72. DTOL_FAI = (f"{BASE}/yahs/soff250801_sanger.hic/"
  73. "soff250801_sanger.hic_scaffolds_final.fa.fai")
  75. # ── Run options ──────────────────────────────────────────────────────────────
  76. VIEW_WINDOW = 200_000 # bp shown in terminal read / coverage tracks
  77. MIN_MAPQ = 10
  78. OUTPUT_SPANNING = "spanning_reads.pdf"
  81. # %% [markdown]
  82. # ### 1.2 Constants, breakpoint definitions, colours
  84. # %%
  86. # ── Font sizes (edit here to scale all labels together) ─────────────
  87. FONT = {
  88. "title": 10, # scaffold name on density track
  89. "label": 9, # axis labels (ylabel, xlabel)
  90. "tick": 8, # tick labels
  91. "annot": 7, # small annotations ("200 kb view" etc.)
  92. "gap": 9, # gap-column text (pair counts, spanning reads)
  93. "legend": 9, # figure legend
  94. "boundary": 7, # "scaffold end / start" rotated label
  95. }
  97. MB = 1_000_000
  98. COV_BIN = 1000
  99. DENS_BIN = 500_000
  100. DOT_Y = -0.7
  101. MAX_ARCS = 200
  103. BREAKPOINTS_SPANNING = [
  104. ("mpibr", "scaffold_40", "scaffold_44", "BP1_MPIBR_40-44", "MPIBR_40", "MPIBR_44"),
  105. ("dtol", "scaffold_31", "scaffold_40", "BP2_DToL_31-40", "DToL_31", "DToL_40"),
  106. ("dtol", "scaffold_41", "scaffold_46", "BP3_DToL_41-46", "DToL_41", "DToL_46"),
  107. ("dtol", "scaffold_44", "scaffold_45", "BP4_DToL_44-45", "DToL_44", "DToL_45"),
  108. ]
  110. # Correct single scaffold from the other assembly for each breakpoint
  111. # (the chromosome that the broken assembly should have kept intact)
  112. CORRECT_SCAFFOLDS = [
  113. ("dtol", "scaffold_5", "BP1", "DToL_5"), # BP1: MPIBR split → DToL kept intact
  114. ("mpibr", "scaffold_2", "BP2", "MPIBR_2"), # BP2: DToL split → MPIBR kept intact
  115. ("mpibr", "scaffold_6", "BP3", "MPIBR_6"), # BP3: DToL split → MPIBR kept intact
  116. ("mpibr", "scaffold_7", "BP4", "MPIBR_7"), # BP4: DToL split → MPIBR kept intact
  117. ]
  119. COLOURS_SPANNING = {
  120. "mpibr": "#2166AC",
  121. "dtol": "#e65c00",
  122. "spanning": "#2CA02C",
  123. "background": "#CCCCCC",
  124. "hic_contact": "#9467BD",
  125. "boundary": "#444444",
  126. "coverage": {"mpibr": "#AEC6E8", "dtol": "#F4A582"},
  127. "density": {"mpibr": "#6BAED6", "dtol": "#FC8D59"},
  128. }
  131. # %% [markdown]
  132. # ### 1.3 Data-loading helpers
  134. # %%
  135. def load_fai(fai_path):
  136. lengths = {}
  137. with open(fai_path) as fh:
  138. for line in fh:
  139. parts = line.strip().split("\t")
  140. if len(parts) >= 2:
  141. lengths[parts[0]] = int(parts[1])
  142. return lengths
  145. def find_trans_tsv(trans_dir, label, assembly):
  146. return os.path.join(trans_dir, f"{label}_{assembly}_trans.tsv")
  149. def find_ctrl_trans_tsv(trans_dir, assembly, ctrl_scaf):
  150. ctrl_n = ctrl_scaf.replace("scaffold_", "")
  151. prefix = f"Ctrl_{assembly.upper()}_{ctrl_n}-"
  152. for fname in os.listdir(trans_dir):
  153. if fname.startswith(prefix) and fname.endswith("_trans.tsv"):
  154. return os.path.join(trans_dir, fname)
  155. return ""
  158. def load_trans_pairs(tsv_path):
  159. if not tsv_path or not os.path.isfile(tsv_path):
  160. return pd.DataFrame()
  161. df = pd.read_csv(tsv_path, sep="\t")
  162. return df if not df.empty else pd.DataFrame()
  165. def compute_density(positions, scaffold_len, bin_size=DENS_BIN):
  166. n_bins = max(1, scaffold_len // bin_size)
  167. counts, edges = np.histogram(positions, bins=n_bins, range=(0, scaffold_len))
  168. mids = (edges[:-1] + edges[1:]) / 2
  169. return mids, counts.astype(float)
  172. def normalised_contact_rate(n_pairs, len_A, len_B):
  173. la, lb = len_A / MB, len_B / MB
  174. return 0.0 if la == 0 or lb == 0 else n_pairs / (la * lb)
  177. def get_reads_in_window(bam_path, scaffold, win_start, win_end,
  178. spanning_names, min_mapq=10):
  179. reads = []
  180. try:
  181. bam = pysam.AlignmentFile(bam_path, "rb")
  182. for read in bam.fetch(scaffold, win_start, win_end):
  183. if read.is_unmapped or read.is_secondary:
  184. continue
  185. if read.mapping_quality < min_mapq:
  186. continue
  187. reads.append({
  188. "read_name": read.query_name,
  189. "start": max(read.reference_start, win_start),
  190. "end": min(read.reference_end, win_end),
  191. "strand": "-" if read.is_reverse else "+",
  192. "is_spanning": read.query_name in spanning_names,
  193. })
  194. bam.close()
  195. except (ValueError, KeyError) as e:
  196. print(f" WARNING: fetch failed {scaffold}:{win_start}-{win_end}: {e}")
  197. return reads
  200. def get_coverage(bam_path, scaffold, win_start, win_end,
  201. min_mapq=10, bin_size=COV_BIN):
  202. try:
  203. bam = pysam.AlignmentFile(bam_path, "rb")
  204. cov_arrays = bam.count_coverage(
  205. scaffold, win_start, win_end,
  206. quality_threshold=0,
  207. read_callback=lambda r: (
  208. not r.is_unmapped and not r.is_secondary and
  209. r.mapping_quality >= min_mapq
  210. )
  211. )
  212. bam.close()
  213. except (ValueError, KeyError) as e:
  214. print(f" WARNING: coverage failed {scaffold}:{win_start}-{win_end}: {e}")
  215. n = (win_end - win_start) // bin_size + 1
  216. return np.zeros(n), np.zeros(n)
  217. depth = sum(np.array(a) for a in cov_arrays)
  218. n_bins = len(depth) // bin_size
  219. if n_bins == 0:
  220. return np.array([len(depth) / 2]), np.array([depth.mean()])
  221. binned = depth[:n_bins * bin_size].reshape(n_bins, bin_size).mean(axis=1)
  222. mids = np.arange(n_bins) * bin_size + bin_size / 2
  223. return mids, binned
  226. def stack_reads(reads):
  227. reads_sorted = sorted(reads, key=lambda r: (not r["is_spanning"], r["start"]))
  228. track_ends = []
  229. for read in reads_sorted:
  230. placed = False
  231. for i, end in enumerate(track_ends):
  232. if read["start"] > end + 500:
  233. track_ends[i] = read["end"]
  234. read["track"] = i
  235. placed = True
  236. break
  237. if not placed:
  238. read["track"] = len(track_ends)
  239. track_ends.append(read["end"])
  240. return reads_sorted
  242. # %% [markdown]
  243. # ### 1.4 Drawing helpers
  245. # %%
  246. def draw_density_track(ax, mids, counts, scaffold_len, assembly,
  247. global_ymax, title, view_win, boundary_x,
  248. show_boundary_side, view_win_other=None):
  249. colour = COLOURS_SPANNING["density"][assembly]
  250. line_c = COLOURS_SPANNING[assembly]
  251. ax.fill_between(mids / MB, 0, counts,
  252. color=colour, alpha=0.7, linewidth=0, step="mid", zorder=1)
  253. ax.step(mids / MB, counts, color=line_c, linewidth=0.5, where="mid", zorder=2)
  254. # Other-end window (pale gold, drawn first)
  255. if view_win_other is not None:
  256. vw2_start, vw2_end = view_win_other
  257. ax.axvspan(vw2_start / MB, vw2_end / MB, color="gold", alpha=0.4, zorder=3, linewidth=0)
  258. inner_edge2 = vw2_end / MB if vw2_start == 0 else vw2_start / MB
  259. ax.axvline(inner_edge2, color="goldenrod", linewidth=0.8, linestyle="-", zorder=6)
  260. mid_vw2 = (vw2_start + vw2_end) / 2 / MB
  261. # ax.annotate("other\nend", xy=(mid_vw2, global_ymax * 1.05),
  262. # xycoords="data", ha="center", va="bottom",
  263. # fontsize=FONT["annot"], color="goldenrod", zorder=7, annotation_clip=False)
  264. # Junction-facing window (bright gold)
  265. vw_start, vw_end = view_win
  266. ax.axvspan(vw_start / MB, vw_end / MB, color="gold", alpha=0.9, zorder=3, linewidth=0)
  267. inner_edge = vw_start / MB if boundary_x > 0 else vw_end / MB
  268. ax.axvline(inner_edge, color="goldenrod", linewidth=1.2, linestyle="-", zorder=6)
  269. mid_vw = (vw_start + vw_end) / 2 / MB
  270. # ax.annotate("junction\nend", xy=(mid_vw, global_ymax * 1.05),
  271. # xycoords="data", ha="center", va="bottom",
  272. # fontsize=FONT["annot"], color="goldenrod", zorder=7, annotation_clip=False)
  273. ax.annotate("200 kb\nview", xy=(mid_vw, global_ymax * 1.05),
  274. xycoords="data", ha="center", va="bottom",
  275. fontsize=FONT["annot"], color="goldenrod", zorder=7, annotation_clip=False)
  276. ax.axvline(boundary_x / MB, color=COLOURS_SPANNING["boundary"],
  277. linewidth=1.0, linestyle="--", zorder=5)
  278. ax.set_xlim(0, scaffold_len / MB)
  279. ax.set_ylim(0, global_ymax * 1.08)
  280. ax.yaxis.set_major_locator(MaxNLocator(nbins=3, min_n_ticks=2, integer=True))
  281. ax.tick_params(axis="both", labelsize=FONT["tick"], pad=1)
  282. ax.tick_params(axis="x", labelbottom=True)
  283. ax.spines["top"].set_visible(False)
  284. ax.spines["right"].set_visible(False)
  285. ax.set_title(title, fontsize=FONT["title"], fontweight="bold",
  286. color=COLOURS_SPANNING[assembly], pad=3)
  289. def draw_read_track(ax, reads, win_start, win_end, trans_df, pos_col,
  290. read_height=0.6):
  291. n_tracks = max((r["track"] for r in reads), default=0) + 1
  292. for read in reads:
  293. x0 = read["start"] - win_start
  294. x1 = read["end"] - win_start
  295. colour = COLOURS_SPANNING["spanning"] if read["is_spanning"] else COLOURS_SPANNING["background"]
  296. alpha = 0.9 if read["is_spanning"] else 0.4
  297. zorder = 3 if read["is_spanning"] else 1
  298. ax.barh(read["track"], x1 - x0, left=x0, height=read_height,
  299. color=colour, alpha=alpha, zorder=zorder, linewidth=0)
  300. if not trans_df.empty:
  301. in_win = trans_df[
  302. (trans_df[pos_col] >= win_start) &
  303. (trans_df[pos_col] < win_end)].copy()
  304. if not in_win.empty:
  305. ax.scatter(in_win[pos_col] - win_start, [DOT_Y] * len(in_win),
  306. s=4, color=COLOURS_SPANNING["hic_contact"],
  307. alpha=0.6, zorder=4, linewidths=0, clip_on=False)
  308. ax.set_xlim(0, win_end - win_start)
  309. ax.set_ylim(DOT_Y - 0.5, n_tracks + 1)
  310. ax.set_yticks([])
  311. ax.set_xticks([])
  312. for spine in ["top", "right", "left", "bottom"]:
  313. ax.spines[spine].set_visible(False)
  314. return n_tracks
  317. def draw_hic_arcs(fig, ax_A, ax_B, trans_df, win_A, win_B):
  318. if trans_df.empty:
  319. return
  320. win_A_start, win_A_end = win_A
  321. win_B_start, win_B_end = win_B
  322. in_win = trans_df[
  323. (trans_df["pos1"] >= win_A_start) & (trans_df["pos1"] < win_A_end) &
  324. (trans_df["pos2"] >= win_B_start) & (trans_df["pos2"] < win_B_end)]
  325. if in_win.empty:
  326. return
  327. df_plot = in_win.sample(min(MAX_ARCS, len(in_win)), random_state=42)
  328. for _, row in df_plot.iterrows():
  329. con = ConnectionPatch(
  330. xyA=(row["pos1"] - win_A_start, DOT_Y), coordsA="data", axesA=ax_A,
  331. xyB=(row["pos2"] - win_B_start, DOT_Y), coordsB="data", axesB=ax_B,
  332. color=COLOURS_SPANNING["hic_contact"],
  333. linewidth=0.35, alpha=0.25, clip_on=False, zorder=5)
  334. fig.add_artist(con)
  337. def draw_coverage_track(ax, mids, depth, win_start, win_end, assembly, global_ymax):
  338. colour = COLOURS_SPANNING["coverage"][assembly]
  339. ax.fill_between(mids, 0, depth, color=colour, alpha=0.7, linewidth=0)
  340. ax.plot(mids, depth, color=COLOURS_SPANNING[assembly], linewidth=0.5)
  341. ax.set_xlim(0, win_end - win_start)
  342. ax.set_ylim(0, global_ymax * 1.08)
  343. ax.yaxis.set_major_locator(MaxNLocator(nbins=3, min_n_ticks=2))
  344. ax.tick_params(axis="y", labelsize=FONT["tick"], pad=1)
  345. ax.spines["top"].set_visible(False)
  346. ax.spines["right"].set_visible(False)
  347. tick_pos = np.linspace(0, win_end - win_start, 5)
  348. tick_lbl = [f"{(win_start + p) / MB:.2f}" for p in tick_pos]
  349. ax.set_xticks(tick_pos)
  350. ax.set_xticklabels(tick_lbl, fontsize=FONT["tick"])
  351. ax.set_xlabel("Position (Mb)", fontsize=FONT["label"], labelpad=2)
  354. def draw_boundary_line(axes, x):
  355. for ax in axes:
  356. ax.axvline(x, color=COLOURS_SPANNING["boundary"],
  357. linewidth=1.0, linestyle="--", zorder=10)
  360. def add_boundary_label(ax, x, side):
  361. text = "scaffold end" if side == "right" else "scaffold start"
  362. ax.text(x, ax.get_ylim()[1] * 0.92, text,
  363. ha=side, va="top", fontsize=FONT["boundary"],
  364. color=COLOURS_SPANNING["boundary"], rotation=90)
  366. # %% [markdown]
  367. # ### 1.5 Build row metadata and pre-compute coverage
  369. # %%
  370. lengths = {
  371. "mpibr": load_fai(MPIBR_FAI),
  372. "dtol": load_fai(DTOL_FAI),
  373. }
  374. bams = {"mpibr": MPIBR_BAM, "dtol": DTOL_BAM}
  376. ASM_DISPLAY = {"mpibr": "MPIBR", "dtol": "DToL"}
  377. # Correct-scaffold entries — one per breakpoint, matched in order
  378. controls = [
  379. (asm, scaf, None, f"Correct_{bp_lbl}", lbl, None, True)
  380. for asm, scaf, bp_lbl, lbl in CORRECT_SCAFFOLDS
  381. ]
  383. all_rows = BREAKPOINTS_SPANNING + controls
  384. n_rows = len(all_rows)
  386. print("Pre-computing HiFi coverage...")
  387. cov_store = {}
  389. def get_cov(assembly, scaffold, ws, we, bin_size=COV_BIN):
  390. key = (assembly, scaffold, ws, we)
  391. if key not in cov_store:
  392. mids, depth = get_coverage(bams[assembly], scaffold, ws, we,
  393. min_mapq=MIN_MAPQ, bin_size=bin_size)
  394. cov_store[key] = (mids, depth)
  395. print(f" {assembly} {scaffold}:{ws}-{we} mean={depth.mean():.1f}x")
  396. return cov_store[key]
  398. # ── Trans HiC diagnostic ───────────────────────────────────────────────
  399. _hic_dir_ok = os.path.isdir(TRANS_HIC_DIR)
  400. if _hic_dir_ok:
  401. _hic_files = [f for f in os.listdir(TRANS_HIC_DIR) if f.endswith('_trans.tsv')]
  402. print(f"Trans HiC dir: {TRANS_HIC_DIR} ({len(_hic_files)} TSV files)")
  403. else:
  404. print(f"WARNING: TRANS_HIC_DIR not found: {TRANS_HIC_DIR}")
  405. print(" → Run extract_hic_trans_pairs.sh first, then re-run this cell.")
  407. row_meta = []
  408. for entry in all_rows:
  409. is_ctrl = len(entry) == 7
  410. assembly = entry[0]
  411. scaf_A = entry[1]
  412. scaf_B = entry[2]
  413. label = entry[3]
  414. lbl_A = entry[4]
  415. lbl_B = entry[5]
  416. lens = lengths[assembly]
  417. len_A = lens.get(scaf_A, 0)
  418. end_A, end_B = "end", "start"
  419. spanning_names = set()
  421. if not is_ctrl:
  422. tsv = os.path.join(SPANNING_DIR, f"{label}_{assembly}.tsv")
  423. if os.path.isfile(tsv):
  424. sdf = pd.read_csv(tsv, sep="\t")
  425. spanning_names = set(sdf["read_name"]) if not sdf.empty else set()
  426. if not sdf.empty:
  427. a_rows = sdf[sdf["primary_scaffold"] == scaf_A]
  428. b_rows = sdf[sdf["primary_scaffold"] == scaf_B]
  429. a_supp = sdf[sdf["supp_scaffold"] == scaf_A]
  430. b_supp = sdf[sdf["supp_scaffold"] == scaf_B]
  431. if not a_rows.empty: end_A = a_rows["primary_end_label"].mode()[0]
  432. elif not a_supp.empty: end_A = a_supp["supp_end_label"].mode()[0]
  433. if not b_rows.empty: end_B = b_rows["primary_end_label"].mode()[0]
  434. elif not b_supp.empty: end_B = b_supp["supp_end_label"].mode()[0]
  436. win_A = (max(0, len_A - VIEW_WINDOW), len_A) if end_A == "end" else (0, min(VIEW_WINDOW, len_A))
  437. win_A_other = (0, min(VIEW_WINDOW, len_A)) if end_A == "end" else (max(0, len_A - VIEW_WINDOW), len_A)
  438. win_B = None
  439. win_B_other = None
  440. len_B = 0
  441. if not is_ctrl:
  442. len_B = lens.get(scaf_B, 0)
  443. win_B = (0, min(VIEW_WINDOW, len_B)) if end_B == "start" else (max(0, len_B - VIEW_WINDOW), len_B)
  444. win_B_other = (max(0, len_B - VIEW_WINDOW), len_B) if end_B == "start" else (0, min(VIEW_WINDOW, len_B))
  446. get_cov(assembly, scaf_A, *win_A)
  447. get_cov(assembly, scaf_A, *win_A_other) # both ends for all rows
  448. if win_B:
  449. get_cov(assembly, scaf_B, *win_B)
  450. get_cov(assembly, scaf_B, *win_B_other)
  452. trans_df = pd.DataFrame()
  453. len_partner = 0
  454. if is_ctrl:
  455. tsv_path = find_ctrl_trans_tsv(TRANS_HIC_DIR, assembly, scaf_A)
  456. if tsv_path:
  457. fname = os.path.basename(tsv_path)
  458. partner_n = fname.split("-")[1].split("_")[0]
  459. len_partner = lens.get(f"scaffold_{partner_n}", 0)
  460. trans_df = load_trans_pairs(tsv_path)
  461. if trans_df.empty:
  462. _expected = os.path.basename(tsv_path) if tsv_path else "(no TSV path)"
  463. print(f" WARNING: no trans HiC pairs loaded for {label} "
  464. f"-- file not found or empty: {_expected}")
  465. n_trans = len(trans_df)
  466. rate = normalised_contact_rate(n_trans, len_A, len_partner) \
  467. if len_partner > 0 else 0.0
  468. print(f" {label}: {n_trans} trans HiC pairs ({rate:.2f} pairs/Mb²)")
  469. else:
  470. tsv_path = find_trans_tsv(TRANS_HIC_DIR, label, assembly)
  471. len_partner = len_B
  472. trans_df = load_trans_pairs(tsv_path)
  473. if trans_df.empty:
  474. _expected = os.path.basename(tsv_path) if tsv_path else "(no TSV path)"
  475. print(f" WARNING: no trans HiC pairs loaded for {label} "
  476. f"-- file not found or empty: {_expected}")
  477. n_trans = len(trans_df)
  478. rate = normalised_contact_rate(n_trans, len_A, len_partner) \
  479. if len_partner > 0 else 0.0
  480. print(f" {label}: {n_trans} trans HiC pairs ({rate:.2f} pairs/Mb²)")
  482. row_meta.append({
  483. "assembly": assembly, "scaf_A": scaf_A, "scaf_B": scaf_B,
  484. "label": label, "lbl_A": lbl_A, "lbl_B": lbl_B,
  485. "is_ctrl": is_ctrl, "end_A": end_A, "end_B": end_B,
  486. "len_A": len_A, "len_B": len_B,
  487. "win_A": win_A, "win_A_other": win_A_other,
  488. "win_B": win_B, "win_B_other": win_B_other,
  489. "spanning_names": spanning_names, "trans_df": trans_df,
  490. "contact_rate": rate,
  491. })
  493. def global_cov_ymax(assembly):
  494. vals = [cov_store[k][1] for k in cov_store if k[0] == assembly]
  495. return float(np.quantile(np.concatenate(vals), 0.995)) if vals else 1.0
  497. cov_ymax = {"mpibr": global_cov_ymax("mpibr"), "dtol": global_cov_ymax("dtol")}
  499. dens_ymax = {"mpibr": 1.0, "dtol": 1.0}
  500. if os.path.isdir(TRANS_HIC_DIR):
  501. for asm in ["mpibr", "dtol"]:
  502. all_counts = []
  503. for meta in row_meta:
  504. if meta["assembly"] != asm or meta["trans_df"].empty:
  505. continue
  506. df = meta["trans_df"]
  507. if meta["len_A"] > 0:
  508. _, c = compute_density(df["pos1"].values, meta["len_A"], DENS_BIN)
  509. all_counts.extend(c.tolist())
  510. if meta["len_B"] > 0:
  511. _, c = compute_density(df["pos2"].values, meta["len_B"], DENS_BIN)
  512. all_counts.extend(c.tolist())
  513. if all_counts:
  514. dens_ymax[asm] = float(np.quantile(all_counts, 0.995))
  516. print("\nMetadata ready.")
  519. # %% [markdown]
  520. # ### 1.6 Generate plot
  522. # %%
  523. # ── Breakpoints figure ───────────────────────────────────────────────
  524. bp_meta = [m for m in row_meta if not m["is_ctrl"]]
  525. ctrl_meta = [m for m in row_meta if m["is_ctrl"]]
  527. def _draw_rows(fig, outer_gs, meta_list, bams, cov_store,
  528. cov_ymax, dens_ymax):
  529. """Draw one figure worth of rows; meta_list is bp_meta or ctrl_meta."""
  530. for row, meta in enumerate(meta_list):
  531. assembly = meta["assembly"]
  532. scaf_A = meta["scaf_A"]
  533. scaf_B = meta["scaf_B"]
  534. label = meta["label"]
  535. lbl_A = meta["lbl_A"]
  536. lbl_B = meta["lbl_B"]
  537. is_ctrl = meta["is_ctrl"]
  538. end_A = meta["end_A"]
  539. end_B = meta["end_B"]
  540. len_A = meta["len_A"]
  541. len_B = meta["len_B"]
  542. win_A_start, win_A_end = meta["win_A"]
  543. spanning_names = meta["spanning_names"]
  544. trans_df = meta["trans_df"]
  545. contact_rate = meta["contact_rate"]
  546. bam = bams[assembly]
  547. ymax_cov = cov_ymax[assembly]
  548. ymax_dens = dens_ymax[assembly]
  550. if is_ctrl:
  551. inner_A = GridSpecFromSubplotSpec(
  552. 5, 1, subplot_spec=outer_gs[row, 0],
  553. height_ratios=[2, 2, 1, 2, 1], hspace=0.55)
  554. else:
  555. inner_A = GridSpecFromSubplotSpec(
  556. 5, 1, subplot_spec=outer_gs[row, 0],
  557. height_ratios=[2, 2, 1, 2, 1], hspace=0.55)
  558. inner_B = GridSpecFromSubplotSpec(
  559. 5, 1, subplot_spec=outer_gs[row, 2],
  560. height_ratios=[2, 2, 1, 2, 1], hspace=0.55)
  562. if is_ctrl:
  563. win_A_o_start, win_A_o_end = meta["win_A_other"]
  565. ax_dens_A = fig.add_subplot(inner_A[0])
  566. ax_reads_end = fig.add_subplot(inner_A[1])
  567. ax_cov_end = fig.add_subplot(inner_A[2])
  568. ax_reads_str = fig.add_subplot(inner_A[3])
  569. ax_cov_str = fig.add_subplot(inner_A[4])
  571. # Density track (full scaffold, same as breakpoints)
  572. if not trans_df.empty and len_A > 0:
  573. mids_d, counts_d = compute_density(trans_df["pos1"].values, len_A, DENS_BIN)
  574. else:
  575. mids_d, counts_d = np.array([len_A / 2]), np.array([0.0])
  576. bnd_genome_A = len_A if end_A == "end" else 0
  577. draw_density_track(
  578. ax_dens_A, mids_d, counts_d, len_A, assembly, ymax_dens,
  579. title=lbl_A, view_win=meta["win_A"], boundary_x=bnd_genome_A,
  580. show_boundary_side="right" if end_A == "end" else "left",
  581. view_win_other=meta["win_A_other"])
  582. ax_dens_A.set_ylabel("Trans HiC\ncontacts", fontsize=FONT["label"], labelpad=3)
  585. # End terminal (scaffold end — win_A)
  586. reads_end = stack_reads(get_reads_in_window(
  587. bam, scaf_A, win_A_start, win_A_end, spanning_names, MIN_MAPQ))
  588. mids_end, dep_end = cov_store[(assembly, scaf_A, win_A_start, win_A_end)]
  589. draw_read_track(ax_reads_end, reads_end, win_A_start, win_A_end,
  590. trans_df, "pos1")
  591. draw_coverage_track(ax_cov_end, mids_end, dep_end,
  592. win_A_start, win_A_end, assembly, ymax_cov)
  593. ax_cov_end.set_ylabel("Depth (x)", fontsize=FONT["label"], labelpad=3)
  594. ax_cov_end.set_xlabel("")
  595. ax_cov_end.tick_params(axis="x", labelbottom=False)
  596. draw_boundary_line([ax_reads_end, ax_cov_end], win_A_end - win_A_start)
  598. # Start terminal (scaffold start — win_A_other)
  599. reads_str = stack_reads(get_reads_in_window(
  600. bam, scaf_A, win_A_o_start, win_A_o_end, spanning_names, MIN_MAPQ))
  601. mids_str, dep_str = cov_store[(assembly, scaf_A, win_A_o_start, win_A_o_end)]
  602. draw_read_track(ax_reads_str, reads_str, win_A_o_start, win_A_o_end,
  603. trans_df, "pos1")
  604. draw_coverage_track(ax_cov_str, mids_str, dep_str,
  605. win_A_o_start, win_A_o_end, assembly, ymax_cov)
  606. ax_cov_str.set_ylabel("Depth (x)", fontsize=FONT["label"], labelpad=3)
  607. draw_boundary_line([ax_reads_str, ax_cov_str], 0)
  609. # Turn off right panel
  610. ax_right = fig.add_subplot(outer_gs[row, 2])
  611. ax_right.set_axis_off()
  613. else:
  614. reads_A_j = stack_reads(get_reads_in_window(
  615. bam, scaf_A, win_A_start, win_A_end, spanning_names, MIN_MAPQ))
  616. mids_A_j, depth_A_j = cov_store[(assembly, scaf_A, win_A_start, win_A_end)]
  617. win_A_o_start, win_A_o_end = meta["win_A_other"]
  618. reads_A_o = stack_reads(get_reads_in_window(
  619. bam, scaf_A, win_A_o_start, win_A_o_end, spanning_names, MIN_MAPQ))
  620. mids_A_o, depth_A_o = cov_store[(assembly, scaf_A, win_A_o_start, win_A_o_end)]
  622. ax_dens_A = fig.add_subplot(inner_A[0])
  623. ax_reads_A_j = fig.add_subplot(inner_A[1])
  624. ax_cov_A_j = fig.add_subplot(inner_A[2])
  625. ax_reads_A_o = fig.add_subplot(inner_A[3])
  626. ax_cov_A_o = fig.add_subplot(inner_A[4])
  628. if not trans_df.empty and len_A > 0:
  629. mids_d, counts_d = compute_density(trans_df["pos1"].values, len_A, DENS_BIN)
  630. else:
  631. mids_d, counts_d = np.array([len_A / 2]), np.array([0.0])
  632. bnd_genome_A = len_A if end_A == "end" else 0
  633. draw_density_track(
  634. ax_dens_A, mids_d, counts_d, len_A, assembly, ymax_dens,
  635. title=lbl_A, view_win=meta["win_A"], boundary_x=bnd_genome_A,
  636. show_boundary_side="right" if end_A == "end" else "left",
  637. view_win_other=meta["win_A_other"])
  638. ax_dens_A.set_ylabel("Trans HiC\ncontacts", fontsize=FONT["label"], labelpad=3)
  640. draw_read_track(ax_reads_A_j, reads_A_j, win_A_start, win_A_end, trans_df, "pos1")
  641. draw_coverage_track(ax_cov_A_j, mids_A_j, depth_A_j,
  642. win_A_start, win_A_end, assembly, ymax_cov)
  643. ax_cov_A_j.set_ylabel("Depth (x)", fontsize=FONT["label"], labelpad=3)
  644. ax_cov_A_j.set_xlabel("")
  645. ax_cov_A_j.tick_params(axis="x", labelbottom=False)
  646. bnd_x_A_j = (win_A_end - win_A_start) if end_A == "end" else 0
  647. draw_boundary_line([ax_reads_A_j, ax_cov_A_j], bnd_x_A_j)
  649. end_A_other = "start" if end_A == "end" else "end"
  650. draw_read_track(ax_reads_A_o, reads_A_o, win_A_o_start, win_A_o_end, trans_df, "pos1")
  651. draw_coverage_track(ax_cov_A_o, mids_A_o, depth_A_o,
  652. win_A_o_start, win_A_o_end, assembly, ymax_cov)
  653. ax_cov_A_o.set_ylabel("Depth (x)", fontsize=FONT["label"], labelpad=3)
  654. bnd_x_A_o = 0 if end_A_other == "start" else (win_A_o_end - win_A_o_start)
  655. draw_boundary_line([ax_reads_A_o, ax_cov_A_o], bnd_x_A_o)
  657. ax_gap = fig.add_subplot(outer_gs[row, 1])
  658. ax_gap.set_axis_off()
  660. if is_ctrl:
  661. #bp_ref = meta["label"].replace("Correct_", "")
  662. # ax_gap.text(0.5, 0.75,
  663. # f"{bp_ref}\ncorrect\nassembly",
  664. # ha="center", va="center", fontsize=FONT["gap"],
  665. # color="grey", style="italic",
  666. # transform=ax_gap.transAxes)
  667. ax_gap.text(0.5, 0.45, f"{len(trans_df):,}\ntrans HiC\npairs",
  668. ha="center", va="center", fontsize=FONT["gap"],
  669. color=COLOURS_SPANNING["hic_contact"], fontweight="bold",
  670. transform=ax_gap.transAxes)
  671. ax_gap.text(0.5, 0.20, f"{contact_rate:.2f}\npairs/Mb\u00b2",
  672. ha="center", va="center", fontsize=FONT["gap"],
  673. color=COLOURS_SPANNING["hic_contact"],
  674. transform=ax_gap.transAxes)
  675. else:
  676. win_B_start, win_B_end = meta["win_B"]
  677. win_B_o_start, win_B_o_end = meta["win_B_other"]
  679. reads_B_j = stack_reads(get_reads_in_window(
  680. bam, scaf_B, win_B_start, win_B_end, spanning_names, MIN_MAPQ))
  681. mids_B_j, depth_B_j = cov_store[(assembly, scaf_B, win_B_start, win_B_end)]
  682. reads_B_o = stack_reads(get_reads_in_window(
  683. bam, scaf_B, win_B_o_start, win_B_o_end, spanning_names, MIN_MAPQ))
  684. mids_B_o, depth_B_o = cov_store[(assembly, scaf_B, win_B_o_start, win_B_o_end)]
  686. ax_dens_B = fig.add_subplot(inner_B[0])
  687. ax_reads_B_j = fig.add_subplot(inner_B[1])
  688. ax_cov_B_j = fig.add_subplot(inner_B[2])
  689. ax_reads_B_o = fig.add_subplot(inner_B[3])
  690. ax_cov_B_o = fig.add_subplot(inner_B[4])
  692. if not trans_df.empty and len_B > 0:
  693. mids_dB, counts_dB = compute_density(trans_df["pos2"].values, len_B, DENS_BIN)
  694. else:
  695. mids_dB, counts_dB = np.array([len_B / 2]), np.array([0.0])
  696. bnd_genome_B = 0 if end_B == "start" else len_B
  697. draw_density_track(
  698. ax_dens_B, mids_dB, counts_dB, len_B, assembly, ymax_dens,
  699. title=lbl_B, view_win=meta["win_B"], boundary_x=bnd_genome_B,
  700. show_boundary_side="left" if end_B == "start" else "right",
  701. view_win_other=meta["win_B_other"])
  702. ax_dens_B.set_yticklabels([])
  704. draw_read_track(ax_reads_B_j, reads_B_j, win_B_start, win_B_end, trans_df, "pos2")
  705. draw_coverage_track(ax_cov_B_j, mids_B_j, depth_B_j,
  706. win_B_start, win_B_end, assembly, ymax_cov)
  707. ax_cov_B_j.set_yticklabels([])
  708. ax_cov_B_j.set_xlabel("")
  709. ax_cov_B_j.tick_params(axis="x", labelbottom=False)
  710. bnd_x_B_j = 0 if end_B == "start" else (win_B_end - win_B_start)
  711. draw_boundary_line([ax_reads_B_j, ax_cov_B_j], bnd_x_B_j)
  713. end_B_other = "end" if end_B == "start" else "start"
  714. draw_read_track(ax_reads_B_o, reads_B_o, win_B_o_start, win_B_o_end, trans_df, "pos2")
  715. draw_coverage_track(ax_cov_B_o, mids_B_o, depth_B_o,
  716. win_B_o_start, win_B_o_end, assembly, ymax_cov)
  717. ax_cov_B_o.set_yticklabels([])
  718. bnd_x_B_o = (win_B_o_end - win_B_o_start) if end_B_other == "end" else 0
  719. draw_boundary_line([ax_reads_B_o, ax_cov_B_o], bnd_x_B_o)
  721. if not trans_df.empty:
  722. draw_hic_arcs(fig, ax_reads_A_j, ax_reads_B_j, trans_df,
  723. meta["win_A"], meta["win_B"])
  725. n_span = len(spanning_names)
  726. ax_gap.text(0.5, 0.88, f"{len(trans_df):,}\ntrans HiC\npairs",
  727. ha="center", va="center", fontsize=FONT["gap"],
  728. color=COLOURS_SPANNING["hic_contact"], fontweight="bold",
  729. transform=ax_gap.transAxes)
  730. ax_gap.text(0.5, 0.70, f"{contact_rate:.2f}\npairs/Mb\u00b2",
  731. ha="center", va="center", fontsize=FONT["gap"],
  732. color=COLOURS_SPANNING["hic_contact"], transform=ax_gap.transAxes)
  733. ax_gap.text(0.5, 0.50, f"{n_span}\nspanning\nread{'s' if n_span != 1 else ''}",
  734. ha="center", va="center", fontsize=FONT["gap"],
  735. color=COLOURS_SPANNING["spanning"], fontweight="bold",
  736. transform=ax_gap.transAxes)
  739. handles_leg = [
  740. mpatches.Patch(color=COLOURS_SPANNING["spanning"], alpha=0.9, label="HiFi spanning read"),
  741. mpatches.Patch(color=COLOURS_SPANNING["background"], alpha=0.5, label="HiFi non-spanning read"),
  742. mpatches.Patch(color=COLOURS_SPANNING["hic_contact"], alpha=0.7, label="Trans HiC contact"),
  743. mpatches.Patch(color=COLOURS_SPANNING["coverage"]["mpibr"], alpha=0.7, label="HiFi depth (MPIBR)"),
  744. mpatches.Patch(color=COLOURS_SPANNING["coverage"]["dtol"], alpha=0.7, label="HiFi depth (DToL)"),
  745. plt.Line2D([0], [0], color=COLOURS_SPANNING["boundary"], linewidth=0.9,
  746. linestyle="--", label="Scaffold boundary"),
  747. mpatches.Patch(color="yellow", alpha=0.4, label="200 kb view window"),
  748. ]
  750. # ── Breakpoints figure ────────────────────────────────────────────────
  751. n_bp = len(bp_meta)
  752. fig_bp = plt.figure(figsize=(18.0 / 2.54, 14)) #old height:(15.0 * n_bp + 2.5) / 2.54
  753. outer_gs_bp = GridSpec(
  754. n_bp, 3, figure=fig_bp,
  755. width_ratios=[1, 0.20, 1],
  756. height_ratios=[15.0] * n_bp,
  757. hspace=0.25, wspace=0.08,
  758. left=0.11, right=0.97, top=0.97, bottom=0.07,
  759. )
  760. _draw_rows(fig_bp, outer_gs_bp, bp_meta, bams, cov_store, cov_ymax, dens_ymax)
  761. fig_bp.legend(handles=handles_leg, loc="lower center",
  762. ncol=4, fontsize=FONT["legend"], frameon=False,
  763. bbox_to_anchor=(0.5, 0.0))
  764. fig_bp.savefig(OUTPUT_SPANNING, dpi=300, bbox_inches="tight")
  765. plt.show()
  766. print(f"Saved: {OUTPUT_SPANNING}")
  768. # ── Controls figure ───────────────────────────────────────────────────
  769. output_ctrl = OUTPUT_SPANNING.replace(".pdf", "_controls.pdf")
  770. n_ctrl = len(ctrl_meta)
  771. fig_ctrl = plt.figure(figsize=(4, 13))
  772. outer_gs_ctrl = GridSpec(
  773. n_ctrl, 3, figure=fig_ctrl,
  774. width_ratios=[1, 0.20, 0.20],
  775. height_ratios=[15.0] * n_ctrl,
  776. hspace=0.25, wspace=0.08,
  777. left=0.11, right=0.97, top=0.97, bottom=0.04,
  778. )
  779. _draw_rows(fig_ctrl, outer_gs_ctrl, ctrl_meta, bams, cov_store, cov_ymax, dens_ymax)
  780. fig_ctrl.savefig(output_ctrl, dpi=300, bbox_inches="tight")
  781. plt.show()
  782. print(f"Saved: {output_ctrl}")
  785. # %%
  786. # ── Breakpoints figure ───────────────────────────────────────────────
  787. bp_meta = [m for m in row_meta if not m["is_ctrl"]]
  788. ctrl_meta = [m for m in row_meta if m["is_ctrl"]]
  790. def _draw_rows(fig, outer_gs, meta_list, bams, cov_store,
  791. cov_ymax, dens_ymax):
  792. """Draw one figure worth of rows; meta_list is bp_meta or ctrl_meta."""
  793. for row, meta in enumerate(meta_list):
  794. assembly = meta["assembly"]
  795. scaf_A = meta["scaf_A"]
  796. scaf_B = meta["scaf_B"]
  797. label = meta["label"]
  798. lbl_A = meta["lbl_A"]
  799. lbl_B = meta["lbl_B"]
  800. is_ctrl = meta["is_ctrl"]
  801. end_A = meta["end_A"]
  802. end_B = meta["end_B"]
  803. len_A = meta["len_A"]
  804. len_B = meta["len_B"]
  805. win_A_start, win_A_end = meta["win_A"]
  806. spanning_names = meta["spanning_names"]
  807. trans_df = meta["trans_df"]
  808. contact_rate = meta["contact_rate"]
  809. bam = bams[assembly]
  810. ymax_cov = cov_ymax[assembly]
  811. ymax_dens = dens_ymax[assembly]
  813. if is_ctrl:
  814. inner_A = GridSpecFromSubplotSpec(
  815. 3, 1, subplot_spec=outer_gs[row, 0],
  816. height_ratios=[2, 3, 1], hspace=0.55)
  817. else:
  818. inner_A = GridSpecFromSubplotSpec(
  819. 5, 1, subplot_spec=outer_gs[row, 0],
  820. height_ratios=[2, 2, 1, 2, 1], hspace=0.55)
  821. inner_B = GridSpecFromSubplotSpec(
  822. 5, 1, subplot_spec=outer_gs[row, 2],
  823. height_ratios=[2, 2, 1, 2, 1], hspace=0.55)
  825. if is_ctrl:
  826. reads_A = stack_reads(get_reads_in_window(
  827. bam, scaf_A, win_A_start, win_A_end, spanning_names, MIN_MAPQ))
  828. mids_A, depth_A = cov_store[(assembly, scaf_A, win_A_start, win_A_end)]
  829. ax_dens_A = fig.add_subplot(inner_A[0])
  830. ax_reads_A = fig.add_subplot(inner_A[1])
  831. ax_cov_A = fig.add_subplot(inner_A[2])
  832. if not trans_df.empty and len_A > 0:
  833. mids_d, counts_d = compute_density(trans_df["pos1"].values, len_A, DENS_BIN)
  834. else:
  835. mids_d, counts_d = np.array([len_A / 2]), np.array([0.0])
  836. bnd_genome_A = len_A if end_A == "end" else 0
  837. draw_density_track(
  838. ax_dens_A, mids_d, counts_d, len_A, assembly, ymax_dens,
  839. title=lbl_A, view_win=meta["win_A"], boundary_x=bnd_genome_A,
  840. show_boundary_side="right" if end_A == "end" else "left")
  841. ax_dens_A.set_ylabel("Trans HiC\ncontacts", fontsize=FONT["label"], labelpad=3)
  842. draw_read_track(ax_reads_A, reads_A, win_A_start, win_A_end, trans_df, "pos1")
  843. draw_coverage_track(ax_cov_A, mids_A, depth_A,
  844. win_A_start, win_A_end, assembly, ymax_cov)
  845. ax_cov_A.set_ylabel("Depth (x)", fontsize=FONT["label"], labelpad=3)
  846. bnd_x_A = (win_A_end - win_A_start) if end_A == "end" else 0
  847. draw_boundary_line([ax_reads_A, ax_cov_A], bnd_x_A)
  849. else:
  850. reads_A_j = stack_reads(get_reads_in_window(
  851. bam, scaf_A, win_A_start, win_A_end, spanning_names, MIN_MAPQ))
  852. mids_A_j, depth_A_j = cov_store[(assembly, scaf_A, win_A_start, win_A_end)]
  853. win_A_o_start, win_A_o_end = meta["win_A_other"]
  854. reads_A_o = stack_reads(get_reads_in_window(
  855. bam, scaf_A, win_A_o_start, win_A_o_end, spanning_names, MIN_MAPQ))
  856. mids_A_o, depth_A_o = cov_store[(assembly, scaf_A, win_A_o_start, win_A_o_end)]
  858. ax_dens_A = fig.add_subplot(inner_A[0])
  859. ax_reads_A_j = fig.add_subplot(inner_A[1])
  860. ax_cov_A_j = fig.add_subplot(inner_A[2])
  861. ax_reads_A_o = fig.add_subplot(inner_A[3])
  862. ax_cov_A_o = fig.add_subplot(inner_A[4])
  864. if not trans_df.empty and len_A > 0:
  865. mids_d, counts_d = compute_density(trans_df["pos1"].values, len_A, DENS_BIN)
  866. else:
  867. mids_d, counts_d = np.array([len_A / 2]), np.array([0.0])
  868. bnd_genome_A = len_A if end_A == "end" else 0
  869. draw_density_track(
  870. ax_dens_A, mids_d, counts_d, len_A, assembly, ymax_dens,
  871. title=lbl_A, view_win=meta["win_A"], boundary_x=bnd_genome_A,
  872. show_boundary_side="right" if end_A == "end" else "left",
  873. view_win_other=meta["win_A_other"])
  874. ax_dens_A.set_ylabel("Trans HiC\ncontacts", fontsize=FONT["label"], labelpad=3)
  876. draw_read_track(ax_reads_A_j, reads_A_j, win_A_start, win_A_end, trans_df, "pos1")
  877. draw_coverage_track(ax_cov_A_j, mids_A_j, depth_A_j,
  878. win_A_start, win_A_end, assembly, ymax_cov)
  879. ax_cov_A_j.set_ylabel("Depth (x)", fontsize=FONT["label"], labelpad=3)
  880. ax_cov_A_j.set_xlabel("")
  881. ax_cov_A_j.tick_params(axis="x", labelbottom=True)
  882. bnd_x_A_j = (win_A_end - win_A_start) if end_A == "end" else 0
  883. draw_boundary_line([ax_reads_A_j, ax_cov_A_j], bnd_x_A_j)
  885. end_A_other = "start" if end_A == "end" else "end"
  886. draw_read_track(ax_reads_A_o, reads_A_o, win_A_o_start, win_A_o_end, trans_df, "pos1")
  887. draw_coverage_track(ax_cov_A_o, mids_A_o, depth_A_o,
  888. win_A_o_start, win_A_o_end, assembly, ymax_cov)
  889. ax_cov_A_o.set_ylabel("Depth (x)", fontsize=FONT["label"], labelpad=3)
  890. bnd_x_A_o = 0 if end_A_other == "start" else (win_A_o_end - win_A_o_start)
  891. draw_boundary_line([ax_reads_A_o, ax_cov_A_o], bnd_x_A_o)
  893. ax_gap = fig.add_subplot(outer_gs[row, 1])
  894. ax_gap.set_axis_off()
  896. if is_ctrl:
  897. ax_gap.text(0.5, 0.65, f"{len(trans_df):,}\ntrans HiC\npairs",
  898. ha="center", va="center", fontsize=FONT["gap"],
  899. color=COLOURS_SPANNING["hic_contact"], fontweight="bold",
  900. transform=ax_gap.transAxes)
  901. ax_gap.text(0.5, 0.28, f"{contact_rate:.2f}\npairs/Mb\u00b2",
  902. ha="center", va="center", fontsize=FONT["gap"],
  903. color=COLOURS_SPANNING["hic_contact"], transform=ax_gap.transAxes)
  904. else:
  905. win_B_start, win_B_end = meta["win_B"]
  906. win_B_o_start, win_B_o_end = meta["win_B_other"]
  908. reads_B_j = stack_reads(get_reads_in_window(
  909. bam, scaf_B, win_B_start, win_B_end, spanning_names, MIN_MAPQ))
  910. mids_B_j, depth_B_j = cov_store[(assembly, scaf_B, win_B_start, win_B_end)]
  911. reads_B_o = stack_reads(get_reads_in_window(
  912. bam, scaf_B, win_B_o_start, win_B_o_end, spanning_names, MIN_MAPQ))
  913. mids_B_o, depth_B_o = cov_store[(assembly, scaf_B, win_B_o_start, win_B_o_end)]
  915. ax_dens_B = fig.add_subplot(inner_B[0])
  916. ax_reads_B_j = fig.add_subplot(inner_B[1])
  917. ax_cov_B_j = fig.add_subplot(inner_B[2])
  918. ax_reads_B_o = fig.add_subplot(inner_B[3])
  919. ax_cov_B_o = fig.add_subplot(inner_B[4])
  921. if not trans_df.empty and len_B > 0:
  922. mids_dB, counts_dB = compute_density(trans_df["pos2"].values, len_B, DENS_BIN)
  923. else:
  924. mids_dB, counts_dB = np.array([len_B / 2]), np.array([0.0])
  925. bnd_genome_B = 0 if end_B == "start" else len_B
  926. draw_density_track(
  927. ax_dens_B, mids_dB, counts_dB, len_B, assembly, ymax_dens,
  928. title=lbl_B, view_win=meta["win_B"], boundary_x=bnd_genome_B,
  929. show_boundary_side="left" if end_B == "start" else "right",
  930. view_win_other=meta["win_B_other"])
  931. ax_dens_B.set_yticklabels([])
  933. draw_read_track(ax_reads_B_j, reads_B_j, win_B_start, win_B_end, trans_df, "pos2")
  934. draw_coverage_track(ax_cov_B_j, mids_B_j, depth_B_j,
  935. win_B_start, win_B_end, assembly, ymax_cov)
  936. ax_cov_B_j.set_yticklabels([])
  937. ax_cov_B_j.set_xlabel("")
  938. ax_cov_B_j.tick_params(axis="x", labelbottom=True)
  939. bnd_x_B_j = 0 if end_B == "start" else (win_B_end - win_B_start)
  940. draw_boundary_line([ax_reads_B_j, ax_cov_B_j], bnd_x_B_j)
  942. end_B_other = "end" if end_B == "start" else "start"
  943. draw_read_track(ax_reads_B_o, reads_B_o, win_B_o_start, win_B_o_end, trans_df, "pos2")
  944. draw_coverage_track(ax_cov_B_o, mids_B_o, depth_B_o,
  945. win_B_o_start, win_B_o_end, assembly, ymax_cov)
  946. ax_cov_B_o.set_yticklabels([])
  947. bnd_x_B_o = (win_B_o_end - win_B_o_start) if end_B_other == "end" else 0
  948. draw_boundary_line([ax_reads_B_o, ax_cov_B_o], bnd_x_B_o)
  950. if not trans_df.empty:
  951. draw_hic_arcs(fig, ax_reads_A_j, ax_reads_B_j, trans_df,
  952. meta["win_A"], meta["win_B"])
  954. n_span = len(spanning_names)
  955. ax_gap.text(0.5, 0.88, f"{len(trans_df):,}\ntrans HiC\npairs",
  956. ha="center", va="center", fontsize=FONT["gap"],
  957. color=COLOURS_SPANNING["hic_contact"], fontweight="bold",
  958. transform=ax_gap.transAxes)
  959. ax_gap.text(0.5, 0.70, f"{contact_rate:.2f}\npairs/Mb\u00b2",
  960. ha="center", va="center", fontsize=FONT["gap"],
  961. color=COLOURS_SPANNING["hic_contact"], transform=ax_gap.transAxes)
  962. ax_gap.text(0.5, 0.50, f"{n_span}\nspanning\nread{'s' if n_span != 1 else ''}",
  963. ha="center", va="center", fontsize=FONT["gap"],
  964. color=COLOURS_SPANNING["spanning"], fontweight="bold",
  965. transform=ax_gap.transAxes)
  968. handles_leg = [
  969. mpatches.Patch(color=COLOURS_SPANNING["spanning"], alpha=0.9, label="HiFi spanning read"),
  970. mpatches.Patch(color=COLOURS_SPANNING["background"], alpha=0.5, label="HiFi non-spanning read"),
  971. mpatches.Patch(color=COLOURS_SPANNING["hic_contact"], alpha=0.7, label="Trans HiC contact"),
  972. mpatches.Patch(color=COLOURS_SPANNING["coverage"]["mpibr"], alpha=0.7, label="HiFi depth (MPIBR)"),
  973. mpatches.Patch(color=COLOURS_SPANNING["coverage"]["dtol"], alpha=0.7, label="HiFi depth (DToL)"),
  974. plt.Line2D([0], [0], color=COLOURS_SPANNING["boundary"], linewidth=0.9,
  975. linestyle="--", label="Scaffold boundary"),
  976. mpatches.Patch(color="yellow", alpha=0.4, label="200 kb view window"),
  977. ]
  979. # ── Breakpoints figure ────────────────────────────────────────────────
  980. n_bp = len(bp_meta)
  981. fig_bp = plt.figure(figsize=(18.0 / 2.54, 14)) #old height:(15.0 * n_bp + 2.5) / 2.54
  982. outer_gs_bp = GridSpec(
  983. n_bp, 3, figure=fig_bp,
  984. width_ratios=[1, 0.20, 1],
  985. height_ratios=[15.0] * n_bp,
  986. hspace=0.25, wspace=0.08,
  987. left=0.11, right=0.97, top=0.97, bottom=0.07,
  988. )
  989. _draw_rows(fig_bp, outer_gs_bp, bp_meta, bams, cov_store, cov_ymax, dens_ymax)
  990. fig_bp.legend(handles=handles_leg, loc="lower center",
  991. ncol=4, fontsize=FONT["legend"], frameon=False,
  992. bbox_to_anchor=(0.5, 0.0))
  993. fig_bp.savefig(OUTPUT_SPANNING, dpi=300, bbox_inches="tight")
  994. plt.show()
  995. print(f"Saved: {OUTPUT_SPANNING}")
  997. # ── Controls figure ───────────────────────────────────────────────────
  998. output_ctrl = OUTPUT_SPANNING.replace(".pdf", "_controls.pdf")
  999. n_ctrl = len(ctrl_meta)
  1000. fig_ctrl = plt.figure(figsize=(6, 13))
  1001. outer_gs_ctrl = GridSpec(
  1002. n_ctrl, 3, figure=fig_ctrl,
  1003. width_ratios=[1, 0.20, 1],
  1004. height_ratios=[15.0] * n_ctrl,
  1005. hspace=0.25, wspace=0.08,
  1006. left=0.11, right=0.97, top=0.97, bottom=0.04,
  1007. )
  1008. _draw_rows(fig_ctrl, outer_gs_ctrl, ctrl_meta, bams, cov_store, cov_ymax, dens_ymax)
  1009. fig_ctrl.savefig(output_ctrl, dpi=300, bbox_inches="tight")
  1010. plt.show()
  1011. print(f"Saved: {output_ctrl}")
  1014. # %% [markdown]
  1015. # ---
  1016. # ## 2. Trans HiC contact rate statistics
  1017. #
  1018. # Tests whether the normalised trans contact rate (pairs per Mb²) between breakpoint scaffold pairs is significantly elevated above the background distribution of all trans scaffold pairs.
  1019. #
  1020. # - **Empirical p-value**: fraction of background pairs with rate ≥ observed
  1021. # - **Wilcoxon rank-sum** (one-tailed): DToL only (3 breakpoints)
  1022. # - **Plot A**: violin/strip plot of background, positive control, and breakpoint rates
  1023. # - **Plot B**: histogram of background distribution with breakpoint lines
  1025. # %% [markdown]
  1026. # ### 2.1 Parameters
  1028. # %%
  1029. MPIBR_TRANS = f"{WD}/hic_trans_background/mpibr_trans_counts.tsv"
  1030. MPIBR_INTRA = f"{WD}/hic_trans_background/mpibr_intra_counts.tsv"
  1031. DTOL_TRANS = f"{WD}/hic_trans_background/dtol_trans_counts.tsv"
  1032. DTOL_INTRA = f"{WD}/hic_trans_background/dtol_intra_counts.tsv"
  1034. OUTPUT_STATS = "hic_trans_stats" # prefix; .pdf and _distribution.pdf produced
  1037. # %% [markdown]
  1038. # ### 2.2 Breakpoint definitions and colours
  1040. # %%
  1041. BREAKPOINTS_STATS = {
  1042. "mpibr": [("BP1", "scaffold_40", "scaffold_44")],
  1043. "dtol": [
  1044. ("BP2", "scaffold_31", "scaffold_40"),
  1045. ("BP3", "scaffold_41", "scaffold_46"),
  1046. ("BP4", "scaffold_44", "scaffold_45"),
  1047. ],
  1048. }
  1050. COLOURS_STATS = {
  1051. "negative": "#AAAAAA", # grey — background trans
  1052. "positive": "#4DAF4A", # green — intra-scaffold
  1053. "breakpoint": {
  1054. "mpibr": "#2166AC", # blue — MPIBR breakpoints
  1055. "dtol": "#FF6600", # orange — DToL breakpoints
  1056. },
  1057. }
  1059. # %% [markdown]
  1060. # ### 2.3 Helper functions
  1062. # %%
  1063. def load_trans_bg(path):
  1064. df = pd.read_csv(path, sep="\t")
  1065. df.columns = df.columns.str.strip()
  1066. df["rate_per_mb2"] = pd.to_numeric(df["rate_per_mb2"], errors="coerce")
  1067. return df.dropna(subset=["rate_per_mb2"])
  1070. def load_intra_bg(path):
  1071. df = pd.read_csv(path, sep="\t")
  1072. df.columns = df.columns.str.strip()
  1073. df["rate_per_mb2"] = pd.to_numeric(df["rate_per_mb2"], errors="coerce")
  1074. return df.dropna(subset=["rate_per_mb2"])
  1077. def get_breakpoint_rates(trans_df, assembly):
  1078. rows = []
  1079. for label, sa, sb in BREAKPOINTS_STATS.get(assembly, []):
  1080. key_a, key_b = sorted([sa, sb])
  1081. match = trans_df[
  1082. (trans_df["scaffold_A"] == key_a) &
  1083. (trans_df["scaffold_B"] == key_b)]
  1084. if match.empty:
  1085. print(f" WARNING: {sa}/{sb} not found in trans table")
  1086. continue
  1087. rate = float(match["rate_per_mb2"].values[0])
  1088. rows.append({"label": label, "scaffold_A": sa, "scaffold_B": sb,
  1089. "rate_per_mb2": rate})
  1090. return pd.DataFrame(rows)
  1093. def empirical_pvalue(observed, background):
  1094. n = len(background)
  1095. count = np.sum(background >= observed)
  1096. return (count + 1) / (n + 1)
  1099. def wilcoxon_ranksum(bp_rates, background_rates):
  1100. stat, p = stats.mannwhitneyu(bp_rates, background_rates, alternative="greater")
  1101. return stat, p
  1104. def jitter(n, width=0.15, seed=42):
  1105. rng = np.random.default_rng(seed)
  1106. return rng.uniform(-width, width, n)
  1109. def draw_distribution(ax, rates, x_pos, colour, label, show_violin=True):
  1110. if len(rates) > 1 and show_violin:
  1111. parts = ax.violinplot([rates], positions=[x_pos],
  1112. widths=0.6, showmedians=True, showextrema=False)
  1113. for pc in parts["bodies"]:
  1114. pc.set_facecolor(colour)
  1115. pc.set_alpha(0.35)
  1116. pc.set_edgecolor(colour)
  1117. parts["cmedians"].set_color(colour)
  1118. parts["cmedians"].set_linewidth(1.5)
  1119. ax.scatter(np.full(len(rates), x_pos) + jitter(len(rates)),
  1120. rates, s=8, color=colour, alpha=0.5, linewidths=0, zorder=3)
  1121. ax.errorbar(x_pos, rates.mean(), yerr=rates.std(),
  1122. fmt="D", color=colour, markersize=4, capsize=3,
  1123. linewidth=1.0, zorder=5)
  1126. def draw_breakpoint_lines(ax, bp_df, colour, x_range):
  1127. y_span = ax.get_ylim()[1] - ax.get_ylim()[0]
  1128. min_y_gap = y_span * 0.02
  1129. nudge = y_span * 0.02
  1130. sorted_df = bp_df.sort_values("rate_per_mb2")
  1131. prev_y_label = -np.inf
  1132. for _, row in sorted_df.iterrows():
  1133. ax.axhline(row["rate_per_mb2"], color=colour, linewidth=0.5,
  1134. linestyle="--", alpha=0.8, zorder=4)
  1135. y_label = row["rate_per_mb2"]
  1136. if y_label - prev_y_label < min_y_gap:
  1137. y_label = prev_y_label + nudge
  1138. ax.text(x_range[1] - 0.05, y_label, f" {row['label']}",
  1139. va="center", ha="left", fontsize=8, color=colour)
  1140. prev_y_label = y_label
  1143. # %% [markdown]
  1144. # ### 2.4 Run analysis and generate plots
  1146. # %%
  1147. assemblies = [("mpibr", MPIBR_TRANS, MPIBR_INTRA),
  1148. ("dtol", DTOL_TRANS, DTOL_INTRA)]
  1150. n_asm = len(assemblies)
  1151. fig_w = 6 #16.0 / 2.54
  1152. fig_h = 6 # 8.0 / 2.54 * n_asm
  1153. fig, axes = plt.subplots(1, n_asm, figsize=(fig_w, fig_h), sharey=False)
  1154. if n_asm == 1:
  1155. axes = [axes]
  1157. all_results = []
  1158. dist_store = {}
  1160. for ax, (assembly, trans_path, intra_path) in zip(axes, assemblies):
  1161. print(f"\n{'='*60}\nAssembly: {assembly.upper()}\n{'='*60}")
  1162. trans_df = load_trans_bg(trans_path)
  1163. intra_df = load_intra_bg(intra_path)
  1164. print(f" Negative control pairs: {len(trans_df)}")
  1165. print(f" Positive control scaffolds: {len(intra_df)}")
  1167. bp_df = get_breakpoint_rates(trans_df, assembly)
  1168. if bp_df.empty:
  1169. print(" No breakpoint rates found — skipping.")
  1170. ax.set_visible(False)
  1171. continue
  1173. bp_keys = set()
  1174. for _, sa, sb in BREAKPOINTS_STATS.get(assembly, []):
  1175. bp_keys.add(tuple(sorted([sa, sb])))
  1176. bg_mask = ~trans_df.apply(
  1177. lambda r: tuple(sorted([r["scaffold_A"], r["scaffold_B"]])) in bp_keys, axis=1)
  1178. background = trans_df[bg_mask]["rate_per_mb2"].values
  1179. intra_rates = intra_df["rate_per_mb2"].values
  1180. bp_rates = bp_df["rate_per_mb2"].values
  1182. print(f"\n Background: mean={background.mean():.3f} "
  1183. f"median={np.median(background):.3f} sd={background.std():.3f} n={len(background)}")
  1184. print(f" Intra-scaffold: mean={intra_rates.mean():.3f} "
  1185. f"median={np.median(intra_rates):.3f} sd={intra_rates.std():.3f} n={len(intra_rates)}")
  1187. emp_pvals = []
  1188. print(f"\n Breakpoint rates (n={len(background)} background pairs):")
  1189. for _, row in bp_df.iterrows():
  1190. p_emp = empirical_pvalue(row["rate_per_mb2"], background)
  1191. emp_pvals.append(p_emp)
  1192. print(f" {row['label']} ({row['scaffold_A']} / {row['scaffold_B']}): "
  1193. f"rate={row['rate_per_mb2']:.3f} empirical p={p_emp:.4f}")
  1194. all_results.append({
  1195. "assembly": assembly, "label": row["label"],
  1196. "scaffold_A": row["scaffold_A"], "scaffold_B": row["scaffold_B"],
  1197. "rate_per_mb2": row["rate_per_mb2"],
  1198. "bg_mean": background.mean(), "bg_median": np.median(background),
  1199. "bg_sd": background.std(), "bg_n": len(background),
  1200. "empirical_p": p_emp, "wilcoxon_U": np.nan, "wilcoxon_p": np.nan,
  1201. })
  1203. colour_bp = COLOURS_STATS["breakpoint"][assembly]
  1204. colour_asm = "#2166AC" if assembly == "mpibr" else "#e65c00"
  1205. x_neg, x_pos = 1.0, 2.0
  1206. x_range = (0.3, 2.9)
  1208. draw_distribution(ax, background, x_neg, COLOURS_STATS["negative"], "Background trans")
  1209. draw_distribution(ax, intra_rates, x_pos, COLOURS_STATS["positive"], "Intra-scaffold")
  1210. draw_breakpoint_lines(ax, bp_df, colour_bp, x_range)
  1211. for _, row in bp_df.iterrows():
  1212. ax.scatter(x_neg, row["rate_per_mb2"], s=40, color=colour_bp,
  1213. zorder=6, marker="*", linewidths=0)
  1215. ax.set_xlim(*x_range)
  1216. ax.set_xticks([x_neg, x_pos])
  1217. ax.set_xticklabels(["Background\ntrans pairs", "Intra-scaffold\n(positive ctrl)"],
  1218. fontsize=8)
  1219. ax.set_ylabel("Contact rate (pairs / Mb²)", fontsize=9)
  1220. ax.set_title(assembly.upper(), fontsize=10, fontweight="bold", color=colour_asm)
  1221. ax.spines["top"].set_visible(False)
  1222. ax.spines["right"].set_visible(False)
  1224. _wilcoxon_p = np.nan
  1225. if len(bp_rates) > 1:
  1226. stat, _wilcoxon_p = wilcoxon_ranksum(bp_rates, background)
  1227. print(f"\n Wilcoxon (joint) n={len(bp_rates)} vs n={len(background)}: "
  1228. f"U={stat:.1f} p={_wilcoxon_p:.4f} (one-tailed, greater)")
  1229. all_results[-1]["wilcoxon_U"] = stat
  1230. all_results[-1]["wilcoxon_p"] = _wilcoxon_p
  1231. else:
  1232. print(f"\n Wilcoxon skipped for {assembly.upper()} (single breakpoint)")
  1234. dist_store[assembly] = {
  1235. "background": background, "bp_df": bp_df,
  1236. "emp_pvals": emp_pvals, "wilcoxon_p": _wilcoxon_p,
  1237. "colour_bp": colour_bp, "colour_asm": colour_asm,
  1238. }
  1240. handles_stats = [
  1241. mpatches.Patch(color=COLOURS_STATS["negative"], alpha=0.7, label="Background trans pairs"),
  1242. mpatches.Patch(color=COLOURS_STATS["positive"], alpha=0.7, label="Intra-scaffold long-range"),
  1243. plt.Line2D([0], [0], color=COLOURS_STATS["breakpoint"]["mpibr"],
  1244. linewidth=1.0, linestyle="--", label="MPIBR breakpoint rate"),
  1245. plt.Line2D([0], [0], color=COLOURS_STATS["breakpoint"]["dtol"],
  1246. linewidth=1.0, linestyle="--", label="DToL breakpoint rate"),
  1247. plt.Line2D([0], [0], marker="*", color="w",
  1248. markerfacecolor=COLOURS_STATS["breakpoint"]["mpibr"],
  1249. markersize=8, label="Breakpoint observed"),
  1250. ]
  1251. plt.tight_layout()
  1252. fig.subplots_adjust(bottom=0.18)
  1253. fig.legend(handles=handles_stats, loc="lower center", ncol=3, fontsize=9,
  1254. frameon=False, bbox_to_anchor=(0.5, 0.01))
  1255. plot_path = f"{OUTPUT_STATS}.pdf"
  1256. fig.savefig(plot_path, dpi=300, bbox_inches="tight")
  1257. plt.show()
  1258. print(f"\nSaved: {plot_path}")
  1261. # %% [markdown]
  1262. # ### 2.5 Background distribution histogram
  1264. # %%
  1265. n_dist = len(dist_store)
  1266. fig2, axes2 = plt.subplots(1, n_dist,
  1267. figsize=(6, 3),
  1268. sharey=False)
  1269. if n_dist == 1:
  1270. axes2 = [axes2]
  1272. for ax2, (assembly, store) in zip(axes2, dist_store.items()):
  1273. bg = store["background"]
  1274. bp_df = store["bp_df"]
  1275. emp_pvals = store["emp_pvals"]
  1276. colour_bp = store["colour_bp"]
  1277. colour_asm= store["colour_asm"]
  1279. counts, _, _ = ax2.hist(bg, bins=40, color=COLOURS_STATS["negative"],
  1280. alpha=0.75, edgecolor="white", linewidth=0.3,
  1281. zorder=1, label="Background trans pairs")
  1282. y_max = counts.max() * 1.15
  1283. ax2.set_ylim(0, y_max)
  1285. bp_sorted = sorted(zip(bp_df.iterrows(), emp_pvals),
  1286. key=lambda t: t[0][1]["rate_per_mb2"])
  1287. min_x_gap = (bg.max() - bg.min()) * 0.06
  1288. label_step = y_max * 0.18
  1289. prev_x, nudge_level = -np.inf, 0
  1291. for (_, row), p_emp in bp_sorted:
  1292. x_val = row["rate_per_mb2"]
  1293. ax2.axvline(x_val, color=colour_bp, linewidth=0.5, linestyle="--", zorder=3)
  1294. p_str = f"p = {p_emp:.4f}" if p_emp >= 0.0001 else "p < 0.0001"
  1295. if x_val - prev_x < min_x_gap:
  1296. nudge_level += 1
  1297. else:
  1298. nudge_level = 0
  1299. y_label = y_max * 0.97 - nudge_level * label_step
  1300. ax2.text(x_val, y_label, f" {row['label']}\n {p_str}",
  1301. va="top", ha="left", fontsize=8, color=colour_bp,
  1302. style="italic", zorder=4)
  1303. prev_x = x_val
  1305. wilcoxon_p = store["wilcoxon_p"]
  1306. if not np.isnan(wilcoxon_p):
  1307. p_str = f"p = {wilcoxon_p:.4f}" if wilcoxon_p >= 0.0001 else "p < 0.0001"
  1308. ax2.text(0.99, 0.15, f"Wilcoxon (joint)\n{p_str}",
  1309. transform=ax2.transAxes, va="bottom", ha="right",
  1310. fontsize=8, color="#FF6600",
  1311. bbox=dict(boxstyle="round,pad=0.2", fc="white",
  1312. ec="#FF6600", linewidth=0.6, alpha=0.6))
  1314. ax2.set_xlabel("Contact rate (pairs / Mb²)", fontsize=9)
  1315. ax2.set_ylabel("Number of scaffold pairs", fontsize=9)
  1316. ax2.set_title(assembly.upper(), fontsize=10, fontweight="bold", color=colour_asm)
  1317. ax2.spines["top"].set_visible(False)
  1318. ax2.spines["right"].set_visible(False)
  1320. plt.tight_layout()
  1321. dist_path = f"{OUTPUT_STATS}_distribution.pdf"
  1322. fig2.savefig(dist_path, dpi=300, bbox_inches="tight")
  1323. plt.show()
  1324. print(f"Saved: {dist_path}")
  1327. # %% [markdown]
  1328. # ### 2.6 Results table
  1330. # %%
  1331. if all_results:
  1332. results_df = pd.DataFrame(all_results)
  1333. tsv_path = f"{OUTPUT_STATS}_results.tsv"
  1334. results_df.to_csv(tsv_path, sep="\t", index=False, float_format="%.4f")
  1335. print(f"Saved: {tsv_path}\n")
  1336. display(results_df)
  1339. # %% [markdown]
  1340. # ---
  1341. # ## 3. Repeat content at junction windows
  1342. #
  1343. # Parses RepeatMasker `.out` files for junction and control scaffold windows and produces a stacked bar chart showing repeat class fractions per window.
  1345. # %% [markdown]
  1346. # ### 3.1 Parameters
  1348. # %%
  1349. REPEAT_OUTPUT_DIR = f"{BASE}/junction_repeats"
  1350. OUTPUT_REPEATS = "junction_repeat_content" # prefix; .tsv and .pdf produced
  1351. WINDOW_SIZE = 200_000
  1354. # %% [markdown]
  1355. # ### 3.2 Repeat class colours and label metadata
  1357. # %%
  1358. CLASS_COLOURS_RM = {
  1359. "SINE": "#EE6677",
  1360. "LINE": "#377EB8",
  1361. "DNA": "#4DAF4A",
  1362. "LTR": "#984EA3",
  1363. "Simple_repeat": "#FF7F00",
  1364. "Low_complexity": "#A65628",
  1365. "Unknown": "#999999",
  1366. "Other": "#CCCCCC",
  1367. }
  1369. LABEL_META_RM = {
  1370. "BP1_MPIBR_40_end": ("mpibr", False),
  1371. "BP1_MPIBR_44_start": ("mpibr", False),
  1372. "BP2_DToL_31_end": ("dtol", False),
  1373. "BP2_DToL_40_start": ("dtol", False),
  1374. "BP3_DToL_41_end": ("dtol", False),
  1375. "BP3_DToL_46_start": ("dtol", False),
  1376. "BP4_DToL_44_end": ("dtol", False),
  1377. "BP4_DToL_45_start": ("dtol", False),
  1378. "Ctrl_MPIBR_3_end": ("mpibr", True),
  1379. "Ctrl_MPIBR_3_start": ("mpibr", True),
  1380. "Ctrl_MPIBR_5_end": ("mpibr", True),
  1381. "Ctrl_MPIBR_5_start": ("mpibr", True),
  1382. "Ctrl_DToL_1_end": ("dtol", True),
  1383. "Ctrl_DToL_1_start": ("dtol", True),
  1384. "Ctrl_DToL_2_end": ("dtol", True),
  1385. "Ctrl_DToL_2_start": ("dtol", True),
  1386. }
  1387. LABEL_ASSEMBLY_RM = {k: v[0] for k, v in LABEL_META_RM.items()}
  1389. ASSEMBLY_COLOURS_RM = {"mpibr": "#2166AC", "dtol": "#e65c00"}
  1391. ROW_ORDER_RM = [
  1392. "BP1_MPIBR_40_end", "BP1_MPIBR_44_start",
  1393. "BP2_DToL_31_end", "BP2_DToL_40_start",
  1394. "BP3_DToL_41_end", "BP3_DToL_46_start",
  1395. "BP4_DToL_44_end", "BP4_DToL_45_start",
  1396. "Ctrl_MPIBR_3_end", "Ctrl_MPIBR_3_start",
  1397. "Ctrl_MPIBR_5_end", "Ctrl_MPIBR_5_start",
  1398. "Ctrl_DToL_1_end", "Ctrl_DToL_1_start",
  1399. "Ctrl_DToL_2_end", "Ctrl_DToL_2_start",
  1400. ]
  1403. # %% [markdown]
  1404. # ### 3.3 Parsing helpers
  1406. # %%
  1407. def normalise_class_rm(repeat_class):
  1408. rc = repeat_class.upper()
  1409. if rc.startswith("SINE"): return "SINE"
  1410. if rc.startswith("LINE"): return "LINE"
  1411. if rc.startswith("DNA"): return "DNA"
  1412. if rc.startswith("LTR"): return "LTR"
  1413. if "SIMPLE" in rc: return "Simple_repeat"
  1414. if "LOW_COMPLEX" in rc: return "Low_complexity"
  1415. if rc in ("UNKNOWN", "UNSPECIFIED"): return "Unknown"
  1416. return "Other"
  1419. def parse_repeatmasker_out(out_path):
  1420. rows = []
  1421. if not os.path.isfile(out_path):
  1422. print(f" WARNING: .out file not found: {out_path}")
  1423. return pd.DataFrame()
  1424. with open(out_path) as fh:
  1425. for _ in range(3):
  1426. fh.readline()
  1427. for line in fh:
  1428. line = line.rstrip("\n")
  1429. if not line.strip():
  1430. continue
  1431. parts = line.split()
  1432. if len(parts) < 11:
  1433. continue
  1434. try:
  1435. seq_name = parts[4]
  1436. start = int(parts[5])
  1437. end = int(parts[6])
  1438. rep_class = parts[10]
  1439. rep_length = end - start + 1
  1440. rows.append({
  1441. "seq_name": seq_name,
  1442. "start": start,
  1443. "end": end,
  1444. "rep_class": normalise_class_rm(rep_class),
  1445. "rep_length": rep_length,
  1446. })
  1447. except (ValueError, IndexError):
  1448. continue
  1449. return pd.DataFrame(rows)
  1452. def window_length_from_fasta(fa_path):
  1453. lengths = {}
  1454. current, length = None, 0
  1455. with open(fa_path) as fh:
  1456. for line in fh:
  1457. line = line.strip()
  1458. if line.startswith(">"):
  1459. if current:
  1460. lengths[current] = length
  1461. current = line[1:].split()[0]
  1462. length = 0
  1463. else:
  1464. length += len(line)
  1465. if current:
  1466. lengths[current] = length
  1467. return lengths
  1470. def summarise_repeats_rm(rm_df, seq_lengths):
  1471. rows = []
  1472. all_classes = list(CLASS_COLOURS_RM.keys())
  1473. for seq_name, seq_len in seq_lengths.items():
  1474. sub = rm_df[rm_df["seq_name"] == seq_name] if not rm_df.empty else pd.DataFrame()
  1475. row = {"label": seq_name, "seq_len": seq_len}
  1476. total_masked = 0
  1477. for cls in all_classes:
  1478. cls_bp = 0 if sub.empty else int(sub[sub["rep_class"] == cls]["rep_length"].sum())
  1479. row[f"{cls}_bp"] = cls_bp
  1480. row[f"{cls}_fraction"] = cls_bp / seq_len if seq_len > 0 else 0.0
  1481. total_masked += cls_bp
  1482. row["total_masked_bp"] = total_masked
  1483. row["total_repeat_fraction"] = total_masked / seq_len if seq_len > 0 else 0.0
  1484. rows.append(row)
  1485. return pd.DataFrame(rows)
  1488. # %% [markdown]
  1489. # ### 3.4 Parse RepeatMasker output
  1491. # %%
  1492. all_summaries = []
  1494. for asm in ["mpibr", "dtol"]:
  1495. rm_dir = os.path.join(REPEAT_OUTPUT_DIR, f"wd/{asm}")
  1496. fa_path = os.path.join(REPEAT_OUTPUT_DIR, f"fa/{asm}_junction_windows.fa")
  1497. out_file = os.path.join(rm_dir, f"{asm}_junction_windows.fa.out")
  1499. print(f"\n{asm.upper()}")
  1500. if not os.path.isfile(fa_path):
  1501. print(f" WARNING: window FASTA not found: {fa_path} — skipping")
  1502. continue
  1504. seq_lengths = window_length_from_fasta(fa_path)
  1505. rm_df = parse_repeatmasker_out(out_file)
  1506. summary = summarise_repeats_rm(rm_df, seq_lengths)
  1507. summary["assembly"] = asm
  1508. all_summaries.append(summary)
  1510. print(f" {'Label':<30} {'Total%':>7} {'SINE%':>6} {'LINE%':>6} "
  1511. f"{'DNA%':>6} {'LTR%':>6} {'Simple%':>8} {'Unknown%':>9}")
  1512. for _, row in summary.iterrows():
  1513. print(f" {row['label']:<30} {row['total_repeat_fraction']*100:>7.1f} "
  1514. f"{row['SINE_fraction']*100:>6.1f} {row['LINE_fraction']*100:>6.1f} "
  1515. f"{row['DNA_fraction']*100:>6.1f} {row['LTR_fraction']*100:>6.1f} "
  1516. f"{row['Simple_repeat_fraction']*100:>8.1f} "
  1517. f"{row['Unknown_fraction']*100:>9.1f}")
  1519. combined_rm = pd.concat(all_summaries, ignore_index=True)
  1520. combined_rm["_order"] = combined_rm["label"].map(
  1521. {lbl: i for i, lbl in enumerate(ROW_ORDER_RM)})
  1522. combined_rm = combined_rm.sort_values("_order").drop(columns="_order")
  1524. tsv_rm = f"{OUTPUT_REPEATS}.tsv"
  1525. frac_cols = ["label", "assembly", "seq_len", "total_repeat_fraction"] + [f"{c}_fraction" for c in CLASS_COLOURS_RM]
  1526. combined_rm[frac_cols].to_csv(tsv_rm, sep="\t", index=False, float_format="%.4f")
  1527. print(f"\nSummary TSV saved: {tsv_rm}")
  1530. # %% [markdown]
  1531. # ### 3.5 Plot repeat content
  1533. # %%
  1534. labels_rm = combined_rm["label"].tolist()
  1535. n_rm = len(labels_rm)
  1536. classes = list(CLASS_COLOURS_RM.keys())
  1538. fig_w = 6 #14.0 / 2.54
  1539. fig_h = 4.5 #max(6.0, n_rm * 0.55 + 2.0) / 2.54
  1540. fig3, ax3 = plt.subplots(figsize=(fig_w, fig_h))
  1542. y_positions = np.arange(n_rm)
  1543. bar_height = 0.55
  1544. lefts = np.zeros(n_rm)
  1546. for cls in classes:
  1547. fracs = combined_rm[f"{cls}_fraction"].values * 100
  1548. ax3.barh(y_positions, fracs, left=lefts,
  1549. height=bar_height, color=CLASS_COLOURS_RM[cls],
  1550. label=cls, linewidth=0)
  1551. lefts += fracs
  1553. for i, row in enumerate(combined_rm.itertuples()):
  1554. pct = row.total_repeat_fraction * 100
  1555. ax3.text(lefts[i] + 0.3, i, f"{pct:.1f}%", va="center", ha="left", fontsize=8)
  1557. ax3.set_yticks(y_positions)
  1558. ax3.set_yticklabels([])
  1559. ax3.tick_params(axis="y", length=0)
  1560. first_ctrl_idx = None
  1561. for i, lbl in enumerate(labels_rm):
  1562. meta_rm = LABEL_META_RM.get(lbl, ("mpibr", False))
  1563. asm, is_ctrl = meta_rm
  1564. colour = ASSEMBLY_COLOURS_RM[asm]
  1565. style = "italic" if is_ctrl else "normal"
  1566. prefix = "" if is_ctrl else ""
  1567. ax3.text(-0.5, i, prefix + lbl.replace("_", " "),
  1568. va="center", ha="right", fontsize=8,
  1569. color=colour, style=style)
  1570. if is_ctrl and first_ctrl_idx is None:
  1571. first_ctrl_idx = i
  1573. if first_ctrl_idx is not None:
  1574. ax3.axhline(first_ctrl_idx - 0.5, color="#444444",
  1575. linewidth=0.8, linestyle="--", zorder=5)
  1576. # ax3.text(max(lefts) * 0.5, first_ctrl_idx - 0.55,
  1577. # "breakpoints above | controls below",
  1578. # va="top", ha="center", fontsize=5,
  1579. # color="#444444", style="italic")
  1581. ax3.set_xlabel("Repeat content (%)", fontsize=9)
  1582. x_max = max(lefts) * 1.15
  1583. ax3.set_xlim(-0.5, x_max)
  1584. ax3.spines["top"].set_visible(False)
  1585. ax3.spines["right"].set_visible(False)
  1586. ax3.spines["left"].set_visible(False)
  1588. handles_rm = [mpatches.Patch(color=CLASS_COLOURS_RM[c], label=c) for c in classes]
  1589. plt.tight_layout()
  1590. fig3.subplots_adjust(bottom=0.20)
  1591. fig3.legend(handles=handles_rm, loc="lower center", fontsize=9,
  1592. frameon=False, ncol=4, bbox_to_anchor=(0.5, 0.01))
  1594. plot_rm = f"{OUTPUT_REPEATS}.pdf"
  1595. fig3.savefig(plot_rm, dpi=300, bbox_inches="tight")
  1596. plt.show()
  1597. print(f"Saved: {plot_rm}")

breakpoint_analysis.ipynb at commit efe218e, under MIT · at the source

Overview

Authors: Simone Daniela Rencken1,2, Georgi Tushev1, David Hain1,3, Elena Ciirdaeva1, Oleg Simakov4, Gilles Laurent1
ORCID iDs: Simone Daniela Rencken, Georgi Tushev, David Hain, Elena Ciirdaeva, Gilles Laurent
  1. Max Planck Institute for Brain Research Frankfurt am Main Germany
  2. Radboud University, Donders Institute for Brain, Cognition and Behaviour Nijmegen Netherlands
  3. Faculty of Biological Sciences, Goethe University Frankfurt am Main Germany
  4. Department of Neuroscience and Developmental Biology, University of Vienna Vienna Austria
Institutions: Max Planck Institute for Brain Research (Germany); Donders Institute for Brain, Cognition and Behaviour (Netherlands); Goethe University Frankfurt (Germany); University of Vienna (Austria)
Journal: eLife, volume 14, article RP107393
Dates: published online 28 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.107393 · PMID 42206967 · PMCID PMC13218726 · OpenAlex W4412827090
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.7554/elife.107393.2
Status: code verified
Categories: other (organism), cellular / molecular (subfield)
Methods: Statistics, Preprocessing, fMRI & imaging
Keywords: Sepia officinalis, cephalopod, genome assembly, Other
MeSH: Chromosomes*, Genome*, Sepia*, Animals (* major topic)
Journal subjects: Genetics and Genomics
Topic: Cephalopods and Marine Biology (Ecology, Evolution, Behavior and Systematics, Agricultural and Biological Sciences), according to OpenAlex
Funding: European Research Council (Advanced Grant CAMOUFLAGE 101141501, Horizon 2020 European Union Research and Innovation Programme grant No. 945026)
Citations: cited by 2 papers (Europe PMC); 195 references in the paper
Research resources: R v4.4.2 RRID:SCR_001905, STAR v2.7.11b RRID:SCR_004463, InterProScan v5.73–104 RRID:SCR_005829, bedtools v2.30 RRID:SCR_006646, featureCounts (Subread v2.0.8) RRID:SCR_012919, RepeatMasker v4.1.7-p1 RRID:SCR_012954, BUSCO v5.5.0 RRID:SCR_015008, RepeatModeler v2.0.6 RRID:SCR_015027, DESeq2 v1.42.0 RRID:SCR_015687, DIAMOND2 RRID:SCR_016071, StringTie v3.0.0 RRID:SCR_016323, VecScreen RRID:SCR_016577, clusterProfiler v4.12.6 RRID:SCR_016884, GenomeScope 2.0 RRID:SCR_017014, OrthoFinder v2.5 RRID:SCR_017118, ape v5.8.1 (R package) RRID:SCR_017343, TransDecoder v5.7.0 RRID:SCR_017647, minimap2 RRID:SCR_018550, CAFE5 v5.1.1 RRID:SCR_018924, seqtk RRID:SCR_018927, mosdepth RRID:SCR_018929, BRAKER3 (incl. TSEBRA) RRID:SCR_018964, pysam v0.22.1 RRID:SCR_021017, hifiasm RRID:SCR_021069, eggNOG-mapper v2.1.12 RRID:SCR_021165, MCScanX RRID:SCR_022067, bwa-mem2 v2.3 RRID:SCR_022192, YAHS RRID:SCR_022965, pairtools v1.1.0 RRID:SCR_023038, Winnowmap2 RRID:SCR_025349, apeglm RRID:SCR_026951

Abstract

Coleoid cephalopods, a subclass of mollusks that includes octopuses, cuttlefish, and squid, exhibit sophisticated biological features, such as dynamic and neurally driven camouflage behavior, inter-individual communication, single-lens camera-like eyes, the largest brains among invertebrates, and a distinctive embryonic development. The common cuttlefish Sepia officinalis has served as a model organism in various research fields, spanning biophysics, neurobiology, behavior, evolution, ecology, and biomechanics. More recently, it has become a model to investigate the neural mechanisms underlying cephalopod camouflage, using quantitative behavioral approaches alongside molecular techniques to characterize the identity, evolution, and development of neuronal cell types. Despite significant interest in this animal, a high-quality, annotated genome of this species is still lacking. To address this, we sequenced and assembled a chromosome-scale genome for S. officinalis. Our assembly spans 5.68 billion base pairs and comprises 1n=47 repeat-rich chromosome scaffolds. This was unexpected because the haploid karyotypes of other decapods indicate 46 chromosomes. Detailed comparisons of our data to those from published decapod genome assemblies and to another recent genome assembly of S. officinalis (itself suggesting 1n=49 chromosomes) in fact revealed clear homologies between 46 scaffolds across all the datasets. In-depth comparison of datasets reveals highly repetitive regions at discordant scaffold boundaries and suggests that the true karyotype of S. officinalis is probably 1n=46 chromosomes, a likely ancestral and if true, conserved decapod karyotype. Our results include a comprehensive gene annotation and full-length transcript prediction, which we used to characterize orthologous gene families across mollusks. We identified several large-scale expansions specific to cephalopods, with many genes specific to neural or non-neural tissues of adult S. officinalis. In summary, this genome should provide a valuable resource for future research on the evolution, brain organization, information processing, development, and behavior in this important clade.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 40 matches between paragraphs and lines of code.

gitlab.mpcdf.mpg.de/mpibr/laur/cuttlefishomics/soffgenome

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: efe218e0bf3c8318c5c729cbab94478dbf5b5502, 28 April 2026
Languages: Shell (50), R (8), Perl (5), Python (2), Jupyter (1)
Size: 80 files, 66 scripts
Software Heritage: archived
Found in: the references
Holds: README, license file, environment (envs/environment.yml), 3 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: SAMtools (15 files), tidyverse (8 files), STAR (6 files), circlize (3 files), ComplexHeatmap (3 files), BEDTools (2 files), clusterProfiler (2 files), cowplot (2 files), DESeq2 (2 files), ggplot2 (2 files), pandas (2 files), pheatmap (2 files), pysam (2 files), Subread (featureCounts) (2 files), data.table (1 file), Matplotlib (1 file), NumPy (1 file), SciPy (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
68 files
At the source: archive.softwareheritage.org/swh:1:dir:55d0af83f088bd734c66ca681fe48855e35eae83;origin=https://gitlab.mpcdf.mpg.de/mpibr/laur/cuttlefishomics/soffgenome;visit=swh:1:snp:ef236282e660b36085bf722edd6916303f50ce33;anchor=swh:1:rev:efe218e0bf3c8318c5c729cbab94478dbf5b5502

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 66 scripts, each with its path and the digest of its content;
  • 40 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

The genome assembly and raw data can be found at the BioProject PRJNA1091451 on NCBI. Raw sequencing reads are deposited at SRA (study accession SRP570862). The code for the genome assembly and annotation is available at https://gitlab.mpcdf.mpg.de/mpibr/laur/cuttlefishomics/soffgenome (copy archived at Tushev, 2026). Genome annotation files are deposited at https://doi.org/10.17617/3.CGO7QG.

The following datasets were generated:

RenckenS TushevG HainD CiirdaevaE SimakovO LaurentG 2025Sepia officinalis isolate:GLC-03058 (common cuttlefish)NCBI BioProjectPRJNA109145110.7554/eLife.107393PMC1321872642206967

RenckenS TushevG HainD CiirdaevaE SimakovO LaurentG 2025Chromosome-scale genome assembly of the European common cuttlefish Sepia officinalisEdmond - the Open Research Data Repository of the Max Planck Society10.17617/3.CGO7QGPMC1321872642206967

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 6 authors, 4 keywords, 4 MeSH terms, 1 funder, 183 references, 31 RRIDs.

Cite

This paper

Rencken, S. D., Tushev, G., Hain, D., Ciirdaeva, E., Simakov, O., & Laurent, G. (2026). Chromosome-scale genome assembly of the European common cuttlefish &lt;i&gt;Sepia officinalis&lt;/i&gt;. eLife, 14, RP107393. https://doi.org/10.7554/elife.107393

BibTeX

@article{rencken2026chromosome,
author = {Rencken, Simone Daniela and Tushev, Georgi and Hain, David and Ciirdaeva, Elena and Simakov, Oleg and Laurent, Gilles},
title = {{Chromosome-scale genome assembly of the European common cuttlefish \&lt;i\&gt;Sepia officinalis\&lt;/i\&gt;}},
journal = {eLife},
year = {2026},
month = may,
volume = {14},
pages = {RP107393},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.107393},
url = {https://doi.org/10.7554/elife.107393},
pmid = {42206967},
pmcid = {PMC13218726}
}

RIS

TY - JOUR
AU - Rencken, Simone Daniela
AU - Tushev, Georgi
AU - Hain, David
AU - Ciirdaeva, Elena
AU - Simakov, Oleg
AU - Laurent, Gilles
TI - Chromosome-scale genome assembly of the European common cuttlefish &lt;i&gt;Sepia officinalis&lt;/i&gt;
T2 - eLife
J2 - Elife
PY - 2026
DA - 2026/05/28
VL - 14
SP - RP107393
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.107393
UR - https://doi.org/10.7554/elife.107393
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.107393",
"type": "article-journal",
"title": "Chromosome-scale genome assembly of the European common cuttlefish &lt;i&gt;Sepia officinalis&lt;/i&gt;",
"container-title": "eLife",
"author": [
{
"family": "Rencken",
"given": "Simone Daniela"
},
{
"family": "Tushev",
"given": "Georgi"
},
{
"family": "Hain",
"given": "David"
},
{
"family": "Ciirdaeva",
"given": "Elena"
},
{
"family": "Simakov",
"given": "Oleg"
},
{
"family": "Laurent",
"given": "Gilles"
}
],
"container-title-short": "Elife",
"volume": "14",
"page": "RP107393",
"DOI": "10.7554/elife.107393",
"PMID": "42206967",
"PMCID": "PMC13218726",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.107393",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
28
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41586-026-10512-9 [code]
Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.
Journal: Nature
In common: Subread (featureCounts), STAR, pysam, 15 other tools, cellular / molecular, 4 references
[2] doi:10.21203/rs.3.rs-9927928/v1 [code]
Genome-wide and allele-resolved maps of the radial architecture of the mouse genome
Journal: Research Square (preprint)
In common: Subread (featureCounts), STAR, pysam, 13 other tools, 4 references
[3] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: Subread (featureCounts), STAR, SAMtools, 13 other tools, 4 references
[4] doi:10.1038/s41467-026-69944-6 [code]
Multi-modal dissection of cell-type specific TDP-43 pathology in the motor cortex.
Journal: Nature communications
In common: pysam, BEDTools, SAMtools, 12 other tools, 5 references
[5] doi:10.1038/s42003-026-10402-w [code]
Temporal orchestration of transcriptional and epigenomic programming underlying maternal embryonic diapause in a cricket model.
Journal: Communications biology
In common: SAMtools, DESeq2, pheatmap, 2 other tools, other, 12 references
[6] doi:10.1126/sciadv.aed2952 [code]
Activation of transposable elements is linked to a region- and cell type-specific interferon response in Parkinson's disease.
Journal: Science advances
In common: Subread (featureCounts), STAR, pysam, 12 other tools, cellular / molecular, 3 references
[7] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: Subread (featureCounts), pysam, BEDTools, 12 other tools, cellular / molecular, 2 references
[8] doi:10.1038/s41586-026-10629-x [code]
Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: BEDTools, circlize, DESeq2, 9 other tools, other, cellular / molecular, 5 references
[9] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: STAR, pysam, BEDTools, 11 other tools, 2 references
[10] doi:10.1038/s41467-026-71790-5 [code]
Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification.
Journal: Nature communications
In common: pysam, BEDTools, SAMtools, 12 other tools, cellular / molecular

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.