diff --git a/src/cell_metrics.py b/src/cell_metrics.py index e59d9a3..4e1ce07 100644 --- a/src/cell_metrics.py +++ b/src/cell_metrics.py @@ -207,6 +207,10 @@ def weighted_median(): summary['Annotated_genes'] = distinct(gene_code, anno_mask) summary['Novel_genes'] = distinct(gene_code, novel_gene) + # FL (read) counts from annotated vs novel genes -- numerators the report + # used to re-derive by re-exploding CB/FL (Annotated/Novel reads % plot). + summary['Annotated_genes_reads'] = msum(anno_mask) + summary['Novel_genes_reads'] = msum(novel_gene) # ---- junctions: per-isoform type counts weighted by per-cell FL ---- junc_types = ['known_canonical', 'known_non_canonical', 'novel_canonical', 'novel_non_canonical'] @@ -249,13 +253,16 @@ def weighted_median(): numer = msum(m_cat(cat) & (subcat == lv)) summary[f"{cat_to_tag[cat]}_{lv.replace('-', '_')}_prop"] = safe_prop(numer, denom).fillna(0) - # ---- gene read-count bins (per gene type; boundaries 1 / 2-4 / 5-9 / 10+) ---- + # ---- gene read-count bins (genes per bin, by gene type) ---- + # Store raw gene COUNTS per bin (not percentages): the report derives every + # percentage it needs (per-type grouped, all-genes-combined, annotated-only) + # from these and no longer re-explodes CB/FL. Boundaries: 1 / 2-5 / 6-9 / >=10. def gene_bins(gene_type_mask, prefix): - labels = [(f"{prefix}_bin1_perc", None), (f"{prefix}_bin2_4_perc", None), - (f"{prefix}_bin5_9_perc", None), (f"{prefix}_bin10plus_perc", None)] + bin_labels = [f"{prefix}_gene_reads_bin_1", f"{prefix}_gene_reads_bin_2_5", + f"{prefix}_gene_reads_bin_6_9", f"{prefix}_gene_reads_bin_10plus"] sel = gene_type_mask[row] & (gene_code[row] >= 0) if not sel.any(): - for lab, _ in labels: + for lab in bin_labels: summary[lab] = 0 return g = gene_code[row][sel].astype(np.int64) @@ -269,16 +276,14 @@ def gene_bins(gene_type_mask, prefix): idx_start = np.nonzero(first)[0] sums = np.add.reduceat(w, idx_start) cc = (key[idx_start] % Cn).astype(np.int64) - total_genes = pd.Series(np.bincount(cc, minlength=Cn).astype(float), index=cells_index) - def binperc(binmask, lab): - b = pd.Series(np.bincount(cc[binmask], minlength=Cn).astype(float), index=cells_index) - summary[lab] = safe_prop(b, total_genes).fillna(0) + def bincount_col(binmask, lab): + summary[lab] = pd.Series(np.bincount(cc[binmask], minlength=Cn).astype(float), index=cells_index) - binperc(sums == 1, f"{prefix}_bin1_perc") - binperc((sums >= 2) & (sums <= 4), f"{prefix}_bin2_4_perc") - binperc((sums >= 5) & (sums <= 9), f"{prefix}_bin5_9_perc") - binperc(sums >= 10, f"{prefix}_bin10plus_perc") + bincount_col(sums == 1, bin_labels[0]) + bincount_col((sums >= 2) & (sums <= 5), bin_labels[1]) + bincount_col((sums >= 6) & (sums <= 9), bin_labels[2]) + bincount_col(sums >= 10, bin_labels[3]) gene_bins(anno_mask, 'anno') gene_bins(novel_gene, 'novel') @@ -507,7 +512,12 @@ def lenbin_masks(base_mask): # frame. Produces the same columns as the reads-mode branch below. summary = _isoforms_summary(cls, junc) summary = summary.reset_index() - summary = summary.rename(columns={'total_reads': 'Transcripts_in_cell', 'total_reads_no_monoexon': 'total_transcripts_no_monoexon', 'MT_reads_count': 'MT_transcripts_count'}) + # In isoforms mode the per-cell counter is transcript-level (sum of FL), + # so the gene read-count columns follow the tool's convention and are + # named *_transcripts / *_gene_transcripts_bin_* (cf. Transcripts_in_cell). + gene_tx = {c: c.replace('_reads', '_transcripts') for c in summary.columns + if 'genes_reads' in c or 'gene_reads_bin' in c} + summary = summary.rename(columns={'total_reads': 'Transcripts_in_cell', 'total_reads_no_monoexon': 'total_transcripts_no_monoexon', 'MT_reads_count': 'MT_transcripts_count', **gene_tx}) for c in summary.columns[1:]: summary[c] = pd.to_numeric(summary[c], errors='coerce').fillna(0) try: @@ -587,6 +597,9 @@ def lenbin_masks(base_mask): anno = (~cls_valid['associated_gene'].fillna('').str.startswith('novel')) summary['Annotated_genes'] = cls_valid[anno].groupby('CB')['associated_gene'].nunique().reindex(summary.index, fill_value=0) summary['Novel_genes'] = cls_valid[~anno].groupby('CB')['associated_gene'].nunique().reindex(summary.index, fill_value=0) + # Read counts from annotated vs novel genes (report reads these directly). + summary['Annotated_genes_reads'] = cls_valid[anno].groupby('CB')['_count'].sum().reindex(summary.index, fill_value=0) + summary['Novel_genes_reads'] = cls_valid[~anno].groupby('CB')['_count'].sum().reindex(summary.index, fill_value=0) if not junc.empty: junc_types = ['known_canonical', 'known_non_canonical', 'novel_canonical', 'novel_non_canonical'] @@ -672,25 +685,28 @@ def subkey(cat, lv): numer = tbl.get(lv, pd.Series(0, index=summary.index)).reindex(summary.index, fill_value=0) summary[subkey(cat, lv.replace('-', '_'))] = safe_prop(numer, denom).fillna(0) + # Gene read-count bins as raw gene COUNTS (boundaries 1 / 2-5 / 6-9 / >=10); + # the report reads these directly and derives whatever percentages each plot needs. gene_counts = cls_valid.groupby(['CB','associated_gene'])['_count'].sum().rename('read_count').reset_index() gene_counts['gene_type'] = np.where(gene_counts['associated_gene'].fillna('').str.startswith('novel'), 'novel', 'annotated') + bin_suffixes = ['gene_reads_bin_1', 'gene_reads_bin_2_5', 'gene_reads_bin_6_9', 'gene_reads_bin_10plus'] bins = gene_counts.groupby(['CB','gene_type']).agg( - bin1_count=('read_count', lambda s: (s == 1).sum()), - bin2_4_count=('read_count', lambda s: ((s >= 2) & (s <= 4)).sum()), - bin5_9_count=('read_count', lambda s: ((s >= 5) & (s <= 9)).sum()), - bin10plus_count=('read_count', lambda s: (s >= 10).sum()), - total_genes_in_type=('associated_gene','nunique') + gene_reads_bin_1=('read_count', lambda s: (s == 1).sum()), + gene_reads_bin_2_5=('read_count', lambda s: ((s >= 2) & (s <= 5)).sum()), + gene_reads_bin_6_9=('read_count', lambda s: ((s >= 6) & (s <= 9)).sum()), + gene_reads_bin_10plus=('read_count', lambda s: (s >= 10).sum()), ).reset_index() - def bin_props(df, gene_kind, out_prefix): + def bin_counts(df, gene_kind, out_prefix): out = pd.DataFrame(index=summary.index) keyed = df[df['gene_type'] == gene_kind].set_index('CB') if not df.empty else pd.DataFrame(index=summary.index) - for label, src in [(f"{out_prefix}_bin1_perc", 'bin1_count'), (f"{out_prefix}_bin2_4_perc", 'bin2_4_count'), (f"{out_prefix}_bin5_9_perc", 'bin5_9_count'), (f"{out_prefix}_bin10plus_perc", 'bin10plus_count')]: - if not keyed.empty and src in keyed.columns: - out[label] = safe_prop(keyed[src].reindex(summary.index, fill_value=0), keyed['total_genes_in_type'].reindex(summary.index, fill_value=0)).fillna(0) + for suffix in bin_suffixes: + col = f"{out_prefix}_{suffix}" + if not keyed.empty and suffix in keyed.columns: + out[col] = keyed[suffix].reindex(summary.index, fill_value=0) else: - out[label] = 0 + out[col] = 0 return out - summary = summary.join(bin_props(bins, 'annotated', 'anno')).join(bin_props(bins, 'novel', 'novel')) + summary = summary.join(bin_counts(bins, 'annotated', 'anno')).join(bin_counts(bins, 'novel', 'novel')) if args.mode != 'isoforms': gene_ujc = cls_valid[cls_valid['exons'] > 1].groupby(['CB','associated_gene'])['jxn_string'].nunique().rename('ujc_count').reset_index() @@ -977,7 +993,9 @@ def compute_lenbins_by_cb(df_group): summary = summary.reset_index() if args.mode == 'isoforms': - summary = summary.rename(columns={'total_reads': 'Transcripts_in_cell', 'total_reads_no_monoexon': 'total_transcripts_no_monoexon', 'MT_reads_count': 'MT_transcripts_count'}) + gene_tx = {c: c.replace('_reads', '_transcripts') for c in summary.columns + if 'genes_reads' in c or 'gene_reads_bin' in c} + summary = summary.rename(columns={'total_reads': 'Transcripts_in_cell', 'total_reads_no_monoexon': 'total_transcripts_no_monoexon', 'MT_reads_count': 'MT_transcripts_count', **gene_tx}) else: summary = summary.rename(columns={'total_reads': 'Reads_in_cell', 'total_UMI': 'UMIs_in_cell'}) diff --git a/src/report_assets/SQANTI-sc_report.R b/src/report_assets/SQANTI-sc_report.R index 41f60c5..b53ba56 100644 --- a/src/report_assets/SQANTI-sc_report.R +++ b/src/report_assets/SQANTI-sc_report.R @@ -1392,68 +1392,14 @@ generate_sqantisc_plots <- function(SQANTI_cell_summary, Classification_file, Ju classification_valid <- Classification_file[Classification_file$CB != "unassigned" & !is.na(Classification_file$CB), ] if (nrow(classification_valid) > 0) { - # Function to expand FL and CB columns into a long format for correct counting per cell - expand_isoform_counts <- function(df, mode) { - if (mode == "reads") { - return(df %>% group_by(CB) %>% summarise(count = n(), .groups = "drop")) - } else { - # Isoforms mode: Each row has comma-separated FL (counts) and CB (barcodes) - # We need to split them and sum counts per barcode - - # Initialize lists to store expanded data - all_cbs <- character() - all_counts <- numeric() - - # Iterate through rows (this might be slow for huge files, but safe) - # A vectorised approach would be better if possible, but strsplit returns list - fl_list <- strsplit(as.character(df$FL), ",") - cb_list <- strsplit(as.character(df$CB), ",") - - # Check if lengths match (they should) - if (length(fl_list) != length(cb_list)) { - stop("Mismatch in row counts between FL and CB columns") - } - - # Use mapply to create a data frame of all counts - # This creates a list of data frames, one per isoform - expanded_list <- mapply(function(fl, cb) { - if (length(fl) != length(cb)) { - # Warning or skip? For now, we assume they match as per SQANTI specs - return(NULL) - } - data.frame(CB = cb, count = as.numeric(fl), stringsAsFactors = FALSE) - }, fl_list, cb_list, SIMPLIFY = FALSE) - - # Bind all tiny data frames - long_df <- do.call(rbind, expanded_list) - - # Now group by CB and sum - return(long_df %>% group_by(CB) %>% summarise(count = sum(count, na.rm = TRUE), .groups = "drop")) - } - } - - annotated_reads_per_cell <- classification_valid %>% - filter(!grepl("^novel", associated_gene)) - - annotated_reads_per_cell <- expand_isoform_counts(annotated_reads_per_cell, mode) %>% - rename(Annotated_genes_reads = count) - - novel_reads_per_cell <- classification_valid %>% - filter(grepl("^novel", associated_gene)) - - novel_reads_per_cell <- expand_isoform_counts(novel_reads_per_cell, mode) %>% - rename(Novel_genes_reads = count) - - SQANTI_cell_summary <- SQANTI_cell_summary %>% - left_join(annotated_reads_per_cell, by = "CB") %>% - left_join(novel_reads_per_cell, by = "CB") - - SQANTI_cell_summary$Annotated_genes_reads[is.na(SQANTI_cell_summary$Annotated_genes_reads)] <- 0 - SQANTI_cell_summary$Novel_genes_reads[is.na(SQANTI_cell_summary$Novel_genes_reads)] <- 0 - - # Revert to original denominator (Total Transcripts in Cell) now that numerators are correct - SQANTI_cell_summary$Annotated_reads_perc <- 100 * SQANTI_cell_summary$Annotated_genes_reads / SQANTI_cell_summary[[count_col]] - SQANTI_cell_summary$Novel_reads_perc <- 100 * SQANTI_cell_summary$Novel_genes_reads / SQANTI_cell_summary[[count_col]] + # Annotated/Novel per-cell counts are provided directly by cell_metrics, so we + # no longer re-explode CB/FL here -- we just read the columns and take the ratio. + # Column names follow the mode convention: Annotated_genes_reads (reads mode) / + # Annotated_genes_transcripts (isoforms mode), via entity_label_plural_lower. + anno_genes_col <- paste0("Annotated_genes_", entity_label_plural_lower) + novel_genes_col <- paste0("Novel_genes_", entity_label_plural_lower) + SQANTI_cell_summary$Annotated_reads_perc <- 100 * SQANTI_cell_summary[[anno_genes_col]] / SQANTI_cell_summary[[count_col]] + SQANTI_cell_summary$Novel_reads_perc <- 100 * SQANTI_cell_summary[[novel_genes_col]] / SQANTI_cell_summary[[count_col]] SQANTI_cell_summary$Annotated_reads_perc <- ifelse(is.na(SQANTI_cell_summary$Annotated_reads_perc) | is.infinite(SQANTI_cell_summary$Annotated_reads_perc), 0, SQANTI_cell_summary$Annotated_reads_perc) SQANTI_cell_summary$Novel_reads_perc <- ifelse(is.na(SQANTI_cell_summary$Novel_reads_perc) | is.infinite(SQANTI_cell_summary$Novel_reads_perc), 0, SQANTI_cell_summary$Novel_reads_perc) @@ -1497,68 +1443,35 @@ generate_sqantisc_plots <- function(SQANTI_cell_summary, Classification_file, Ju #################################################################### # Define gene read-count bins and labels - gene_bin_label <- function(n) { - if (is.na(n)) { - return(NA_character_) - } - if (n == 1) { - return("1") - } - if (n >= 2 && n <= 5) { - return("2-5") - } - if (n >= 6 && n <= 9) { - return("6-9") - } - return(">=10") - } gene_bin_levels <- c("1", "2-5", "6-9", ">=10") - # Build per-cell per-gene read counts from classification. - # In isoforms mode, must explode the comma-separated CB/FL columns so each - # (cell, isoform) pair is weighted by its FL count, then sum per (CB, gene). - # In reads mode, each row is one read so n() is correct. - if (mode == "isoforms" && "FL" %in% colnames(Classification_file) && "CB" %in% colnames(Classification_file)) { - genes_by_cb_base <- Classification_file %>% - filter(!is.na(CB), CB != "unassigned", !is.na(associated_gene)) %>% - select(CB, FL, associated_gene) - - genes_by_cb_base$CB_raw <- as.character(genes_by_cb_base$CB) - genes_by_cb_base$FL_raw <- as.character(genes_by_cb_base$FL) - genes_by_cb_base <- tidyr::separate_rows(genes_by_cb_base, CB_raw, FL_raw, sep = ",") - genes_by_cb_base$FL_num <- suppressWarnings(as.numeric(trimws(genes_by_cb_base$FL_raw))) - genes_by_cb_base$FL_num[is.na(genes_by_cb_base$FL_num) | genes_by_cb_base$FL_num < 0] <- 0 - genes_by_cb_base$CB_clean <- trimws(genes_by_cb_base$CB_raw) - - genes_by_cb <- genes_by_cb_base %>% - filter(CB_clean != "" & CB_clean != "unassigned" & FL_num > 0) %>% - group_by(CB = CB_clean, associated_gene) %>% - summarise(reads_per_gene = sum(FL_num), .groups = "drop") %>% - mutate( - gene_type = ifelse(grepl("^novel", associated_gene), "Novel", "Annotated"), - bin = vapply(reads_per_gene, gene_bin_label, character(1)) - ) %>% - filter(!is.na(bin)) - } else { - genes_by_cb <- Classification_file %>% - filter(!is.na(CB), CB != "unassigned", !is.na(associated_gene)) %>% - group_by(CB, associated_gene) %>% - summarise(reads_per_gene = n(), .groups = "drop") %>% - mutate( - gene_type = ifelse(grepl("^novel", associated_gene), "Novel", "Annotated"), - bin = vapply(reads_per_gene, gene_bin_label, character(1)) - ) %>% - filter(!is.na(bin)) - } + # Gene read-count bins come directly from cell_metrics (both modes): the summary + # holds per-cell gene counts per bin. Column names follow the mode convention -- + # anno/novel_gene_reads_bin_* (reads) or _gene_transcripts_bin_* (isoforms), via + # entity_label_plural_lower. Build a per-(CB, gene_type, bin) num_genes table from + # those columns -- no CB/FL explosion. + bin_suffix <- paste0("_gene_", entity_label_plural_lower, "_bin_") + bin_defs <- list( + c("1", paste0(bin_suffix, "1")), c("2-5", paste0(bin_suffix, "2_5")), + c("6-9", paste0(bin_suffix, "6_9")), c(">=10", paste0(bin_suffix, "10plus")) + ) + read_bins_counts <- do.call(rbind, lapply(list(c("Annotated", "anno"), c("Novel", "novel")), function(gp) { + do.call(rbind, lapply(bin_defs, function(bd) { + data.frame( + CB = SQANTI_cell_summary$CB, + gene_type = gp[1], + bin = bd[1], + num_genes = as.numeric(SQANTI_cell_summary[[paste0(gp[2], bd[2])]]), + stringsAsFactors = FALSE + ) + })) + })) # Percent of genes per bin within each CB and gene type - read_bins_data <- genes_by_cb %>% - group_by(CB, gene_type, bin) %>% - summarise(num_genes = n(), .groups = "drop") %>% + read_bins_data <- read_bins_counts %>% group_by(CB, gene_type) %>% - mutate(percentage = 100 * num_genes / sum(num_genes)) %>% - ungroup() %>% - tidyr::complete(CB, gene_type, bin = gene_bin_levels, fill = list(num_genes = 0, percentage = 0)) + mutate(percentage = ifelse(sum(num_genes) > 0, 100 * num_genes / sum(num_genes), 0)) %>% + ungroup() read_bins_data$bin <- factor(read_bins_data$bin, levels = gene_bin_levels) read_bins_data$gene_type <- factor(read_bins_data$gene_type, levels = c("Annotated", "Novel")) @@ -1582,28 +1495,24 @@ generate_sqantisc_plots <- function(SQANTI_cell_summary, Classification_file, Ju ) } - # Combined (all genes together): one violin per bin + # Combined: one violin per bin. Reads mode shows annotated genes only; isoforms + # mode shows all genes (annotated + novel) summed per bin. if (mode == "reads") { - # Filter for Annotated genes only - read_bins_all <- genes_by_cb %>% + read_bins_all <- read_bins_counts %>% filter(gene_type == "Annotated") %>% - group_by(CB, bin) %>% - summarise(num_genes = n(), .groups = "drop") %>% group_by(CB) %>% - mutate(percentage = 100 * num_genes / sum(num_genes)) %>% + mutate(percentage = ifelse(sum(num_genes) > 0, 100 * num_genes / sum(num_genes), 0)) %>% ungroup() %>% - tidyr::complete(CB, bin = gene_bin_levels, fill = list(num_genes = 0, percentage = 0)) + select(CB, bin, num_genes, percentage) plot_title_all <- paste("Distribution of Annotated Genes by", entity_label, "Count Bins Across Cells") } else { - # All genes (Annotated + Novel) - read_bins_all <- genes_by_cb %>% + read_bins_all <- read_bins_counts %>% group_by(CB, bin) %>% - summarise(num_genes = n(), .groups = "drop") %>% + summarise(num_genes = sum(num_genes), .groups = "drop") %>% group_by(CB) %>% - mutate(percentage = 100 * num_genes / sum(num_genes)) %>% - ungroup() %>% - tidyr::complete(CB, bin = gene_bin_levels, fill = list(num_genes = 0, percentage = 0)) + mutate(percentage = ifelse(sum(num_genes) > 0, 100 * num_genes / sum(num_genes), 0)) %>% + ungroup() plot_title_all <- paste("Distribution of Genes by", entity_label, "Count Bins Across Cells") } @@ -1635,14 +1544,12 @@ generate_sqantisc_plots <- function(SQANTI_cell_summary, Classification_file, Ju # New plot: Distribution of Known Genes by Unique Isoform Count Bins Across Cells (Isoforms mode) if (mode == "isoforms") { - iso_bins_annot <- genes_by_cb %>% + iso_bins_annot <- read_bins_counts %>% filter(gene_type == "Annotated") %>% - group_by(CB, bin) %>% - summarise(num_genes = n(), .groups = "drop") %>% group_by(CB) %>% - mutate(percentage = 100 * num_genes / sum(num_genes)) %>% + mutate(percentage = ifelse(sum(num_genes) > 0, 100 * num_genes / sum(num_genes), 0)) %>% ungroup() %>% - tidyr::complete(CB, bin = gene_bin_levels, fill = list(num_genes = 0, percentage = 0)) + select(CB, bin, num_genes, percentage) iso_bins_annot$bin <- factor(iso_bins_annot$bin, levels = gene_bin_levels) diff --git a/tests/sqanti_sc_test.py b/tests/sqanti_sc_test.py index 2ddc6a3..03cdb61 100644 --- a/tests/sqanti_sc_test.py +++ b/tests/sqanti_sc_test.py @@ -1229,6 +1229,66 @@ def test_fl_weighting_is_per_cell_independent(self, mock_args, tmpdir): assert abs(cb1["FSM_prop"] - 100.0) < 0.01 assert abs(cb2["FSM_prop"] - 100.0) < 0.01 + # ------------------------------------------------------------------ + # Test 5b — annotated/novel READ counts are FL-weighted, split by gene + # annotation (associated_gene 'novel' prefix, not category) + # ------------------------------------------------------------------ + + def test_annotated_novel_reads_are_fl_weighted(self, mock_args, tmpdir): + """ + Annotated_genes_transcripts / Novel_genes_transcripts are per-cell FL sums split by + whether associated_gene starts with 'novel'. + CB1: geneA FL=3 (annotated) + novelGene FL=1 (novel) + CB2: geneA FL=10 + geneA FL=4 (annotated) + novelGene FL=6 (novel) + """ + cls_rows = [ + self._cls_row("iso1", "CB1,CB2", "3,10", "full-splice_match", associated_gene="geneA"), + self._cls_row("iso2", "CB2", "4", "full-splice_match", associated_gene="geneA"), + self._cls_row("iso3", "CB1,CB2", "1,6", "novel_in_catalog", associated_gene="novelGene_1"), + ] + summary = self._run(mock_args, tmpdir, cls_rows) + cb1 = summary[summary["CB"] == "CB1"].iloc[0] + cb2 = summary[summary["CB"] == "CB2"].iloc[0] + + assert cb1["Annotated_genes_transcripts"] == 3 + assert cb1["Novel_genes_transcripts"] == 1 + assert cb2["Annotated_genes_transcripts"] == 14 # 10 + 4 + assert cb2["Novel_genes_transcripts"] == 6 + + # ------------------------------------------------------------------ + # Test 5c — gene read-count bins store gene COUNTS in the report's + # boundaries (1 / 2-5 / 6-9 / >=10), FL-summed per gene + # ------------------------------------------------------------------ + + def test_gene_read_count_bins_are_fl_weighted_counts(self, mock_args, tmpdir): + """ + Each bin column counts genes whose per-cell FL sum falls in that bin. + CB1: geneA=3+2=5 (anno -> 2-5), novelGene=1 (novel -> 1) + CB2: geneA=10 (anno -> >=10), geneB=4 (anno -> 2-5), novelGene=6 (novel -> 6-9) + """ + cls_rows = [ + self._cls_row("iso1", "CB1,CB2", "3,10", "full-splice_match", associated_gene="geneA"), + self._cls_row("iso2", "CB1", "2", "full-splice_match", associated_gene="geneA"), + self._cls_row("iso3", "CB2", "4", "full-splice_match", associated_gene="geneB"), + self._cls_row("iso4", "CB1,CB2", "1,6", "novel_in_catalog", associated_gene="novelGene_1"), + ] + summary = self._run(mock_args, tmpdir, cls_rows) + cb1 = summary[summary["CB"] == "CB1"].iloc[0] + cb2 = summary[summary["CB"] == "CB2"].iloc[0] + + # CB1 annotated: geneA=5 -> bin 2-5 only + assert cb1["anno_gene_transcripts_bin_1"] == 0 + assert cb1["anno_gene_transcripts_bin_2_5"] == 1 + assert cb1["anno_gene_transcripts_bin_6_9"] == 0 + assert cb1["anno_gene_transcripts_bin_10plus"] == 0 + # CB1 novel: novelGene=1 -> bin 1 + assert cb1["novel_gene_transcripts_bin_1"] == 1 + # CB2 annotated: geneA=10 -> >=10, geneB=4 -> 2-5 + assert cb2["anno_gene_transcripts_bin_2_5"] == 1 + assert cb2["anno_gene_transcripts_bin_10plus"] == 1 + # CB2 novel: novelGene=6 -> 6-9 + assert cb2["novel_gene_transcripts_bin_6_9"] == 1 + # ------------------------------------------------------------------ # Test 6 — reads mode uses 1 count per row regardless of FL column # ------------------------------------------------------------------ @@ -1282,6 +1342,14 @@ def test_reads_mode_uses_one_count_per_row(self, mock_args, tmpdir): f"Reads mode: expected FSM_prop≈66.7 %, got {cb1['FSM_prop']}" ) + # Gene-annotation columns are gene-based, not category-based: all 3 reads + # belong to annotated geneA (even the novel-category read), 3 reads -> bin 2-5. + assert cb1["Annotated_genes_reads"] == 3 + assert cb1["Novel_genes_reads"] == 0 + assert cb1["anno_gene_reads_bin_1"] == 0 + assert cb1["anno_gene_reads_bin_2_5"] == 1 + assert cb1["novel_gene_reads_bin_2_5"] == 0 + # ============================================================================== # Export Scanpy / Seurat Tests