<jupyter_start><jupyter_text>Mouse celseq1 organoid data were collected from GSE62270 and processed with kallisto buskallisto bus -i index/Mus_musculus.GRCm38.cdna.all.fa.kidx -x CELSeq -t 8 -o kallisto_Organoid_out/L001/ \ celseq1/Organoid1_GSE62270_L001_R1.fastq celseq1/Organoid1_GSE62270_L001_R2.fastqbustools sort -t 8 -o kallisto_Organoid_out/L001/output.sort.bus kallisto_Organoid_out/L001/output.busbustools dump -o kallisto_Organoid_out/L001/output.sort.txt kallisto_Organoid_out/L001/output.sort.busChange L001 to L00X, X in {1,7} to process data from all 96 wells.<jupyter_code>import matplotlib
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline  
import sys, collections<jupyter_output><empty_output><jupyter_text>Load celseq1 kallisto bus results (L001 to L007) to generate gene count tables. Change L007 to L00X x in {1, 7} to process the bus results.<jupyter_code>#load mouse transcript

tr2g = {}
trlist = []
with open('/home/fgao/Data_single_cell/kallisto_benchmark/index/tr2g_mouse') as f:
    for line in f:
        l = line.split()
        tr2g[l[0]] = l[1]
        trlist.append(l[0])

genes = list(set(tr2g[t] for t in tr2g))

# load equivalence classes
ecs = {}
with open('/home/fgao/Data_single_cell/kallisto_benchmark/kallisto_Organoid_out/L007/matrix.ec') as f:
    for line in f:
        l = line.split()
        ec = int(l[0])
        trs = [int(x) for x in l[1].split(',')]
        ecs[ec] = trs
        
def ec2g(ec):
    if ec in ecs:
        return list(set(tr2g[trlist[t]] for t in ecs[ec]))        
    else:
        return []
# load test dataset

cell_gene = collections.defaultdict(lambda: collections.defaultdict(float))
pbar=None
pumi=None
with open('/home/fgao/Data_single_cell/kallisto_benchmark/kallisto_Organoid_out/L007/output.sort.txt') as f:
    gs = set()
    gl_sum=[]
    for line in f:
        l = line.split()
        barcode,umi,ec,count = line.split()
        ec = int(ec)
        
        if barcode == pbar:
            # same barcode
            if umi == pumi:
                # same UMI, let's add new genelist of existing genelist
                gl = ec2g(ec)
                gl_sum = gl_sum + gl
            else:
                # new UMI, process the previous gene set
                gs = set(gl_sum)
                g_max = max(gl_sum, key=gl_sum.count)
                for i in set(gl_sum):
                    if (gl_sum.count(i) < gl_sum.count(g_max)):
                        gs.remove(i)
                for g in gs:
                    cell_gene[barcode][g] += 1.0/len(gs)
                # record new umi, reset gene set
                pumi = umi
                gl_sum = ec2g(ec)
        else:
            # bypass the first bus record
            if (gl_sum == []):
                pbar = barcode
                pumi = umi
                gl_sum = ec2g(ec)
            else:
                # work with previous gene list
                gs = set(gl_sum)
                g_max = max(gl_sum, key=gl_sum.count)
                for i in set(gl_sum):
                    if (gl_sum.count(i) < gl_sum.count(g_max)):
                        gs.remove(i)
                for g in gs:
                    cell_gene[pbar][g] += 1.0/len(gs)
            
                if sum(cell_gene[pbar][g] for g in cell_gene[pbar]) < 10:
                    del cell_gene[pbar]
                pbar = barcode
                pumi = umi         
                gl_sum = ec2g(ec)
    #remember the last gene
    gs = set(gl_sum)
    g_max = max(gl_sum, key=gl_sum.count)
    for i in set(gl_sum):
        if (gl_sum.count(i) < gl_sum.count(g_max)):
            gs.remove(i)
    for g in gs:
        cell_gene[pbar][g] += 1.0/len(gs)
        
    if sum(cell_gene[pbar][g] for g in cell_gene[pbar]) < 10:
        del cell_gene[pbar]


barcode_hist = collections.defaultdict(int)
for barcode in cell_gene:
    cg = cell_gene[barcode]
    s = len([cg[g] for g in cg])
    barcode_hist[barcode] += s
    
whitelist = set(x.strip() for x in open('/home/fgao/Data_single_cell/kallisto_benchmark/kallisto_Organoid_out/BC.txt'))
bcv = [x for b,x in barcode_hist.items() if x > 0 and b in whitelist]
_ = plt.hist(bcv,bins=100)
print(len(bcv))
outfile = '/home/fgao/Data_single_cell/kallisto_benchmark/kallisto_Organoid_out/L007/matrix.mtx'

gene_to_id = dict((g,i+1) for i,g in enumerate(genes))
barcodes_to_use = [b for b,x in barcode_hist.items() if x > 0 and b in whitelist]

num_entries = 0
for barcode in barcodes_to_use:
    num_entries += len([x for x in cell_gene[barcode].values() if round(x)>0])


with open(outfile, 'w') as of:
    of.write('%%MatrixMarket matrix coordinate real general\n%\n')
    #number of genes
    of.write("%d %d %d\n"%(len(genes), len(barcodes_to_use), num_entries))
    bcid = 0
    for barcode in barcodes_to_use:
        bcid += 1
        cg = cell_gene[barcode]
        gl = [(gene_to_id[g],round(cg[g])) for g in cg if round(cg[g]) > 0]
        gl.sort()
        for x in gl:
            of.write("%d %d %d\n"%(x[0],bcid,x[1]))
    
gene_names = {}
with open('/home/fgao/Data_single_cell/kallisto_benchmark/index/mart_export_mouse.txt') as f:
    f.readline()
    for line in f:
        g,t,gn = line.split()
        gene_names[g] = gn

id_to_genes = dict((i,g) for (g,i) in gene_to_id.items())
gl = []
for i in range(1,len(genes)+1):
    g = id_to_genes[i]
    gid = g[:g.find('.')]
    if gid in gene_names:
        gn = gene_names[gid]
    else:
        gn = ''
    gl.append((g,gn))

with open('/home/fgao/Data_single_cell/kallisto_benchmark/kallisto_Organoid_out/L007/genes.tsv','w') as of:
    for g,gn in gl:
        of.write("%s\t%s\n"%(g,gn))
        
with open('/home/fgao/Data_single_cell/kallisto_benchmark/kallisto_Organoid_out/L007/barcodes.tsv','w') as of:
    of.write('\n'.join(x + '' for x in barcodes_to_use))
    of.write('\n')
import numpy as np
import pandas as pd
import scanpy.api as sc

sc.settings.verbosity = 3  # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_versions()
results_file = '/home/fgao/Data_single_cell/kallisto_benchmark/kallisto_Organoid_out/kallisto_Organoid.h5ad'
sc.settings.set_figure_params(dpi=80)

adata1 = sc.read_10x_mtx('../kallisto_Organoid_out/L001/', var_names='gene_symbols', cache=False)
adata2 = sc.read_10x_mtx('../kallisto_Organoid_out/L002/', var_names='gene_symbols', cache=False)
adata3 = sc.read_10x_mtx('../kallisto_Organoid_out/L003/', var_names='gene_symbols', cache=False)
adata4 = sc.read_10x_mtx('../kallisto_Organoid_out/L004/', var_names='gene_symbols', cache=False)
adata5 = sc.read_10x_mtx('../kallisto_Organoid_out/L005/', var_names='gene_symbols', cache=False)
adata6 = sc.read_10x_mtx('../kallisto_Organoid_out/L006/', var_names='gene_symbols', cache=False)
adata7 = sc.read_10x_mtx('../kallisto_Organoid_out/L007/', var_names='gene_symbols', cache=False)
adata = adata1.concatenate(adata2, adata3, adata4, adata5, adata6, adata7, join='outer')
adata.write_csvs('../kallisto_Organoid_out/', skip_data=False)
adata.var_names_make_unique()
sc.pp.filter_cells(adata, min_genes=100)
sc.pp.filter_genes(adata, min_cells=2)
mito_genes = [name for name in adata.var_names if name.startswith('mt-')]
# for each cell compute fraction of counts in mito genes vs. all genes
# the `.A1` is only necessary as X is sparse to transform to a dense array after summing
adata.obs['percent_mito'] = np.sum(adata[:, mito_genes].X, axis=1) / np.sum(adata.X, axis=1)
# add the total counts per cell as observations-annotation to adata
adata.obs['n_counts'] = adata.X.sum(axis=1)
adata = adata[adata.obs['n_genes'] < 2500, :]
adata
adata = adata[adata.obs['percent_mito'] < 0.2, :]
adata
#adata.raw = adata
sc.pp.log1p(adata)
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, log=True)
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=30)
sc.tl.umap(adata)
sc.tl.louvain(adata, resolution=1)
sc.pl.umap(adata, color=['louvain'], size=25)
sc.tl.tsne(adata, n_pcs=30)
sc.pl.tsne(adata, color=['louvain'], size=25)
sc.tl.rank_genes_groups(adata, 'louvain', method='t-test', corr_method='benjamini-hochberg')
result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
#pd.DataFrame({group + '_' + key: result[key][group]
#    for group in groups for key in ['names', 'scores','logfoldchanges','pvals','pvals_adj']}).head(5)
marker_genes = ['Alpi']
#marker_genes[0]
sc.pl.umap(adata, color=marker_genes[0], size=25)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13439716312056738, 'fraction_numerical': 0.025295508274231677, 'mean_word_length': 2.9985822306238186, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 13, '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': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '367453', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 3405, 'n_tokens_neox': 3203, 'n_words': 672, 'path': 'dataset-notebooks/celseq1_GSE62270_python/celseq1_organoid.ipynb', 'repo_name': 'Jx-b/BUS_notebooks_python', 'size': '8460'}