diff --git a/src/cell_metrics.py b/src/cell_metrics.py index eeba999..7aa1d69 100644 --- a/src/cell_metrics.py +++ b/src/cell_metrics.py @@ -45,6 +45,30 @@ def safe_prop(numer, denom): arr = np.where(denom > 0, (numer / denom) * 100.0, 0.0) return pd.Series(arr, index=numer.index) + # ---- RT-switching junction column naming (shared by both mode paths) ---- + # The four splice-junction types, and the per-cell columns that feed the + # report's two RT-switching violins. Totals already exist as + # `{Known,Novel}_{canonical,non_canonical}_junctions`; here we add the + # RTS-flagged numerators (all-junction plot) plus the distinct-intron totals + # and RTS numerators (unique-junction plot). Column presence is the report's + # capability signal: emitted only when the junctions file carries the needed + # inputs (RTS_junction for the all-junction plot; RTS_junction + genomic + # coords for the unique-junction plot), so the report skips a plot exactly + # when its columns are absent -- mirroring the old runtime guards. + junc_type_order = ['known_canonical', 'known_non_canonical', 'novel_canonical', 'novel_non_canonical'] + junc_total_name = { + 'known_canonical': 'Known_canonical_junctions', + 'known_non_canonical': 'Known_non_canonical_junctions', + 'novel_canonical': 'Novel_canonical_junctions', + 'novel_non_canonical': 'Novel_non_canonical_junctions', + } + rts_junc_name = {jc: 'RTS_' + junc_total_name[jc] for jc in junc_type_order} + unique_junc_name = {jc: 'unique_' + junc_total_name[jc] for jc in junc_type_order} + rts_unique_junc_name = {jc: 'RTS_unique_' + junc_total_name[jc] for jc in junc_type_order} + + def _rts_true(series): + return series.astype(str).str.lower().isin(['true', 't', '1', 'yes']) + warnings.simplefilter('ignore', PerformanceWarning) def _isoforms_summary(cls, junc): @@ -244,6 +268,82 @@ def weighted_median(): for jc in junc_types: base = per_tx[jc].fillna(0).to_numpy() if jc in per_tx.columns else np.zeros(T) summary[f"{tag}_{jc}_junctions"] = wsum(base * cmask) + + # ---- RT-switching numerators (all-junction violin) ---- + # RTS-flagged per-transcript junction-type counts, FL-weighted like the + # totals above. Emitted only if the junctions file has RTS_junction; the + # denominators are the {Known,Novel}_*_junctions totals. The report + # derives % = 100 * RTS / total per SJ type per cell. + if 'RTS_junction' in junc.columns: + jt_rts = junc.loc[_rts_true(junc['RTS_junction']), + ['isoform', 'junction_category', 'canonical']].copy() + jt_rts['jtype'] = jt_rts['junction_category'].astype(str) + '_' + jt_rts['canonical'].astype(str) + per_tx_rts = (jt_rts.groupby(['isoform', 'jtype']).size().unstack(fill_value=0) + .reindex(cls['isoform'].to_numpy())) + for jc in junc_types: + w_T = per_tx_rts[jc].fillna(0).to_numpy() if jc in per_tx_rts.columns else np.zeros(T) + summary[rts_junc_name[jc]] = wsum(w_T) + + # ---- unique (distinct-intron) junction counts per cell ---- + # A junction is "present" in a cell iff ANY transcript carrying it is + # expressed there. Collapsing the SAME genomic intron across transcripts + # makes this a set-union, not an FL-weighted sum, so it needs the + # junction x transcript incidence times the transcript x cell presence: + # P = Jt @ Xp, binarized. Distinct-per-cell counts (and their RTS subset) + # are then column sums of P over each SJ type's rows. Kept fully sparse so + # the (junction, cell) product is never densified. Emitted only when the + # junctions file has RTS_junction + genomic coordinates. + _chrom_col = 'chrom' if 'chrom' in junc.columns else ('chr' if 'chr' in junc.columns else None) + have_coords = ('RTS_junction' in junc.columns and _chrom_col is not None + and {'genomic_start_coord', 'genomic_end_coord'}.issubset(junc.columns)) + if have_coords: + from scipy import sparse + iso_to_idx = pd.Series(np.arange(T), index=cls['isoform'].astype(str).to_numpy()) + jc_df = junc[['isoform', 'junction_category', 'canonical', _chrom_col, + 'genomic_start_coord', 'genomic_end_coord', 'RTS_junction']].copy() + if 'strand' in junc.columns: + jc_df['strand'] = junc['strand'].astype(str) + else: + jc_df['strand'] = '' + t_idx = iso_to_idx.reindex(jc_df['isoform'].astype(str).to_numpy()).to_numpy() + gs = pd.to_numeric(jc_df['genomic_start_coord'], errors='coerce').to_numpy() + ge = pd.to_numeric(jc_df['genomic_end_coord'], errors='coerce').to_numpy() + chrom_s = jc_df[_chrom_col].astype(str).to_numpy() + keep = (~np.isnan(t_idx)) & (~np.isnan(gs)) & (~np.isnan(ge)) & (chrom_s != '') & (chrom_s != 'nan') + if keep.any(): + t_idx = t_idx[keep].astype(np.int64) + jkey = (chrom_s[keep] + ':' + gs[keep].astype(np.int64).astype(str) + ':' + + ge[keep].astype(np.int64).astype(str) + ':' + jc_df['strand'].to_numpy()[keep]) + jkey_code, _uniq = pd.factorize(jkey) + nJ = len(_uniq) + jtype_arr = (jc_df['junction_category'].astype(str) + '_' + + jc_df['canonical'].astype(str)).to_numpy()[keep] + rts_arr = _rts_true(jc_df['RTS_junction']).to_numpy()[keep] + # per distinct junction: its SJ type (consistent across rows) and + # whether ANY supporting row is RTS. + jtype_J = np.empty(nJ, dtype=object) + jtype_J[jkey_code] = jtype_arr + rts_J = np.zeros(nJ, dtype=bool) + np.logical_or.at(rts_J, jkey_code, rts_arr) + # Jt: distinct junction x transcript incidence (binary) + Jt = sparse.coo_matrix((np.ones(len(jkey_code)), (jkey_code, t_idx)), + shape=(nJ, T)).tocsr() + Jt.data[:] = 1.0 + # Xp: transcript x cell presence (binary) from the parsed COO + Xp = sparse.coo_matrix((np.ones(len(row)), (row, col)), shape=(T, Cn)).tocsr() + Xp.data[:] = 1.0 + P = Jt.dot(Xp) # (nJ x Cn): supporting transcripts per (junction, cell) + P.data[:] = 1.0 # binarize -> present / absent + for jc in junc_types: + rmask = (jtype_J == jc) + if rmask.any(): + u_tot = np.asarray(P[rmask, :].sum(axis=0)).ravel() + u_rts = np.asarray(P[rmask & rts_J, :].sum(axis=0)).ravel() + else: + u_tot = np.zeros(Cn) + u_rts = np.zeros(Cn) + summary[unique_junc_name[jc]] = pd.Series(u_tot, index=cells_index) + summary[rts_unique_junc_name[jc]] = pd.Series(u_rts, index=cells_index) else: for cn in list(junc_rename.values()) + ['total_junctions'] + [f"{v}_prop" for v in junc_rename.values()]: summary[cn] = 0 @@ -498,8 +598,14 @@ def lenbin_masks(base_mask): # the classification table. That CB column is the bulk of the per-cell-enriched # junctions file (~70GB for Isosceles, dense EM output), so excluding it from the # read is the difference between OOM and a few GB. Reads mode still needs it. + # RTS_junction + genomic coordinates feed the two RT-switching-by-SJ-type + # violins. They are small next to the CB column that dominates the isoforms + # junctions file (~70GB for dense EM output), so adding them keeps the read + # cheap while letting the metrics be computed here instead of re-exploding + # junctions x cells in the report. _junc_usecols = {'isoform','readID','read_id','ID','read_name','read', - 'junction_category','canonical'} + 'junction_category','canonical','RTS_junction', + 'genomic_start_coord','genomic_end_coord','chrom','chr','strand'} if args.mode != 'isoforms': _junc_usecols.add('CB') try: @@ -720,6 +826,52 @@ def lenbin_masks(base_mask): if jc in wide.columns: summary[f"{tag}_{jc}_junctions"] = wide[jc].reindex(summary.index, fill_value=0) + # ---- RT-switching all-junction / unique-junction counts per cell ---- + # Reads mode: each junction row already carries its read's CB (merged from + # the classification table when the junctions file lacks CB). The + # all-junction numerators count RTS rows per (CB, SJ type); the unique-junction + # metric collapses identical genomic introns per cell first (a junction is + # RTS if ANY supporting row is). Columns emitted only when RTS_junction (and, + # for the unique plot, genomic coordinates) are present, so the report skips a + # plot exactly when its columns are absent. + if (not junc.empty and {'junction_category', 'canonical', 'RTS_junction'}.issubset(junc.columns)): + jrt = junc.copy() + iso_col_j = next((c for c in ['isoform', 'readID', 'read_id', 'ID', 'read_name', 'read'] + if c in jrt.columns and c in cls_valid.columns), None) + if ('CB' not in jrt.columns or (jrt['CB'].fillna('') == '').all()) and iso_col_j is not None: + jrt = jrt.merge(cls_valid[[iso_col_j, 'CB']].dropna().drop_duplicates(iso_col_j), + on=iso_col_j, how='left') + if 'CB' in jrt.columns: + jrt = jrt[(jrt['CB'].notna()) & (jrt['CB'] != '')].copy() + if not jrt.empty: + jrt['jtype'] = jrt['junction_category'].astype(str) + '_' + jrt['canonical'].astype(str) + jrt['rts'] = _rts_true(jrt['RTS_junction']) + # all-junction RTS numerators (denominators = *_junctions totals) + rts_all = (jrt[jrt['rts']].groupby(['CB', 'jtype']).size() + .unstack(fill_value=0)) + for jc in junc_type_order: + summary[rts_junc_name[jc]] = (rts_all[jc].reindex(summary.index, fill_value=0) + if jc in rts_all.columns else 0) + # unique-junction counts: distinct genomic introns per (CB, SJ type) + _chrom_col = 'chrom' if 'chrom' in jrt.columns else ('chr' if 'chr' in jrt.columns else None) + if _chrom_col is not None and {'genomic_start_coord', 'genomic_end_coord'}.issubset(jrt.columns): + strand_s = jrt['strand'].astype(str) if 'strand' in jrt.columns else '' + jrt['jkey'] = (jrt[_chrom_col].astype(str) + ':' + + jrt['genomic_start_coord'].astype(str) + ':' + + jrt['genomic_end_coord'].astype(str) + ':' + strand_s) + valid = jrt['genomic_start_coord'].notna() & jrt['genomic_end_coord'].notna() + ju = jrt[valid] + if not ju.empty: + dd = (ju.groupby(['CB', 'jkey']) + .agg(jtype=('jtype', 'first'), rts=('rts', 'any')).reset_index()) + uniq_tot = dd.groupby(['CB', 'jtype']).size().unstack(fill_value=0) + uniq_rts = dd[dd['rts']].groupby(['CB', 'jtype']).size().unstack(fill_value=0) + for jc in junc_type_order: + summary[unique_junc_name[jc]] = (uniq_tot[jc].reindex(summary.index, fill_value=0) + if jc in uniq_tot.columns else 0) + summary[rts_unique_junc_name[jc]] = (uniq_rts[jc].reindex(summary.index, fill_value=0) + if jc in uniq_rts.columns else 0) + sublevels = { 'full-splice_match': ['alternative_3end','alternative_3end5end','alternative_5end','reference_match','mono-exon'], 'incomplete-splice_match': ['3prime_fragment','internal_fragment','5prime_fragment','intron_retention','mono-exon'], diff --git a/src/report_assets/SQANTI-sc_report.R b/src/report_assets/SQANTI-sc_report.R index d7d434a..9ad0400 100644 --- a/src/report_assets/SQANTI-sc_report.R +++ b/src/report_assets/SQANTI-sc_report.R @@ -2538,17 +2538,6 @@ generate_sqantisc_plots <- function(SQANTI_cell_summary, Classification_file, Ju # t1 <- ttheme_default(core=list(core = list(fg_params = list(cex = 0.6)), # colhead = list(fg_params = list(cex = 0.7)))) - # Resolve common isoform ID key between Junctions and Classification_file once, - # shared by the RT-switching junction blocks below (junc_rt). The per-category - # junction-type violins no longer need it — they read the cell summary directly. - junc_iso_key <- NULL - for (.k in c("isoform", "readID", "read_id", "ID", "read_name", "read")) { - if (.k %in% colnames(Junctions) && .k %in% colnames(Classification_file)) { - junc_iso_key <- .k - break - } - } - # Build SJ per-type/per-category plots and the all-canonical grouped plot for HTML (and reuse for PDF) # This block creates plot objects regardless of generate_pdf so the Rmd can render them. { @@ -2720,94 +2709,52 @@ generate_sqantisc_plots <- function(SQANTI_cell_summary, Classification_file, Ju ) } - # NEW: RT-switching by splice junction type across cells (all and unique junctions) - if ("RTS_junction" %in% colnames(Junctions)) { - # Normalize boolean - rts_bool_vec <- tolower(as.character(Junctions$RTS_junction)) %in% c("true", "t", "1", "yes") - - # Optional genomic coordinates for collapsing duplicate junctions (same chr + positions) - rts_coord_cols <- intersect( - c("chrom", "chr", "genomic_start_coord", "genomic_end_coord", "strand"), - colnames(Junctions) - ) - has_rts_junc_coords <- all(c("genomic_start_coord", "genomic_end_coord") %in% rts_coord_cols) && - ("chrom" %in% rts_coord_cols || "chr" %in% rts_coord_cols) - - # In isoforms mode, Junctions$CB is a comma-separated list — must explode via classification - if (mode == "isoforms" && !is.null(junc_iso_key) && - "FL" %in% colnames(Classification_file) && "CB" %in% colnames(Classification_file) && - junc_iso_key %in% colnames(Junctions)) { - - # Explode Classification_file CB/FL to per-(isoform, cell, FL_count) - cls_exp_rts <- Classification_file %>% - filter(!is.na(CB)) %>% - select(all_of(c(junc_iso_key, "CB", "FL"))) %>% - mutate(CB_raw = as.character(CB), FL_raw = as.character(FL)) %>% - tidyr::separate_rows(CB_raw, FL_raw, sep = ",") %>% - mutate(CB_clean = trimws(CB_raw), - FL_num = suppressWarnings(as.numeric(trimws(FL_raw)))) %>% - filter(CB_clean != "", CB_clean != "unassigned", !is.na(FL_num), FL_num > 0) %>% - select(all_of(c(junc_iso_key, "CB_clean", "FL_num"))) - - junc_rt_cols <- unique(c( - junc_iso_key, "junction_category", "canonical", "RTS_junction", - rts_coord_cols - )) - junc_rt_cols <- junc_rt_cols[junc_rt_cols %in% colnames(Junctions)] - - # Join junctions (isoform, junction type, RTS flag) to exploded classification - junc_rt <- Junctions %>% - select(all_of(junc_rt_cols)) %>% - mutate(SJ_type = paste(junction_category, canonical, sep = "_"), - RTS_bool = tolower(as.character(RTS_junction)) %in% c("true", "t", "1", "yes")) %>% - inner_join(cls_exp_rts, by = junc_iso_key) %>% - rename(CB = CB_clean, count = FL_num) - - } else { - # Reads mode: each junction already has one CB, count = 1 - junc_rt_cols <- unique(c( - "CB", "junction_category", "canonical", "RTS_junction", - rts_coord_cols - )) - junc_rt_cols <- junc_rt_cols[junc_rt_cols %in% colnames(Junctions)] - - junc_rt <- Junctions %>% - select(all_of(junc_rt_cols)) %>% - dplyr::filter(!is.na(CB) & CB != "" & CB != "unassigned") %>% - mutate(SJ_type = paste(junction_category, canonical, sep = "_"), - RTS_bool = rts_bool_vec, - count = 1) - } - - # Ensure consistent SJ type levels and labels - sj_levels <- c("known_canonical", "known_non_canonical", "novel_canonical", "novel_non_canonical") - sj_labels <- c("Known\nCanonical", "Known\nNon-canonical", "Novel\nCanonical", "Novel\nNon-canonical") - junc_rt$SJ_type <- factor(junc_rt$SJ_type, levels = sj_levels) - - # Color map consistent with other SJ type plots - sj_fill_map <- c( - known_canonical = "#6BAED6", - known_non_canonical = "goldenrod1", - novel_canonical = "#78C679", - novel_non_canonical = "#FC8D59" - ) - - # Per-cell percentages for ALL junctions (FL-weighted in isoforms mode) - all_junc_by_cell <- junc_rt %>% - group_by(CB, SJ_type) %>% - summarise(total = sum(count, na.rm = TRUE), - rts = sum(count[RTS_bool], na.rm = TRUE), .groups = "drop") %>% - # Only include (CB, SJ_type) that actually have data; NA for absent pairs - mutate(perc = ifelse(total > 0, 100 * rts / total, NA_real_)) %>% - tidyr::complete(CB, SJ_type = sj_levels) + # RT-switching by splice junction type across cells (all and unique junctions), + # sourced from cell-summary columns (see cell_metrics.py) instead of re-exploding + # junctions x cells here. The all-junction plot is 100 * RTS__junctions + # over the _junctions totals; the unique plot is + # 100 * RTS_unique__junctions over unique__junctions. Percentage is + # NA when a cell has no junctions of a type, so that (cell, type) drops out of the + # violin. The metrics module emits these columns only when the junctions file + # carried the needed inputs (RTS flag; RTS flag + genomic coords for the unique + # plot), so column presence is the signal to draw each plot. + sj_levels <- c("known_canonical", "known_non_canonical", "novel_canonical", "novel_non_canonical") + sj_labels <- c("Known\nCanonical", "Known\nNon-canonical", "Novel\nCanonical", "Novel\nNon-canonical") + sj_fill_map <- c( + known_canonical = "#6BAED6", + known_non_canonical = "goldenrod1", + novel_canonical = "#78C679", + novel_non_canonical = "#FC8D59" + ) + sj_colbase <- c( + known_canonical = "Known_canonical_junctions", + known_non_canonical = "Known_non_canonical_junctions", + novel_canonical = "Novel_canonical_junctions", + novel_non_canonical = "Novel_non_canonical_junctions" + ) - df_long_all <- data.frame( - Variable = factor(all_junc_by_cell$SJ_type, levels = sj_levels), - Value = all_junc_by_cell$perc - ) + build_rts_sjtype_long <- function(summ, num_cols, den_cols) { + n <- nrow(summ) + getcol <- function(cn) if (cn %in% names(summ)) suppressWarnings(as.numeric(summ[[cn]])) else rep(NA_real_, n) + parts <- lapply(sj_levels, function(sj) { + num <- getcol(num_cols[[sj]]) + den <- getcol(den_cols[[sj]]) + data.frame( + Variable = sj, + Value = ifelse(!is.na(den) & den > 0, 100 * num / den, NA_real_), + stringsAsFactors = FALSE + ) + }) + out <- do.call(rbind, parts) + out$Variable <- factor(out$Variable, levels = sj_levels) + out + } + # All junctions: numerators RTS__junctions over the existing totals. + num_all <- setNames(paste0("RTS_", sj_colbase), names(sj_colbase)) + if (all(num_all %in% colnames(SQANTI_cell_summary))) { gg_rts_all_by_sjtype <<- build_violin_plot( - df_long = df_long_all, + df_long = build_rts_sjtype_long(SQANTI_cell_summary, num_all, sj_colbase), title = "RT-switching All Junctions by Splice Junction Type Across Cells", x_labels = sj_labels, fill_map = sj_fill_map, @@ -2821,69 +2768,31 @@ generate_sqantisc_plots <- function(SQANTI_cell_summary, Classification_file, Ju violin_outline_fill = TRUE, box_outline_default = "grey20" ) + } else { + message("RT-switching All Junctions plot skipped: RTS junction columns not in cell summary.") + } - # Unique junctions per cell: collapse identical genomic introns (reads mode: 1 per junction; - # isoforms mode: 1 per junction per cell, FL not used). RTS if any supporting row is RTS. - if (has_rts_junc_coords) { - uniq_junc_by_cell <- junc_rt %>% - mutate( - junc_chr = dplyr::coalesce( - if ("chrom" %in% names(.)) as.character(.data$chrom) else NA_character_, - if ("chr" %in% names(.)) as.character(.data$chr) else NA_character_ - ), - strand_part = if ("strand" %in% names(.)) as.character(.data$strand) else "", - junc_key = paste( - junc_chr, - as.character(.data$genomic_start_coord), - as.character(.data$genomic_end_coord), - strand_part, - sep = ":" - ) - ) %>% - dplyr::filter( - !is.na(junc_chr), nzchar(junc_chr), - !is.na(.data$genomic_start_coord), !is.na(.data$genomic_end_coord) - ) %>% - group_by(CB, SJ_type, junc_key) %>% - summarise(RTS_any = any(RTS_bool, na.rm = TRUE), .groups = "drop") %>% - group_by(CB, SJ_type) %>% - summarise( - total = dplyr::n(), - rts = sum(RTS_any, na.rm = TRUE), - .groups = "drop" - ) %>% - mutate(perc = ifelse(total > 0, 100 * rts / total, NA_real_)) %>% - tidyr::complete(CB, SJ_type = sj_levels) - - df_long_unique <- data.frame( - Variable = factor(uniq_junc_by_cell$SJ_type, levels = sj_levels), - Value = uniq_junc_by_cell$perc - ) - - gg_rts_unique_by_sjtype <<- build_violin_plot( - df_long = df_long_unique, - title = "RT-switching Unique Junctions by Splice Junction Type Across Cells", - x_labels = sj_labels, - fill_map = sj_fill_map, - y_label = "Junctions, %", - legend = FALSE, - ylim = c(0, 100), - violin_alpha = 0.7, - box_alpha = 0.3, - box_width = 0.05, - x_tickangle = 45, - violin_outline_fill = TRUE, - box_outline_default = "grey20" - ) - } else { - message( - "RT-switching Unique Junctions plot skipped: need chrom (or chr), genomic_start_coord, ", - "and genomic_end_coord in junctions file." - ) - } - + # Unique junctions: distinct-intron numerators/denominators from the cell summary. + num_uniq <- setNames(paste0("RTS_unique_", sj_colbase), names(sj_colbase)) + den_uniq <- setNames(paste0("unique_", sj_colbase), names(sj_colbase)) + if (all(den_uniq %in% colnames(SQANTI_cell_summary))) { + gg_rts_unique_by_sjtype <<- build_violin_plot( + df_long = build_rts_sjtype_long(SQANTI_cell_summary, num_uniq, den_uniq), + title = "RT-switching Unique Junctions by Splice Junction Type Across Cells", + x_labels = sj_labels, + fill_map = sj_fill_map, + y_label = "Junctions, %", + legend = FALSE, + ylim = c(0, 100), + violin_alpha = 0.7, + box_alpha = 0.3, + box_width = 0.05, + x_tickangle = 45, + violin_outline_fill = TRUE, + box_outline_default = "grey20" + ) } else { - message("RTS_junction column not found in Junctions. Skipping RT-switching by SJ type plots.") + message("RT-switching Unique Junctions plot skipped: unique junction columns not in cell summary.") } # Create grouped violins for % reads with all canonical junctions by structural category (HTML) diff --git a/tests/sqanti_sc_test.py b/tests/sqanti_sc_test.py index 61935b6..0803f2e 100644 --- a/tests/sqanti_sc_test.py +++ b/tests/sqanti_sc_test.py @@ -1039,7 +1039,7 @@ def _cls_row(self, isoform, cb, fl, structural_category, } def _junc_row(self, isoform, junction_category="known", canonical="canonical", - rts=False): + rts=False, start="1000", end="2000"): """Return one junction row dict (isoform key only; CB/FL come from cls join).""" return { "isoform": isoform, @@ -1049,8 +1049,8 @@ def _junc_row(self, isoform, junction_category="known", canonical="canonical", "junction_number": "1", "chrom": "chr1", "strand": "+", - "genomic_start_coord": "1000", - "genomic_end_coord": "2000", + "genomic_start_coord": start, + "genomic_end_coord": end, } def _run(self, mock_args, tmpdir, cls_rows, junc_rows=None): @@ -1241,6 +1241,94 @@ def test_multi_junction_isoform_each_junction_gets_fl_count(self, mock_args, tmp f"got {cb1['Known_canonical_junctions']}" ) + # ------------------------------------------------------------------ + # Test 4b — RT-switching all-junction numerators are FL-weighted; the + # denominators are the existing per-SJ-type totals + # ------------------------------------------------------------------ + + def test_rts_all_junction_counts_fl_weighted(self, mock_args, tmpdir): + """ + CB1 has iso1 (FSM, FL=1) with a known_canonical RTS junction and + iso2 (NIC, FL=9) with a novel_non_canonical non-RTS junction. + RTS_Known_canonical_junctions = 1 (FL=1, RTS) + Known_canonical_junctions = 1 (denominator) + RTS_Novel_non_canonical_junctions = 0 (not RTS) + Novel_non_canonical_junctions = 9 (denominator) + """ + cls_rows = [ + self._cls_row("iso1", "CB1", "1", "full-splice_match"), + self._cls_row("iso2", "CB1", "9", "novel_in_catalog"), + ] + junc_rows = [ + self._junc_row("iso1", "known", "canonical", rts=True), + self._junc_row("iso2", "novel", "non_canonical", rts=False, start="3000", end="4000"), + ] + summary = self._run(mock_args, tmpdir, cls_rows, junc_rows) + cb1 = summary[summary["CB"] == "CB1"].iloc[0] + + assert cb1["RTS_Known_canonical_junctions"] == 1, ( + f"Expected RTS_Known_canonical_junctions=1 (FL=1, RTS), got {cb1['RTS_Known_canonical_junctions']}" + ) + assert cb1["Known_canonical_junctions"] == 1 + assert cb1["RTS_Novel_non_canonical_junctions"] == 0, ( + f"Expected RTS_Novel_non_canonical_junctions=0 (not RTS), got {cb1['RTS_Novel_non_canonical_junctions']}" + ) + assert cb1["Novel_non_canonical_junctions"] == 9 + + # ------------------------------------------------------------------ + # Test 4c — unique-junction counts collapse the SAME genomic intron + # shared across transcripts (set union, NOT FL-weighted sum) + # ------------------------------------------------------------------ + + def test_unique_junction_counts_collapse_shared_introns(self, mock_args, tmpdir): + """ + iso1 (FL=3) and iso2 (FL=5) both carry the SAME intron chr1:1000-2000, + both expressed in CB1. The FL-weighted total counts it twice (3+5=8); + the distinct-intron count collapses it to a single junction. RTS if ANY + supporting row is RTS. + """ + cls_rows = [ + self._cls_row("iso1", "CB1", "3", "full-splice_match"), + self._cls_row("iso2", "CB1", "5", "full-splice_match"), + ] + junc_rows = [ + self._junc_row("iso1", "known", "canonical", rts=True), # chr1:1000-2000 + self._junc_row("iso2", "known", "canonical", rts=False), # chr1:1000-2000 (same intron) + ] + summary = self._run(mock_args, tmpdir, cls_rows, junc_rows) + cb1 = summary[summary["CB"] == "CB1"].iloc[0] + + assert cb1["Known_canonical_junctions"] == 8, ( + f"FL-weighted total should double-count the shared intron (3+5), got {cb1['Known_canonical_junctions']}" + ) + assert cb1["unique_Known_canonical_junctions"] == 1, ( + f"Distinct-intron count should collapse the shared intron to 1, got {cb1['unique_Known_canonical_junctions']}" + ) + assert cb1["RTS_unique_Known_canonical_junctions"] == 1, ( + f"Distinct intron is RTS via iso1, got {cb1['RTS_unique_Known_canonical_junctions']}" + ) + + def test_unique_junction_counts_distinct_introns(self, mock_args, tmpdir): + """ + iso1 (FL=4, 3 exons) has two DIFFERENT known_canonical introns; both are + distinct so unique_Known_canonical_junctions=2 while the FL-weighted total + is 2 junctions x FL=4 = 8. + """ + cls_rows = [ + self._cls_row("iso1", "CB1", "4", "full-splice_match", exons=3), + ] + junc_rows = [ + self._junc_row("iso1", "known", "canonical", start="1000", end="2000"), + self._junc_row("iso1", "known", "canonical", start="3000", end="4000"), + ] + summary = self._run(mock_args, tmpdir, cls_rows, junc_rows) + cb1 = summary[summary["CB"] == "CB1"].iloc[0] + + assert cb1["unique_Known_canonical_junctions"] == 2, ( + f"Two distinct introns should count as 2, got {cb1['unique_Known_canonical_junctions']}" + ) + assert cb1["Known_canonical_junctions"] == 8 + # ------------------------------------------------------------------ # Test 5 — same isoform appearing in multiple cells is weighted # independently per cell @@ -1476,6 +1564,64 @@ def test_reads_mode_uses_one_count_per_row(self, mock_args, tmpdir): # Novel-gene bin columns are not produced at all. assert not [c for c in summary.columns if c.startswith("novel_gene_reads_bin_")] + # ------------------------------------------------------------------ + # Test 7 — reads-mode RT-switching all/unique junction counts. Each read's + # junction carries a CB (here merged from the classification table). + # ------------------------------------------------------------------ + + def test_reads_mode_rts_and_unique_junctions(self, mock_args, tmpdir): + """ + CB1 has 3 reads; r1 and r2 share intron chr1:1000-2000 (r1 RTS), r3 has a + distinct intron chr1:3000-4000. All known_canonical. + Known_canonical_junctions = 3 (one per junction row) + RTS_Known_canonical_junctions = 1 (r1) + unique_Known_canonical_junctions = 2 (two distinct introns) + RTS_unique_Known_canonical_junctions = 1 (the shared intron is RTS via r1) + """ + out_dir = str(tmpdir.join("output_reads_junc")) + file_acc, sampleID = "f1", "s1" + sample_dir = os.path.join(out_dir, file_acc) + os.makedirs(sample_dir, exist_ok=True) + prefix = os.path.join(sample_dir, sampleID) + + mock_args.mode = "reads" + mock_args.out_dir = out_dir + + cls = pd.DataFrame([ + {"isoform": rid, "CB": "CB1", "structural_category": "full-splice_match", + "associated_gene": "geneA", "associated_transcript": "txA", + "exons": 2, "length": 500, "ref_length": 600, "chrom": "chr1"} + for rid in ("r1", "r2", "r3") + ]) + cls.to_csv(f"{prefix}_classification.txt", sep="\t", index=False) + + # Junctions file has no CB column -> merged from classification by isoform. + junc = pd.DataFrame([ + {"isoform": "r1", "junction_category": "known", "canonical": "canonical", + "RTS_junction": "True", "chrom": "chr1", "strand": "+", + "genomic_start_coord": "1000", "genomic_end_coord": "2000"}, + {"isoform": "r2", "junction_category": "known", "canonical": "canonical", + "RTS_junction": "False", "chrom": "chr1", "strand": "+", + "genomic_start_coord": "1000", "genomic_end_coord": "2000"}, + {"isoform": "r3", "junction_category": "known", "canonical": "canonical", + "RTS_junction": "False", "chrom": "chr1", "strand": "+", + "genomic_start_coord": "3000", "genomic_end_coord": "4000"}, + ]) + junc.to_csv(f"{prefix}_junctions.txt", sep="\t", index=False) + + design_df = pd.DataFrame({"sampleID": [sampleID], "file_acc": [file_acc]}) + calculate_metrics_per_cell(mock_args, design_df) + + summary = pd.read_csv( + f"{prefix}_SQANTI_cell_summary.txt.gz", sep="\t", compression="gzip" + ) + cb1 = summary[summary["CB"] == "CB1"].iloc[0] + + assert cb1["Known_canonical_junctions"] == 3 + assert cb1["RTS_Known_canonical_junctions"] == 1 + assert cb1["unique_Known_canonical_junctions"] == 2 + assert cb1["RTS_unique_Known_canonical_junctions"] == 1 + # ============================================================================== # Export Scanpy / Seurat Tests