diff --git a/get_training_test_looms.py b/get_training_test_looms.py index a099fce..0bb817e 100644 --- a/get_training_test_looms.py +++ b/get_training_test_looms.py @@ -1,44 +1,74 @@ #! /usr/bin/python -import loompy -import numpy as np -import argparse +import argparse, loompy, sys, numpy as np from numpy.random import choice - def parse_user_input(): parser = argparse.ArgumentParser() - parser.add_argument('-l','--loom',required=True,help='Path to loom.') - parser.add_argument('-a','--attribute',required=True,help='Column attribute for sampling.') - parser.add_argument('-p','--prefix',required=True,help='Output prefix.') - parser.add_argument('-n','--n-test-cells',type=int,required=True,help='Number of test cells for each sample.') + parser.add_argument('-l', '--loom', required=True, + help='Path to loom file.') + parser.add_argument('-a', '--attribute', required=True, + help='Column attribute for sampling.') + parser.add_argument('-p', '--prefix', required=True, + help='Output prefix.') + parser.add_argument('-n', '--n-test-cells', type=int, required=False, + help='Number of test cells for each sample.', + default=None) + parser.add_argument('-c', '--n-test-cells-pct', type=float, required=False, + help='Percent of test cells for each sample.', + default=None) + parser.add_argument('-s', '--seed', type=int, required=False, + help='Seed. Default is 0.', default=0) return parser parser = parse_user_input() ui = parser.parse_args() +np.random.seed(ui.seed) -data = loompy.connect(ui.loom,validate=False) -attr = data.ca[ui.attribute] - -train_ix=[] -test_ix=[] -for a in set(attr): - ix = np.where(attr==a)[0] - test = choice(ix,ui.n_test_cells,replace=False) - test_ix.extend(list(test)) - train_ix.extend([i for i in ix if i not in test]) -train_ix=np.array(train_ix) -test_ix = np.array(test_ix) +if ui.n_test_cells is not None and ui.n_test_cells_pct is not None: + parser.error("Options --n-test-cells and --n-test-cells-pct cannot be " + "supplied together. Please choose only one.") +elif ui.n_test_cells is None and ui.n_test_cells_pct is None: + parser.error("Please specify one of the options --n-test-cells or " + "--n-test-cells-pct.") -print(len(attr)) -print(len(train_ix)) -print(len(test_ix)) +try: + data = loompy.connect(ui.loom, validate=False) + attr = data.ca[ui.attribute] -train_output = ui.prefix+'.train.loom' -with loompy.new(train_output) as dsout: - for (ix,selection,view) in data.scan(items=train_ix,axis=1): - dsout.add_columns(view.layers,col_attrs=view.ca,row_attrs=view.ra) + train_ix=[] + test_ix=[] + for a in set(attr): + ix = np.where(attr == a)[0] + if ui.n_test_cells_pct is not None and ui.n_test_cells is None: + ix_n = len(ix) + sample_n = int(ix_n * ui.n_test_cells_pct) + if sample_n < 1: + sample_n = 1 + test = choice(ix, sample_n, replace=False) + elif ui.n_test_cells is not None and ui.n_test_cells_pct is None: + test = choice(ix, ui.n_test_cells, replace=False) + else: + raise ValueError("Sampling failed. Check your --n-test-cells or" + "--n-test-cells-pct.") + test_ix.extend(test) + train_ix.extend([i for i in ix if i not in test]) + train_ix = np.array(train_ix) + test_ix = np.array(test_ix) + + print(f'Original number of cells: {len(attr)}') + print(f'Training number of cells: {len(train_ix)}') + print(f'Testing number of cells: {len(test_ix)}') + + train_output = ui.prefix+'.train.loom' + with loompy.new(train_output) as dsout: + for (ix,selection, view) in data.scan(items=train_ix, axis=1): + dsout.add_columns(view.layers, col_attrs=view.ca, row_attrs=view.ra) + + test_output = ui.prefix+'.test.loom' + with loompy.new(test_output) as dsout: + for (ix,selection, view) in data.scan(items=test_ix, axis=1): + dsout.add_columns(view.layers, col_attrs=view.ca, row_attrs=view.ra) +except Exception as e: + print("Error:", e) + sys.exit(1) -test_output = ui.prefix+'.test.loom' -with loompy.new(test_output) as dsout: - for (ix,selection,view) in data.scan(items=test_ix,axis=1): - dsout.add_columns(view.layers,col_attrs=view.ca,row_attrs=view.ra) diff --git a/scHPF_consensus.py b/scHPF_consensus.py index 838a972..2ce23b5 100644 --- a/scHPF_consensus.py +++ b/scHPF_consensus.py @@ -1,54 +1,70 @@ #! /usr/bin/python import argparse import subprocess -from glob import glob -from copy import deepcopy +import warnings import joblib +import schpf +import pandas as pd import numpy as np import igraph as ig -from sklearn import neighbors,metrics +import matplotlib as mpl +from sklearn import neighbors, metrics from scipy.sparse import coo_matrix from scipy.io import mmread -import pandas as pd -import schpf -import matplotlib as mpl +from glob import glob +from copy import deepcopy from matplotlib import pyplot as plt -plt.rcParams['pdf.fonttype'] = 42 -plt.rcParams['ps.fonttype'] =42 def parse_user_input(): parser = argparse.ArgumentParser() - parser.add_argument('-trn','--infile-train',required=True,help='Path to training data output by scHPF prep.') - parser.add_argument('-tst','--infile-test',required=True,help='Path to test data output by scHPF prep-like.') - parser.add_argument('-o','--outdir',required=True,help='Path to output directory.') - parser.add_argument('-p','--prefix',required=True,help='Prefix for output files.') - parser.add_argument('-g','--gene-infile',required=True,help='Path to gene name file.') - parser.add_argument('-k','--k-values',required=True,type=int,nargs='+',help='Values of k to test.') - parser.add_argument('-t','--trials',required=True,type=int,help='Number of trials for each value of k.') - parser.add_argument('-n','--n-models',required=True,type=int,help='Number of models to consider for each value of k.') - parser.add_argument('-m','--min-cluster-size',required=True,type=int,help='Minimum number of factors required for keeping a cluster of factors.') - parser.add_argument('-j','--jobs',required=False,type=int,help='Maximum number of jobs to run in parallel.') + parser.add_argument('-trn', '--infile-train', required=True, + help='Path to training data output by scHPF prep.') + parser.add_argument('-tst', '--infile-test', required=True, + help='Path to test data output by scHPF prep-like.') + parser.add_argument('-o', '--outdir', required=True, + help='Path to output directory.') + parser.add_argument('-p', '--prefix', required=True, + help='Prefix for output files.') + parser.add_argument('-g', '--gene-infile', required=True, + help='Path to gene name file.') + parser.add_argument('-k', '--k-values', required=True, type=int, nargs='+', + help='Values of k to test.') + parser.add_argument('-t', '--trials', required=True, type=int, + help='Number of trials for each value of k.') + parser.add_argument('-n', '--n-models', required=True, type=int, + help='Number of models to consider for each value of k.') + parser.add_argument('-m', '--min-cluster-size', required=True, type=int, + help='Minimum number of factors required for keeping a cluster of factors.') + parser.add_argument('-j', '--jobs', required=False, type=int, + help='Maximum number of jobs to run in parallel.') + parser.add_argument('-r', '--re-write', required=False, default=False, + action='store_true', + help='Re-write .joblib model file. Default False.') + parser.add_argument('-f', '--n-factors', required=False, type=int, + help='Number of factors for consensus scHPF. Default ' + 'is igraph VertexDendrogram.optimal_count().') return parser # function for scHPF refit based on median parameters of factor clusters -def refit_local_params(X, global_params, nfactors, bp, dp, a=0.3, c=0.3, project_kw={}): +def refit_local_params(X, global_params, nfactors, bp, dp, a=0.3, c=0.3, + project_kw={}): """ """ project_defaults = dict(verbose=True, max_iter=50, check_freq=5) eta_shp, eta_rte, beta_shp, beta_rte = global_params - + # make a model eta = schpf.HPF_Gamma(np.ravel(eta_shp), np.ravel(eta_rte)) beta = schpf.HPF_Gamma(beta_shp.T.values, beta_rte.T.values) - model = schpf.scHPF(nfactors, eta=eta, beta=beta, bp=bp, dp=dp, a=a, c=c) - + model = schpf.scHPF(nfactors, eta=eta, beta=beta, bp=bp, dp=dp, a=a, c=c) + # setup projection kwarg - for k,v in project_defaults.items(): + for k, v in project_defaults.items(): if k not in project_kw.keys(): project_kw[k] = v loss = model.project(X, replace=True, **project_kw) model.loss = loss - + return model # utility function for extracting model parameters into pandas @@ -56,24 +72,30 @@ def get_param_dfs(model): eta_shp = pd.Series(np.ravel(model.eta.vi_shape), name=model.name) eta_rte = pd.Series(np.ravel(model.eta.vi_rate), name=model.name) beta_shp = pd.DataFrame(model.beta.vi_shape.T) - beta_shp.index = model.name + ':' + (beta_shp.index + 1).astype(str) + beta_shp.index = (model.name + '_' + str(model.nfactors) + ':' + + (beta_shp.index + 1).astype(str)) beta_rte = pd.DataFrame(model.beta.vi_rate.T, index=beta_shp.index) return eta_shp, eta_rte, beta_shp, beta_rte -# function for converting model parameters into pandas dataframe +# function for converting model parameters into pandas dataframe def get_spectra(models): - eta_shp, eta_rte, beta_shp, beta_rte = zip(*[get_param_dfs(m) for m in models]) - return pd.concat(eta_shp, axis=1).T, pd.concat(eta_rte,axis=1).T, pd.concat(beta_shp), pd.concat(beta_rte) + eta_shp, eta_rte, beta_shp, beta_rte = zip( + *[get_param_dfs(m) for m in models]) + return pd.concat(eta_shp, axis=1).T, pd.concat(eta_rte, + axis=1).T, pd.concat( + beta_shp), pd.concat(beta_rte) # function for extracting gene scores from model object into pandas dataframe def get_genescore_spectra(models): gene_scores = [] for m in models: gs = pd.DataFrame(m.gene_score().T) - gs.index = m.name + ':' + (gs.index + 1).astype(str) + gs.index = (m.name + '_' + str(m.nfactors) + ':' + + (gs.index + 1).astype(str)) gene_scores.append(gs) return pd.concat(gene_scores) +# parse user input parser = parse_user_input() ui = parser.parse_args() @@ -82,126 +104,215 @@ def get_genescore_spectra(models): prefix = ui.prefix trials = ui.trials genes_infile = ui.gene_infile +n_factors = ui.n_factors -if not ui.jobs: # if not restriction on number of parallel jobs, run an scHPF training job for each value of k - procs=[] +# run scHPF training for each value of k +if not ui.jobs: # if not restriction on number of parallel jobs, run an scHPF training job for each value of k + procs = [] for k in ui.k_values: + outfile = outdir + '/' + prefix + (f'.scHPF_K{k}_*' + f'{trials}trials.joblib') + outfile = glob(outfile) + if len(outfile) > 0 and not re_write: + print(f'Skipping existing model for k = {k} at {outfile[0]}.') + continue + print(f'Running scHPF training for k = {k}') cmd = 'scHPF train -i %(infile)s -o %(outdir)s -p %(prefix)s -k %(k)d -t %(trials)d --save-all' % vars() - p = subprocess.Popen(cmd,shell=True) + p = subprocess.Popen(cmd, shell=True) procs.append(p) - p_exit = [p.wait() for p in procs] + p_exit = [p.wait() for p in procs] else: # otherwise, run only ui.jobs k-values at-a-time - st=0 - sp=ui.jobs + st = 0 + sp = ui.jobs while st <= len(ui.k_values): - procs=[] + procs = [] for k in ui.k_values[st:sp]: + outfile = outdir + '/' + prefix + (f'.scHPF_K{k}_*' + f'{trials}trials.joblib') + outfile = glob(outfile) + if len(outfile) > 0 and not ui.re_write: + print(f'Skipping existing model for k = {k} at {outfile[0]}.') + continue + print(f'Running scHPF training for k = {k}') cmd = 'scHPF train -i %(infile)s -o %(outdir)s -p %(prefix)s -k %(k)d -t %(trials)d --save-all' % vars() - p = subprocess.Popen(cmd,shell=True) + p = subprocess.Popen(cmd, shell=True) procs.append(p) p_exit = [p.wait() for p in procs] - st+=ui.jobs - sp+=ui.jobs + st += ui.jobs + sp += ui.jobs # get the model objects for the top ui.n_models models for each value of k -models_str = ui.outdir+'/*scHPF*.joblib' +models_str = ui.outdir + '/*scHPF*.joblib' model_infiles = sorted(glob(models_str)) -top_model_infiles = [model_infile for model_infile in model_infiles if ('reject' not in model_infile or int(model_infile.split('reject')[1][0])0: - gscores = np.concatenate((gscores,model.gene_score()),axis=1) +top_model_infiles = [model_infile for model_infile in model_infiles if ( + 'reject' not in model_infile or int(model_infile.split('reject')[1][ + 0]) < ui.n_models)] +top_model_names = [] +for model_infile in top_model_infiles: + model_name = model_infile.split('/')[-1].rsplit('.', -1)[0] + if 'reject' in model_infile: + rejectn = model_infile.split('/')[-1].split('.')[1].split('_')[-1] + top_model_names.append(f'{model_name}_{rejectn}') else: - gscores=model.gene_score() + top_model_names.append(model_name) +top_model_Ks = [int(model_infile.split('scHPF_K')[1].split('_')[0]) for + model_infile in top_model_infiles] +top_model_dfs = pd.DataFrame(list(zip(top_model_names, top_model_Ks, + top_model_infiles)), + columns=['name', 'K', 'model_file']) +top_models = [joblib.load(model_infile) for model_infile in top_model_infiles] +for model, name in zip(top_models, top_model_names): + model.name = name + +# concatenate factors across models and select highly variable genes gscores = get_genescore_spectra(top_models) -gscore_cvs = (gscores.std()/gscores.mean()) +print(f'\nClustering {gscores.shape[0]} factors across {len(top_models)} ' + f'models with {gscores.shape[1]} genes') +gscore_cvs = (gscores.std() / gscores.mean()) +print(f'Gene coefficient of variation distribution\n' + f'{gscore_cvs.describe()}\n' + f'Selecting top 1000 HVGs for clustering') top_gene_ixs = gscore_cvs.nlargest(1000).index.values gscores = gscores[top_gene_ixs] -print(gscores.shape) - -eta_shp,eta_rte,beta_shp,beta_rte = get_spectra(top_models) -eta_ex = eta_shp/eta_rte -beta_ex = beta_shp/beta_rte -# convert gene score matrix into distance matrix and knn graph -factor_dists = 1.-pd.DataFrame(gscores,index=gscores.index).T.corr().values -n_neighbors = max(5,int(0.25*len(top_model_infiles))) # heuristic for k in knn graph -adj_binary = neighbors.kneighbors_graph(factor_dists,n_neighbors,metric='precomputed') +# convert gene score matrix into knn graph +n_neighbors = max(5, int(0.25 * len( + top_model_infiles))) # heuristic for k in knn graph +adj_binary = neighbors.kneighbors_graph(pd.DataFrame(gscores, + index=gscores.index), + n_neighbors, metric='euclidean') adj = np.zeros(adj_binary.shape) -for i,j in np.stack(adj_binary.nonzero()).T: - adj[i,j] = metrics.jaccard_score(adj_binary[i,:].A[0], adj_binary[j,:].A[0]) -adj=coo_matrix(adj) +for i, j in np.stack(adj_binary.nonzero()).T: + adj[i, j] = metrics.jaccard_score(adj_binary[i, :].A[0], + adj_binary[j, :].A[0]) -vcount = max(adj.shape) -sources,targets = adj.nonzero() -edgelist = list(zip(sources.tolist(),targets.tolist())) +adj = coo_matrix(adj) +sources, targets = adj.nonzero() +edgelist = list(zip(sources.tolist(), targets.tolist())) # perform walktrap clustering on knn graph -knn = ig.Graph(edges=edgelist,directed=False) -knn.vs['label']=gscores.index -knn.es['width']=adj.data -knn.es['weight']=adj.data -cluster_result = knn.community_walktrap(weights=adj.data,steps=4) -nclusters = cluster_result.optimal_count+2 # heuristic for numbr of clusters -cluster_labels = pd.Series(cluster_result.as_clustering(nclusters).membership,index=gscores.index) - -# compute cluster median parameters to initialization scHPF refit +knn = ig.Graph(edges=edgelist, directed=False) +knn.vs['label'] = gscores.index +knn.es['width'] = adj.data +knn.es['weight'] = adj.data +cluster_result = knn.community_walktrap(weights=adj.data, steps=4) + +# Select number of clusters based on modularity +max_k = len(cluster_result.merges) +min_k = 1 # find the minimum valid number of clusters +for i in range(1, max_k + 1): + try: + cluster_result.as_clustering(i) + min_k = i + break + except ig._igraph.InternalError: + continue + +if n_factors is None: # Check for user-supplied number of factors + nclusters = cluster_result.optimal_count +elif 0 < n_factors <= max_k: + nclusters = n_factors +else: + warnings.warn(f'"--n-factors" must be between 1 and {max_k}. ' + f'Using default value.') + nclusters = cluster_result.optimal_count + +if nclusters < min_k: + warnings.warn(f'Optimal number of clusters is {nclusters}, but minimum ' + f'number of clusters where modularity can be calculated is ' + f'{min_k}. Using minimum valid number of clusters.') + nclusters = min_k + +print(f'Number of clusters: {nclusters}') +cluster_labels = pd.Series(cluster_result.as_clustering(nclusters).membership, + index=gscores.index) +min_cluster_size = ui.min_cluster_size +keep = np.where(cluster_labels.value_counts(sort=False) >= min_cluster_size)[0] +cluster_labels = cluster_labels.loc[cluster_labels.isin(keep)] +print(f'Number of clusters with n_factors > {min_cluster_size}: {nclusters}') + +# compute modularity of Walktrap clustering across all posible clusters +ks = np.arange(min_k, max_k + 1) +print(f'\nCalculating modularity for {min_k} - {max_k + 1} k clusters') +modularity = [] +for i in ks: + try: + modularity.append(cluster_result.as_clustering(i).modularity) + except ig._igraph.InternalError as e: + print(f"Error calculating modularity at k {i}: {e}") + +modularity = pd.DataFrame({'k': ks, 'modularity': modularity}) +modularity['diff_modularity'] = modularity['modularity'].diff().fillna(0) +modularity.to_csv(outdir + '/' + prefix + '.walktrap.tsv', sep='\t', + index=False) + +# plot cluster modularity +pdf_outfile = outdir + '/' + prefix + '.walktrap.pdf' +step = min(len(modularity)//10, len(modularity)//3) +plt.figure(figsize=(9, 6)) +plt.rcParams.update({'font.size': 14}) +plt.plot(modularity['k'], modularity['modularity']) +plt.xlim(modularity['k'].min(), modularity['k'].max()) +plt.xticks(modularity['k'].iloc[::step]) +plt.axvline(nclusters, c='r') +plt.text(nclusters + 0.5, plt.ylim()[1] * 0.6, f'Optimal k = {nclusters}', + color='r') +plt.title(f'Modularity of Walktrap across {gscores.shape[0]}\n' + f'factors from n_models ({top_model_dfs["name"].nunique()}) * ' + f'k ({top_model_dfs["K"].nunique()}) models') +plt.xlabel('Number of Clusters (k)') +plt.ylabel('Modularity') +plt.grid(True) +plt.tight_layout() +plt.savefig(pdf_outfile) +plt.close() + +# compute cluster median parameters to initialize scHPF refit +eta_shp, eta_rte, beta_shp, beta_rte = get_spectra(top_models) +eta_ex = eta_shp / eta_rte +beta_ex = beta_shp / beta_rte eta_shp_med = eta_shp.median().values eta_rte_med = eta_rte.median().values beta_shp_med = beta_shp.groupby(cluster_labels).median() beta_rte_med = beta_rte.groupby(cluster_labels).median() -eta_ex_med = eta_shp_med/eta_rte_med -beta_ex_med = beta_shp_med/beta_rte_med -min_cluster_size = ui.min_cluster_size -keep = np.where(cluster_labels.value_counts(sort=False)>=min_cluster_size)[0] -cluster_labels=cluster_labels.loc[cluster_labels.isin(keep)] - -# make cluster modularity plot -pdf_outfile = outdir+'/'+prefix+'.walktrap.pdf' -optimal_count = cluster_result.optimal_count -x=np.arange(-5,10) -modularity = [cluster_result.as_clustering(optimal_count+i).modularity for i in x] -fig,ax=plt.subplots() -ax.plot(x+optimal_count,modularity) -ax.axvline(nclusters,c='r') -ax.set_xlabel('Number of Clusters') -ax.set_ylabel('Modularity') -fig.savefig(pdf_outfile) - -outfile1 = outdir+'/'+prefix+'.consensus.joblib' +eta_ex_med = eta_shp_med / eta_rte_med +beta_ex_med = beta_shp_med / beta_rte_med + +# initializing consensus scHPF by refitting with cluster median parameters +outfile1 = outdir + '/' + prefix + '.consensus.joblib' np.random.seed(0) nfactors = cluster_labels.nunique() -print(nfactors) -a=0.3 -c=0.3 +assert nfactors > 0, 'No valid clusters found.' +print(f'\nComputing consensus scHPF with {nfactors} factors') +a = 0.3 +c = 0.3 for model in top_models: if model.nfactors == nfactors: - a=model.a - c=model.c + a = model.a + c = model.c break -# compute initial consensus scHPF by refitting with initialization from cluster median parameters matrix = mmread(infile) -consensus1 = refit_local_params(matrix, (eta_shp_med, eta_rte_med,beta_shp_med.iloc[keep],beta_rte_med.iloc[keep]), - nfactors,bp=top_models[0].bp,dp=top_models[0].dp,a=a,c=c,project_kw={'max_iter':1}) -joblib.dump(consensus1,outfile1) +consensus1 = refit_local_params(matrix, (eta_shp_med, eta_rte_med, + beta_shp_med.iloc[keep], + beta_rte_med.iloc[keep]), nfactors, + bp=top_models[0].bp, dp=top_models[0].dp, + a=a, c=c, project_kw={'max_iter': 1}) +joblib.dump(consensus1, outfile1) # compare consensus scHPF to randomly initialized model with same k using test data -outfile2 = outdir+'/'+prefix+'.consensus.final.joblib' +outfile2 = outdir + '/' + prefix + '.consensus.final.joblib' test_matrix = mmread(ui.infile_test) -test_loss = schpf.loss.projection_loss_function(schpf.loss.mean_negative_pois_llh,test_matrix,consensus1.nfactors,proj_kwargs={'reinit':False, 'verbose':False}) +test_loss = schpf.loss.projection_loss_function( + schpf.loss.mean_negative_pois_llh, test_matrix, consensus1.nfactors, + proj_kwargs={'reinit': False, 'verbose': False}) # update consensus scHPF training until convergence starting from initial consensus scHFP consensus2 = deepcopy(consensus1) np.random.seed(0) -consensus2.fit(matrix, loss_function=test_loss, reinit=False, verbose=True, max_iter=150) -joblib.dump(consensus2,outfile2) - +consensus2.fit(matrix, loss_function=test_loss, reinit=False, verbose=True, + max_iter=150) +joblib.dump(consensus2, outfile2) +print('Consensus scHPF has finished.') diff --git a/subsample_merged_loom.py b/subsample_merged_loom.py index cf721bc..a411693 100644 --- a/subsample_merged_loom.py +++ b/subsample_merged_loom.py @@ -1,38 +1,66 @@ -#! /usr/bin/python +#!/usr/bin/python import loompy import argparse import numpy as np - def parse_user_input(): - """ - Get and parse user input. - """ - parser = argparse.ArgumentParser() - parser.add_argument('-o','--output-loom',required=True,help='Path to output loom file.') - parser.add_argument('-i','--input-loom',required=True,help='Path to input loom file.') - parser.add_argument('-n','--number-of-cells',required=False,type=int,help='Subsample to this fixed number of cells.') - parser.add_argument('-c','--col-attr',required=False,help='Identify cells for sub-sampling based on this column attribute.') - parser.add_argument('-v','--col-val',nargs='+',required=False,help='Value of column attribute for cells to send for sub-sampling.') - parser.add_argument('-sc','--subsample-col-attr',required=True,help='Column attribute to subset cells for sub-sampling.') - return parser + """ + Get and parse user input. + """ + parser = argparse.ArgumentParser() + parser.add_argument('-o','--output-loom',required=True,help='Path to output loom file.') + parser.add_argument('-i','--input-loom',required=True,help='Path to input loom file.') + parser.add_argument('-n','--number-of-cells',required=False,type=int,help='Subsample to this fixed number of cells.') + parser.add_argument('-p','--pct-of-cells',required=False,type=float,help='Subsample this percentage of cells.') + parser.add_argument('-c','--col-attr',required=False,help='Identify cells for sub-sampling based on this column attribute.') + parser.add_argument('-v','--col-val',nargs='+',required=False,help='Value of column attribute for cells to send for sub-sampling.') + parser.add_argument('-sc','--subsample-col-attr',required=True,help='Column attribute to subset cells for sub-sampling.') + return parser parser = parse_user_input() ui = parser.parse_args() +np.random.seed(0) + with loompy.new(ui.output_loom) as dsout: # Create a new, empty, loom file - with loompy.connect(ui.input_loom,validate=False) as ds: - if ui.col_attr: - cells = np.where(np.isin(ds.ca[ui.col_attr],ui.col_val))[0] - else: - cells = np.arange(ds.shape[1]) - scts = {sample:float(sum(ds.ca[ui.subsample_col_attr][cells]==sample)) for sample in set(ds.ca[ui.subsample_col_attr][cells])} - if not ui.number_of_cells: - mnct = np.min(list(scts.values())) - else: - mnct = ui.number_of_cells - sfracs = {sample:mnct/scts[sample] for sample in set(ds.ca[ui.subsample_col_attr][cells])} - rnd = np.random.rand(ds.shape[1]) - keep = np.array([cell for cell in cells if rnd[cell] < sfracs[ds.ca[ui.subsample_col_attr][cell]]]) - for (ix, selection, view) in ds.scan(items=keep, axis=1, key="Accession"): - dsout.add_columns(view.layers,col_attrs=view.ca,row_attrs=view.ra) + with loompy.connect(ui.input_loom,validate=False) as ds: + print(f'Detected {ds.shape[0]} molecules across {ds.shape[1]} cells.') + if ui.col_attr: + cells = np.where(np.isin(ds.ca[ui.col_attr], ui.col_val))[0] + else: + cells = np.arange(ds.shape[1]) + + scts = {sample: float(sum(ds.ca[ui.subsample_col_attr][cells] == sample)) for sample in set(ds.ca[ui.subsample_col_attr][cells])} + scts_values = np.array(list(scts.values())) + quantiles = np.quantile(scts_values, [0.25, 0.5, 0.75]) + + print(f'Distribution of cells per {ui.subsample_col_attr}') + print(f'count = {len(scts_values)}') + print(f'min = {scts_values.min()}') + print(f'mean = {scts_values.mean():.2f}') + print(f'median = {quantiles[1]:.2f}') + print(f'25th percentile = {quantiles[0]:.2f}') + print(f'50th percentile (median) = {quantiles[1]:.2f}') + print(f'75th percentile = {quantiles[2]:.2f}') + print(f'max = {scts_values.max()}') + + if ui.number_of_cells: + mnct = ui.number_of_cells + elif ui.pct_of_cells and 0 < ui.pct_of_cells < 1: + mnct = ui.pct_of_cells + else: + mnct = scts_values.min() + + if not ui.pct_of_cells: + sfracs = {sample: mnct / scts[sample] for sample in scts} + else: + sfracs = {sample: mnct for sample in scts} + + rnd = np.random.rand(ds.shape[1]) + keep = np.array([cell for cell in cells if rnd[cell] < sfracs[ds.ca[ui.subsample_col_attr][cell]]]) + + for (ix, selection, view) in ds.scan(items=keep, axis=1, key="Accession"): + dsout.add_columns(view.layers, col_attrs=view.ca, row_attrs=view.ra) + + print(f'Subsampled down to {len(keep)} cells.') +