import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import re
import time
import numpy as np
import pandas as pd
from IPython.display import display
import scala

T0 = time.time()
print("dandiset", scala.DANDISET, "version", scala.VERSION)
assets = scala.list_assets()
ids = [re.search(r"_cell-(.+)_icephys", p).group(1).replace("-", "_") for p in assets.path]
ex, it, gl = scala.patchseq_counts()
tab = scala.paper_meta_table()
rnaseqTools = scala.paper_module("rnaseqTools")

dandiset 000008 version 0.211014.0809

print("count tables: genes x cells", ex.shape, "| gene lengths:", gl.shape)
print("dandiset files:", len(ids), "| their cell ids in the count tables:", int(np.isin(ids, ex.columns).sum()),
      "| count-table cells without a dandiset file:", sorted(set(ex.columns) - set(ids)))
print("count values: exon", int(ex.values.min()), "to", int(ex.values.max()), "| intron", int(it.values.min()), "to", int(it.values.max()))
det_tab = ((ex.values + it.values) > 0).sum(axis=0)
print("detected genes per cell equal to the table's 'Number of genes detected':",
      int((det_tab == pd.to_numeric(tab["Number of genes detected"]).reindex(ex.columns).values).sum()), "of", ex.shape[1])
for kind, cnt, col in [("exon", ex, "Exon reads"), ("intron", it, "Intron reads")]:
    ratio = cnt.sum(axis=0) / pd.to_numeric(tab[col]).reindex(cnt.columns)
    print(f"{kind} column sums over the table's '{col}': equal for {int(np.isclose(ratio, 1).sum())} cells; ratio median {ratio.median():.3f}, "
          f"5th to 95th percentile {ratio.quantile(0.05):.3f} to {ratio.quantile(0.95):.3f}")
print("cells in the paper's table order:", bool((ex.columns == tab.index).all()))

count tables: genes x cells (42466, 1329) | gene lengths: (42466, 4)
dandiset files: 1328 | their cell ids in the count tables: 1328 | count-table cells without a dandiset file: ['20190211_sample_10']

count values: exon 0 to 1536853 | intron 0 to 253008

detected genes per cell equal to the table's 'Number of genes detected': 1329 of 1329
exon column sums over the table's 'Exon reads': equal for 14 cells; ratio median 1.035, 5th to 95th percentile 0.384 to 3.078
intron column sums over the table's 'Intron reads': equal for 4 cells; ratio median 1.000, 5th to 95th percentile 0.090 to 10.813
cells in the paper's table order: True

E = ex.values.T.astype(float)
I = it.values.T.astype(float)
ed, idp = E.sum(axis=1), I.sum(axis=1)
det = ((E + I) > 0).sum(axis=1)
print(f"cells: {E.shape[0]}")
print(f"median exonic reads {np.median(ed):,.0f} (log10 mean {np.mean(np.log10(ed)):.2f} +- {np.std(np.log10(ed)):.2f}); "
      f"median intronic reads {np.median(idp):,.0f} (log10 mean {np.mean(np.log10(idp)):.2f} +- {np.std(np.log10(idp)):.2f})")
print(f"detected genes (non-zero exonic or intronic count): mean {det.mean():,.0f}, s.d. {det.std():,.0f}; range {det.min()} to {det.max()}")
d_ids = np.isin(ex.columns, ids)
print(f"the 1,328 cells with a dandiset file: median exonic {np.median(ed[d_ids]):,.0f}, intronic {np.median(idp[d_ids]):,.0f}, "
      f"mean detected genes {det[d_ids].mean():,.0f}")
er, ir = pd.to_numeric(tab["Exon reads"]), pd.to_numeric(tab["Intron reads"])
print(f"variant, the read fields of the paper's table (equal to the dandiset's fields): median exonic {er.median():,.0f} "
      f"(log10 mean {np.mean(np.log10(er)):.2f}), median intronic {ir.median():,.0f} (log10 mean {np.mean(np.log10(ir)):.2f}); "
      f"over the 1,328 dandiset cells {er[ex.columns[d_ids]].median():,.0f} and {ir[ex.columns[d_ids]].median():,.0f}")

cells: 1329
median exonic reads 1,281,492 (log10 mean 5.99 +- 0.59); median intronic reads 687,965 (log10 mean 5.61 +- 0.76)
detected genes (non-zero exonic or intronic count): mean 9,097, s.d. 3,533; range 11 to 24833
the 1,328 cells with a dandiset file: median exonic 1,280,301, intronic 689,454, mean detected genes 9,097
variant, the read fields of the paper's table (equal to the dandiset's fields): median exonic 1,190,890 (log10 mean 5.96), median intronic 688,140 (log10 mean 5.61); over the 1,328 dandiset cells 1,190,524 and 688,002

tasic = scala.reference("tasic2018")
print("Tasic 2018 reference: cells x genes", tasic["counts"].shape, "| clusters:", len(tasic["clusterNames"]))
_, C = rnaseqTools.map_to_clusters(tasic["counts"], tasic["genes"], E, ex.index.values, tasic["clusters"], returnCmeans=True,
                                   newIntronCounts=I, referenceIntronCounts=tasic["intronCounts"],
                                   newExonLengths=gl.exon_bp.values, newIntronLengths=gl.intron_bp.values,
                                   normalizeNew=True, normalizeReference=True)
C[np.isnan(C)] = 0
lowq = np.amax(C, axis=1) < .4
n1 = int(lowq.sum())
e2 = ~lowq & (np.amax(C[:, :117], axis=1) < .4); lowq |= e2
e3 = ~lowq & np.isin(np.argmax(C, axis=1), [115, 116]); lowq |= e3
borders = [0, 29, 60, 115]
maxorder = np.stack([np.amax(C[:, borders[i]:borders[i + 1]], axis=1) for i in range(3)], axis=1)
delta = np.diff(np.sort(maxorder, axis=1)[:, -2:], axis=1).squeeze()
e4 = ~lowq & (delta < 0.02); lowq |= e4
print(f"highest correlation below 0.4: {n1}; highest neural correlation below 0.4 (non-neural): {int(e2.sum())}; "
      f"mapping to Meis2 or CR: {int(e3.sum())}; groups less than 0.02 apart: {int(e4.sum())}; passing: {int((~lowq).sum())}")
excl = tab["Exclusion reasons"].notna().values
print("of the passing cells, with an exclusion reason in the paper's table:", int((~lowq & excl).sum()),
      "| passing and not excluded:", int((~lowq & ~excl).sum()))
same = (lowq == (tab["RNA family"].values == "low quality"))
print("cells whose pass or fail agrees with the paper's table ('low quality' family):", int(same.sum()), "of", len(same))
order = np.argmax(maxorder, axis=1)
order[lowq] = -1

Tasic 2018 reference: cells x genes (23822, 3000) | clusters: 133
Using a common set of 2666 genes.

highest correlation below 0.4: 78; highest neural correlation below 0.4 (non-neural): 14; mapping to Meis2 or CR: 0; groups less than 0.02 apart: 5; passing: 1232
of the passing cells, with an exclusion reason in the paper's table: 5 | passing and not excluded: 1227
cells whose pass or fail agrees with the paper's table ('low quality' family): 1329 of 1329

sets = ["10X_cells_v2_AIBS", "10X_nuclei_v2_AIBS", "10X_cells_v3_AIBS", "10X_nuclei_v3_AIBS",
        "SmartSeq_cells_AIBS", "SmartSeq_nuclei_AIBS", "10X_nuclei_v3_Broad"]
n = E.shape[0]
maxc, assigned = {}, {}
for ds in sets:
    m = scala.reference(ds)
    K = m["neurons"]["clusterNames"].size
    mc, aa = np.full(n, np.nan), np.full(n, np.nan)
    for o, key in enumerate(["viplamp", "pvsst", "exc"]):
        s = order == o
        a_, Cm = rnaseqTools.map_to_clusters(m[key]["counts"], m[key]["genes"], E[s], ex.index.values, m[key]["clusters"],
                                             returnCmeans=True, totalClusters=K, newIntronCounts=I[s],
                                             referenceIntronCounts=m[key]["intronCounts"] if "SmartSeq" in ds else None,
                                             newExonLengths=gl.exon_bp.values, newIntronLengths=gl.intron_bp.values,
                                             normalizeNew=True, normalizeReference="SmartSeq" in ds)
        mc[s] = np.nanmax(Cm, axis=1)
        aa[s] = a_
    maxc[ds], assigned[ds] = mc, aa
    del m
mx = pd.DataFrame(maxc)
summary = pd.DataFrame({"median max correlation": mx.median(), "mean": mx.mean(), "cells": mx.notna().sum()}).sort_values("median max correlation", ascending=False)
display(summary.round(4))
a, b = assigned["SmartSeq_cells_AIBS"], assigned["10X_cells_v2_AIBS"]
ok = ~np.isnan(a) & ~np.isnan(b)
print(f"highest median: {summary.index[0]}; single-cell Smart-seq2 and single-cell 10x v2 assign the same t-type for "
      f"{np.mean(a[ok] == b[ok]):.3f} of {int(ok.sum())} cells")
better = (mx["SmartSeq_nuclei_AIBS"].values[:, None] > mx.drop(columns="SmartSeq_nuclei_AIBS").values).all(axis=1)
print(f"cells whose highest correlation is with the single-nucleus Smart-seq2 reference: {np.mean(better[ok]):.3f}")
print(f"by group: " + ", ".join(f"{g}: {np.nanmedian(mx['SmartSeq_nuclei_AIBS'].values[order == o]):.3f} (snSS) vs "
                                 f"{np.nanmedian(mx['SmartSeq_cells_AIBS'].values[order == o]):.3f} (scSS)"
                                 for o, g in enumerate(["CGE", "MGE", "excitatory"])))

Using a common set of 490 genes.

Using a common set of 493 genes.

Using a common set of 478 genes.

Using a common set of 483 genes.

Using a common set of 484 genes.

Using a common set of 472 genes.

Using a common set of 491 genes.

Using a common set of 492 genes.

Using a common set of 486 genes.

Using a common set of 488 genes.

print(f"run time: {(time.time() - T0) / 60:.1f} min")

run time: 4.6 min

Claim	Tested here
rnaseq-depth	yes
ttype-qc-counts	yes
ed2-smartseq-similarity	yes

Claim	Outcome
rnaseq-depth	reproduced
ttype-qc-counts	reproduced
ed2-smartseq-similarity	reproduced

Scala, Kobak et al. (2021): sequencing claims, on the paper's count tables¶

rnaseq-depth¶

ttype-qc-counts¶

ed2-smartseq-similarity¶

Summary¶

	median max correlation	mean	cells
SmartSeq_nuclei_AIBS	0.7031	0.6892	1232
10X_cells_v2_AIBS	0.6230	0.5861	1232
SmartSeq_cells_AIBS	0.6210	0.6181	1232
10X_cells_v3_AIBS	0.5779	0.5357	1232
10X_nuclei_v3_Broad	0.5562	0.5133	1232
10X_nuclei_v2_AIBS	0.5297	0.4887	1232
10X_nuclei_v3_AIBS	0.4969	0.4554	1232