<jupyter_start><jupyter_text>Make old v young<jupyter_code>import glob
import anndata
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.patches as mpatches
import scanpy as sc
from scipy import stats

from collections import OrderedDict
from sklearn.decomposition import TruncatedSVD
from sklearn.manifold import TSNE
from sklearn.preprocessing import scale

from sklearn.cluster import KMeans
from sklearn.preprocessing import normalize
from sklearn.preprocessing import LabelEncoder
from sklearn.neighbors import NeighborhoodComponentsAnalysis
from matplotlib import cm
from matplotlib.lines import Line2D

def nd(arr):
    return np.asarray(arr).reshape(-1)
def yex(ax):
    lims = [np.min([ax.get_xlim(), ax.get_ylim()]),
            np.max([ax.get_xlim(), ax.get_ylim()])]

    # now plot both limits against eachother
    ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)
    ax.set_aspect('equal')
    ax.set_xlim(lims)
    ax.set_ylim(lims)
    return ax

def trim_axs(axs, N):
    """little helper to massage the axs list to have correct length..."""
    axs = axs.flat
    for ax in axs[N:]:
        ax.remove()
    return axs[:N]

import warnings
warnings.filterwarnings('ignore')

fsize=20

plt.rcParams.update({'font.size': fsize})
%config InlineBackend.figure_format = 'retina'
files = glob.glob("/home/sina/projects/mouse_lung/SRR*/out/counts_filtered/adata.h5ad")
f = files[0]
%%time
ads = []
names = []
for fn, f in enumerate(files):
    names.append(f.split("/")[-4])
    ads.append(anndata.read_h5ad(f))
adata = ads[0].concatenate(*ads[1:], batch_key="SRR", batch_categories=names)
age_map = {
    "SRR8426358": 24,
    "SRR8426359": 24,
    "SRR8426362": 24,
    "SRR8426363": 24,
    "SRR8426364": 24,
    "SRR8426369": 24,
    "SRR8426370": 24,
    "SRR8426360": 3,
    "SRR8426361": 3,
    "SRR8426365": 3,
    "SRR8426366": 3,
    "SRR8426367": 3,
    "SRR8426368": 3,
    "SRR8426371": 3,
    "SRR8426372": 3
}
adata.obs["age"] = adata.obs.SRR.map(age_map)
adata.var["gene_id"] = adata.var.index.values

t2g = pd.read_csv("../reference/t2g.txt", header=None, names=["tid", "gene_id", "gene_name"], sep="\t")
t2g.index = t2g.gene_id
t2g = t2g.loc[~t2g.index.duplicated(keep='first')]

adata.var["gene_name"] = adata.var.gene_id.map(t2g["gene_name"])
adata.var.index = adata.var["gene_name"]
adata<jupyter_output><empty_output><jupyter_text>Populate anndata<jupyter_code># Important constants we will use throughout the notebook

num_TSNE = 2
state = 42
metric = "euclidean"
n_neighbors = 30
num_PCA = 50
num_NCA = 10

# Filtering criteria
cell_threshold = 0
gene_threshold = 3

mito_criteria = 5

n_top_genes = 6000

n_bins = 20

flavor="seurat"
adata.obs["cell_counts"] = adata.X.sum(axis=1)
adata.var["gene_counts"] = nd(adata.X.sum(axis=0))

adata.obs["n_genes"] = nd((adata.X>0).sum(axis=1))
%%time
# n_cells needs to be computed on a per SRR basis
unique = np.unique(adata.obs.SRR.values)
assignments = adata.obs.SRR.values

mtx = adata.X
n_cells = np.zeros((unique.shape[0], mtx.shape[1]))

for sn, s in enumerate(unique):
    mask = s == assignments
    n_cells[sn] = nd((mtx[mask]>0).sum(axis=0))
adata.var["min_cells"] = np.amin(n_cells, axis=0)
mito_genes = adata.var_names.str.startswith('mt-')
adata.obs["percent_mito"] = adata[:,mito_genes].X.sum(axis=1)/adata.X.sum(axis=1)*100<jupyter_output><empty_output><jupyter_text>Knee plot<jupyter_code>expected_num_cells = 35000
knee = np.sort(nd(adata.X.sum(axis=1)))[::-1]

fig, ax = plt.subplots(figsize=(5, 5))

x = knee
y = range(len(knee))

ax.loglog(x, y, linewidth=5, color="g")

ax.axvline(x=knee[expected_num_cells], linewidth=3, color="k")
ax.axhline(y=expected_num_cells, linewidth=3, color="k")

ax.set_xlabel("UMI Counts")
ax.set_ylabel("Set of Barcodes")

plt.show()<jupyter_output><empty_output><jupyter_text>Library saturation<jupyter_code># Create a plot showing genes detected as a function of UMI counts.
fig, ax = plt.subplots(figsize=(5, 5))

x = nd(adata.X.sum(axis=1))
y = nd(np.sum(adata.X>0, axis=1))

ax.scatter(x, y, color="green", alpha=0.25)

ax.set_xlabel("UMI Counts")
ax.set_ylabel("Genes Detected")
ax.set_xscale('log')
ax.set_yscale('log')

ax.set_xlim(1)
ax.set_ylim(1)

plt.show()
fig, ax = plt.subplots(figsize=(5,5))

#histogram definition
bins = [1500, 1500] # number of bins

# histogram the data
hh, locx, locy = np.histogram2d(x, y, bins=bins)

# Sort the points by density, so that the densest points are plotted last
z = np.array([hh[np.argmax(a<=locx[1:]),np.argmax(b<=locy[1:])] for a,b in zip(x,y)])
idx = z.argsort()
x2, y2, z2 = x[idx], y[idx], z[idx]

s = ax.scatter(x2, y2, c=z2, cmap='Greens')  
fig.colorbar(s, ax=ax)

ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel("UMI Counts")
ax.set_ylabel("Genes Detected")

ax.set_xlim(1, 10**5)
ax.set_ylim(1, 10**4)

plt.show()<jupyter_output><empty_output><jupyter_text>Mitochondrial content per cell<jupyter_code>fig, ax = plt.subplots(figsize=(5,5))

x = nd(adata.obs["cell_counts"])
y = nd(adata.obs["percent_mito"])

ax.scatter(x, y, color="green", alpha=0.25)

ax.axhline(y=mito_criteria, linestyle="--", color="k")

ax.set_xlabel("UMI Counts")
ax.set_ylabel("Percent mito")

plt.show()<jupyter_output><empty_output><jupyter_text>Filter the anndata object<jupyter_code>adata.obs["pass_count_filter"] = adata.obs["cell_counts"] > cell_threshold
adata.obs["pass_mito_filter"] = adata.obs.percent_mito < mito_criteria
adata.var["pass_gene_filter"] = adata.var["min_cells"] >= gene_threshold
cell_mask = np.logical_and(adata.obs["pass_count_filter"].values, adata.obs["pass_mito_filter"].values)

gene_mask = adata.var["pass_gene_filter"].values
print("Current Shape: {:,} cells x {:,} genes".format(adata.shape[0], adata.shape[1]))
print("    New shape: {:,} cells x {:,} genes".format(cell_mask.sum(), gene_mask.sum()))
data = adata[cell_mask, gene_mask]
data
data.obs.SRR.value_counts()<jupyter_output><empty_output><jupyter_text>Process matrix for clustering, visualization, and de<jupyter_code>data.layers["raw"] = data.X
scale_num = 1000000
data.layers["norm"] = normalize(data.X, norm="l1", axis=1)*scale_num
data.layers["log1p"] = np.log1p(data.layers["norm"])
data.uns = OrderedDict([("log1p", {"base":None})])
data.X = data.layers["log1p"]
sc.pp.highly_variable_genes(data, n_top_genes=n_top_genes, flavor=flavor, n_bins=n_bins)
hvg_mask = data.var.highly_variable.values
%%time
# Scaling operation cannot be performed on the sparse matrix, so we make a dense matrix
mat = data.layers["log1p"].todense()
data.layers["scale"] = scale(mat, axis=0, with_mean=True, with_std=True, copy=True)
data.X = data.layers["scale"]

del mat<jupyter_output>CPU times: user 7.43 s, sys: 3.94 s, total: 11.4 s
Wall time: 11.4 s<jupyter_text>Generating visualizations  PCA on highly variable genes<jupyter_code>%%time
# PCA
X = data.X[:,hvg_mask]

tsvd = TruncatedSVD(n_components=num_PCA)
data.obsm["X_pca"] = tsvd.fit_transform(X)<jupyter_output>CPU times: user 45.9 s, sys: 1.94 s, total: 47.8 s
Wall time: 3.68 s<jupyter_text>Clustering<jupyter_code>%%time
sc.pp.neighbors(data, n_neighbors=n_neighbors, n_pcs=num_PCA, random_state=state)
%%time
sc.tl.leiden(data, random_state=state)
sc.tl.louvain(data, random_state=state)<jupyter_output><empty_output><jupyter_text>Dimensionality reduction<jupyter_code>%%time
# TSNE
X = data.obsm["X_pca"]

tsne = TSNE(n_components=num_TSNE, metric=metric, random_state=state)
data.obsm["X_pca_tsne"] = tsne.fit_transform(X)
# # NCA
# X = data.X
# y = data.obs.leiden.values # you could use louvain here also
# 
# nca = NeighborhoodComponentsAnalysis(n_components=num_NCA,random_state=state)
# data.obsm["X_nca"] = nca.fit_transform(X, y)
# # TSNE the NCA
# X = data.obsm["X_nca"]
# tsne = TSNE(n_components=num_TSNE, metric=metric, random_state=state)
# data.obsm["X_nca_tsne"] = tsne.fit_transform(X)
fig, ax = plt.subplots(figsize=(7,7))

x = data.obsm["X_pca_tsne"][:,0]
y = data.obsm["X_pca_tsne"][:,1]
c = data.obs["leiden"].astype(int)

ax.scatter(x, y, c = c, cmap='tab20c')

ax.set_axis_off()

plt.tight_layout()
plt.show()
del data.var["gene_name"]
%%time
data.write("write/data.h5ad")
data
def get_xyc(data, gene):
    c = nd(data.layers["raw"].todense()[:,data.var.index.values==gene])
    x = data.obsm["X_pca_tsne"][:,0]
    y = data.obsm["X_pca_tsne"][:,1]
    return (x ,y, c)
def fix_xyc(x,y,c, mask):
    c_f = c[mask]
    argsort = np.argsort(c_f)
    c_f = c_f[argsort]
    x_f = x[mask][argsort]
    y_f = y[mask][argsort]
    
    return (x_f, y_f, c_f)
young_m = nd(data.obs.age.values==3)
old_m =   nd(data.obs.age.values==24)
gene = "Tmprss2"
tmprss2 = get_xyc(data, gene)

gene = "Ace2"
ace2 = get_xyc(data, gene)<jupyter_output><empty_output><jupyter_text>Pass stacked C to Reds and Bluesmask the result by old v young and that is the vector c<jupyter_code>def plot_old_young(x,y,c, ax):
    cmask = c>0
    c = c/c.max()
    x = x
    y = y
    
    young_mask = young_m
    old_mask = old_m
    
    r = cm.Reds(c)
    b = cm.Blues(c)
    
    cc = np.zeros((len(c), 4))
    
    cc[young_mask] = r[young_mask,:]
    cc[old_mask] = b[old_mask,:]
    
    argsort = np.argsort(c)
    
    xx = x[argsort]
    yy = y[argsort]
    cc = cc[argsort]
    
    ax.scatter(xx, yy,c=cc, alpha=1)
    return ax

fig, ax = plt.subplots(figsize=(7,14), nrows=2)


gene = "Ace2"
x,y,c = ace2
plot_old_young(x,y,c,ax[0])
ax[0].set_axis_off()
ax[0].set_title("{} expression".format(gene))

gene = "Tmprss2"
x,y,c = tmprss2
plot_old_young(x,y,c,ax[1])
ax[1].set_axis_off()
ax[1].set_title("{} expression".format(gene))


red = mpatches.Patch(color='red', label='3 month', alpha=0.5)
blue = mpatches.Patch(color='blue', label='24 month', alpha=0.5)
ax[0].legend(bbox_to_anchor=(0, 0.5), loc="center right", handles=[red, blue])

plt.savefig("./figs/pca_ace2_tmprss_old_young.png", bbox_inches='tight',dpi=300)

#plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>DE<jupyter_code>sc.tl.rank_genes_groups(data, 'leiden', method='t-test', corr_method="bonferroni")
genes = pd.DataFrame(data.uns['rank_genes_groups']['names']).to_numpy()
pvals = pd.DataFrame(data.uns['rank_genes_groups']['pvals']).to_numpy()


unique = np.unique(data.obs.leiden.values.astype(int)).astype(str)
markers_gene = pd.DataFrame(index=unique, columns=["gene_name", "p_value"])

for un, u in enumerate(unique):
    g = genes[:,un]
    p = pvals[:,un]
    markers_gene.loc[u]["gene_name"]  = g.tolist()
    markers_gene.loc[u]["p_value"] = p.tolist()
    

markers_gene = markers_gene.apply(pd.Series.explode).reset_index()
markers_gene = markers_gene.rename(columns={"index":'leiden'})
markers_gene.to_csv("write/markers_gene.csv")
markers_gene
data<jupyter_output><empty_output><jupyter_text>Per SRR info<jupyter_code>reads = {"SRR8426358": 109256158,
"SRR8426359": 75353664,
"SRR8426360": 76904856,
"SRR8426361": 62538569,
"SRR8426362": 65131616,
"SRR8426363": 56363976,
"SRR8426364": 61732471,
"SRR8426365": 70225586,
"SRR8426366": 68168108,
"SRR8426367": 60893140,
"SRR8426368": 59394165,
"SRR8426369": 74088609,
"SRR8426370": 68834163,
"SRR8426371": 67565326,
"SRR8426372": 49144954}

pseudo = {"SRR8426358": 75.9,
"SRR8426359": 74.3,
"SRR8426360": 69.0,
"SRR8426361": 79.5,
"SRR8426362": 67.8,
"SRR8426363": 72.7,
"SRR8426364": 57.4,
"SRR8426365": 71.8,
"SRR8426366": 73.8,
"SRR8426367": 71.0,
"SRR8426368": 65.8,
"SRR8426369": 67.3,
"SRR8426370": 66.7,
"SRR8426371": 73.2,
"SRR8426372": 78.0}
df = pd.DataFrame(data.obs.SRR.value_counts().values, index = data.obs.SRR.value_counts().index, columns=["vc"])
df["reads"] = df.index.map(reads).astype(int)
df["pseudo"] = df.index.map(pseudo).astype(float)
df["p_mito_mean"] = df.index.map(data.obs.groupby("SRR")["percent_mito"].mean())
df.sort_values("vc")
df.reads.sum()<jupyter_output><empty_output><jupyter_text>Y chromosome hits (delete)<jupyter_code>"Uba1y", "Kdm5d", "Eif2s3y"
gene = "Kdm5d"
arr = nd(adata.X[:,gene == adata.var.index.values].todense())

df = adata.obs

df[gene] = arr

df.groupby("SRR")[gene].sum()
gene = "Eif2s3y"
arr = nd(adata.X[:,gene == adata.var.index.values].todense())

df = adata.obs

df[gene] = arr

df.groupby("SRR")[gene].sum()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17399476953252696, 'fraction_numerical': 0.057780320366132724, 'mean_word_length': 5.603885590933621, 'pattern_counts': {'":': 47, '<': 58, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '555746', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 5750, 'n_tokens_neox': 5211, 'n_words': 901, 'path': 'notebooks/make_matrix_old_v_young.ipynb', 'repo_name': 'pachterlab/BP_2020', 'size': '12236'}