Hematopoiesis in early mouse embryos [Moignard15]_.

    Returns
    -------
    Annotated data matrix.
961724
def moignard15() -> AnnData:
    
    filename = settings.datasetdir / 'moignard15/nbt.3154-S3.xlsx'
    backup_url = 'http://www.nature.com/nbt/journal/v33/n3/extref/nbt.3154-S3.xlsx'
    adata = sc.read(filename, sheet='dCt_values.txt', backup_url=backup_url)
    # filter out 4 genes as in Haghverdi et al. (2016)
    gene_subset = ~np.in1d(adata.var_names, ['Eif2b1', 'Mrpl19', 'Polr2a', 'Ubc'])
    adata = adata[:, gene_subset]  # retain non-removed genes
    # choose root cell for DPT analysis as in Haghverdi et al. (2016)
    adata.uns['iroot'] = 532  # note that in Matlab/R, counting starts at 1
    # annotate with Moignard et al. (2015) experimental cell groups
    groups_order = ['HF', 'NP', 'PS', '4SG', '4SFG']
    # annotate each observation/cell
    adata.obs['exp_groups'] = [
        next(gname for gname in groups_order if sname.startswith(gname))
        for sname in adata.obs_names]
    # fix the order and colors of names in "groups"
    adata.obs['exp_groups'] = pd.Categorical(adata.obs['exp_groups'],
                                             categories=groups_order)
    adata.uns['exp_groups_colors'] = ['#D7A83E', '#7AAE5D', '#497ABC', '#AF353A', '#765099']
    return adata
1156483
def _cnvkit_prioritize(sample, genes, allele_file, metrics_file):
    
    mdf = pd.read_table(metrics_file)
    mdf.columns = [x.lower() for x in mdf.columns]
    if len(genes) > 0:
        mdf = mdf[mdf["gene"].str.contains("|".join(genes))]
    mdf = mdf[["chromosome", "start", "end", "gene", "log2", "ci_hi", "ci_lo"]]
    adf = pd.read_table(allele_file)
    if len(genes) > 0:
        adf = adf[adf["gene"].str.contains("|".join(genes))]
    if "cn1" in adf.columns and "cn2" in adf.columns:
        adf = adf[["chromosome", "start", "end", "cn", "cn1", "cn2"]]
    else:
        adf = adf[["chromosome", "start", "end", "cn"]]
    df = pd.merge(mdf, adf, on=["chromosome", "start", "end"])
    df = df[df["cn"] != 2]
    if len(df) > 0:
        def passes(row):
            spread = abs(row["ci_hi"] - row["ci_lo"])
            return spread < 0.25
        df["passes"] = df.apply(passes, axis=1)
    df.insert(0, "sample", [sample] * len(df))
    return df
1018496
def preproc_data(data, gene_subset=False, **kwargs):
    
    import uncurl
    from uncurl.preprocessing import log1p, cell_normalize
    from sklearn.decomposition import TruncatedSVD
    data_subset = data
    if gene_subset:
        gene_subset = uncurl.max_variance_genes(data)
        data_subset = data[gene_subset, :]
    tsvd = TruncatedSVD(min(8, data_subset.shape[0] - 1))
    data_tsvd = tsvd.fit_transform(log1p(cell_normalize(data_subset)).T)
    return data_tsvd
1240124
def generate_VJ_junction_transfer_matrices(self):
            
        
        nt2num = {'A': 0, 'C': 1, 'G': 2, 'T': 3}                
        
        #Compute Tvj
        Tvj = {}
        for aa in self.codons_dict.keys():
            current_Tvj = np.zeros((4, 4))
            for init_nt in 'ACGT':
                for codon in self.codons_dict[aa]:
                    current_Tvj[nt2num[codon[2]], nt2num[init_nt]] += self.Rvj[nt2num[codon[2]],nt2num[codon[1]]]*self.Rvj[nt2num[codon[1]],nt2num[codon[0]]] * self.Rvj[nt2num[codon[0]],nt2num[init_nt]]
            Tvj[aa] = current_Tvj
    
        #Compute Svj
        Svj = {}
        for aa in self.codons_dict.keys():
            current_Svj = np.zeros((4, 4))
            for ins_nt in 'ACGT':
                if any([codon.startswith(ins_nt) for codon in self.codons_dict[aa]]):
                    current_Svj[nt2num[ins_nt], :] = self.Rvj[nt2num[ins_nt], :]             
            Svj[aa] = current_Svj
        
        #Compute Dvj               
        Dvj = {}    
        for aa in self.codons_dict.keys():
            current_Dvj = np.zeros((4, 4))
            for init_nt in 'ACGT':
                for codon in self.codons_dict[aa]:
                    current_Dvj[nt2num[codon[2]], nt2num[init_nt]] += self.Rvj[nt2num[codon[1]],nt2num[codon[0]]] * self.Rvj[nt2num[codon[0]],nt2num[init_nt]]
            Dvj[aa] = current_Dvj

        #Compute lTvj
        lTvj = {}
        for aa in self.codons_dict.keys():
            current_lTvj = np.zeros((4, 4))
            for codon in self.codons_dict[aa]:
                current_lTvj[nt2num[codon[2]], nt2num[codon[0]]] += self.Rvj[nt2num[codon[2]],nt2num[codon[1]]]*self.first_nt_bias_insVJ[nt2num[codon[1]]]
            lTvj[aa] = current_lTvj

        #Compute lDvj        
        lDvj = {}    
        for aa in self.codons_dict.keys():
            current_lDvj = np.zeros((4, 4))
            for codon in self.codons_dict[aa]:
                current_lDvj[nt2num[codon[2]], nt2num[codon[0]]] += self.first_nt_bias_insVJ[nt2num[codon[1]]]
            lDvj[aa] = current_lDvj
    

        #Set the attributes
        self.Tvj = Tvj
        self.Svj = Svj
        self.Dvj = Dvj
        self.lTvj = lTvj
        self.lDvj = lDvj
279259
def get_orthology_sources_from_zebrafishmine(self):
        

        # For further documentation you can visit:
        #     http://www.intermine.org/wiki/PythonClient

        # The following two lines will be needed in every python script:
        service = Service("http://zebrafishmine.org/service")

        # Get a new query on the class (table) you will be querying:
        query = service.new_query("Gene")

        # The view specifies the output columns
        query.add_view(
            "primaryIdentifier", "symbol", "homologues.homologue.symbol",
            "homologues.evidence.evidenceCode.abbreviation",
            "homologues.evidence.publications.primaryIdentifier",
            "homologues.evidence.publications.pubMedId",
            "homologues.crossReferences.identifier",
            # only "orthologue" is used
            # "homologues.crossReferences.linkType",
            # only needed if >1 source
            # "homologues.crossReferences.source.name"
        )

        # This query's custom sort order is specified below:
        query.add_sort_order("Gene.name", "ASC")

        # You can edit the constraint values below
        query.add_constraint(
            "homologues.dataSets.name", "=",
            "ZFIN Curated Human, Mouse, Fly, Yeast Orthologue Data Set",
            code="A")
        query.add_constraint(
            "homologues.homologue.organism.name", "=",
            "Homo sapiens", code="B")
        query.add_constraint(
            "homologues.crossReferences.source.name", "=",
            "Gene", code="D")  # NCBIGene
        query.add_constraint("symbol", "=", "*", code="C")

        # Uncomment and edit the code below to specify your own custom logic:
        # query.set_logic("C and A and B and D and D")

        self.files['zmine_ortho_evidence'] = {}
        self.files['zmine_ortho_evidence']['file'] = 'zmine_ortho_evidence.txt'
        file = '/'.join(
            (self.rawdir, self.files['zmine_ortho_evidence']['file']))
        with open(file, 'w', encoding="utf-8", newline='\n') as csvfile:
            filewriter = csv.writer(csvfile, delimiter='\t', quotechar='\"')
            for row in query.rows():
                stuff = [
                    row["primaryIdentifier"],
                    row["symbol"],
                    row["homologues.homologue.symbol"],
                    row["homologues.crossReferences.identifier"],
                    row["homologues.evidence.evidenceCode.abbreviation"],
                    row["homologues.evidence.publications.primaryIdentifier"],
                    row["homologues.evidence.publications.pubMedId"],
                    # row["homologues.crossReferences.linkType"],
                    # row["homologues.crossReferences.source.name"]
                ]
                filewriter.writerow(stuff)

        return
805314
def lms(args):
    
    from random import randint
    from jcvi.graphics.chromosome import HorizontalChromosome

    p = OptionParser(lms.__doc__)
    opts, args, iopts = p.set_image_options(args, figsize="6x6", dpi=300)

    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    # Panel A
    w, h = .7, .35
    ax = fig.add_axes([.15, .6, w, h])

    xdata = [x + randint(-3, 3) for x in range(10, 110, 10)]
    ydata = [x + randint(-3, 3) for x in range(10, 110, 10)]
    ydata[3:7] = ydata[3:7][::-1]
    xydata = zip(xdata, ydata)
    lis = xydata[:3] + [xydata[4]] + xydata[7:]
    lds = xydata[3:7]
    xlis, ylis = zip(*lis)
    xlds, ylds = zip(*lds)
    ax.plot(xlis, ylis, "r-", lw=12, alpha=.3,
                              solid_capstyle="round", solid_joinstyle="round")
    ax.plot(xlds, ylds, "g-", lw=12, alpha=.3,
                              solid_capstyle="round", solid_joinstyle="round")
    ax.plot(xdata, ydata, "k.", mec="k", mfc="w", mew=3, ms=12)
    HorizontalChromosome(root, .15, .15 + w, .57, height=.02, lw=2)
    root.text(.15 + w / 2, .55, "Chromosome location (bp)", ha="center", va="top")

    ax.text(80, 30, "LIS = 7", color="r", ha="center", va="center")
    ax.text(80, 20, "LDS = 4", color="g", ha="center", va="center")
    ax.text(80, 10, "LMS = $max$(LIS, LDS) = 7", ha="center", va="center")
    normalize_lms_axis(ax, xlim=110, ylim=110)

    # Panel B
    w = .37
    p = (0, 45, 75, 110)
    ax = fig.add_axes([.1, .12, w, h])
    xdata = [x for x in range(10, 110, 10)]
    ydata = ydata_orig = [x for x in range(10, 110, 10)]
    ydata = ydata[:4] + ydata[7:] + ydata[4:7][::-1]
    xydata = zip(xdata, ydata)
    lis = xydata[:7]
    xlis, ylis = zip(*lis)
    ax.plot(xlis, ylis, "r-", lw=12, alpha=.3,
                              solid_capstyle="round", solid_joinstyle="round")
    ax.plot(xdata, ydata, "k.", mec="k", mfc="w", mew=3, ms=12)
    ax.vlines(p, 0, 110, colors="beige", lw=3)
    normalize_lms_axis(ax, xlim=110, ylim=110)
    patch = [.1 + w * x / 110. for x in p]
    HorizontalChromosome(root, .1, .1 + w, .09, patch=patch,
                         height=.02, lw=2)
    scaffolds = ("a", "b", "c")
    for i, s in enumerate(scaffolds):
        xx = (patch[i] + patch[i + 1]) / 2
        root.text(xx, .09, s, va="center", ha="center")
    root.text(.1 + w / 2, .04, "LMS($a||b||c$) = 7", ha="center")

    # Panel C
    ax = fig.add_axes([.6, .12, w, h])
    patch = [.6 + w * x / 110. for x in p]
    ydata = ydata_orig
    ax.plot(xdata, ydata, "r-", lw=12, alpha=.3,
                              solid_capstyle="round", solid_joinstyle="round")
    ax.plot(xdata, ydata, "k.", mec="k", mfc="w", mew=3, ms=12)
    ax.vlines(p, [0], [110], colors="beige", lw=3)
    normalize_lms_axis(ax, xlim=110, ylim=110)
    HorizontalChromosome(root, .6, .6 + w, .09, patch=patch,
                         height=.02, lw=2)
    scaffolds = ("a", "-c", "b")
    for i, s in enumerate(scaffolds):
        xx = (patch[i] + patch[i + 1]) / 2
        root.text(xx, .09, s, va="center", ha="center")
    root.text(.6 + w / 2, .04, "LMS($a||-c||b$) = 10", ha="center")

    labels = ((.05, .95, 'A'), (.05, .48, 'B'), (.55, .48, 'C'))
    panel_labels(root, labels)

    normalize_axes(root)

    pf = "lms"
    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
1231775
def poisson_objective(X, m, w):
    
    clusters, cells = w.shape
    genes = X.shape[0]
    #m = m.reshape((X.shape[0], w.shape[0]))
    d = m.dot(w)+eps
    #temp = X/d
    #w_sum = w.sum(1)
    #w2 = w.dot(temp.T)
    #deriv = w_sum - w2.T
    return np.sum(d - X*np.log(d))/genes
670565
def create_dexseq_annotation(gff, count_file):
    
    out_file = count_file + ".ann"
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out:
        with open(tx_out, 'w') as out_handle:
            with open(gff) as in_handle:
                for line in in_handle:
                    cols = line.strip().split("\t")
                    if cols[2] == "exonic_part":
                        exon = [f for f in cols[8].split(";") if f.strip().startswith("exonic_part_number")]
                        gene = [f for f in cols[8].split(";") if f.strip().startswith("gene_id")]
                        exon = exon[0].replace("\"", "").split()[1]
                        gene = gene[0].replace("\"", "").split()[1]
                        length = int(cols[4]) - int(cols[3]) + 1
                        line = "%s:%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (gene, exon, gene,
                                                                    cols[0], cols[3],
                                                                    cols[4],
                                                                    length, cols[6])
                        out_handle.write(line)
1074404
def align_transcriptome(fastq_file, pair_file, ref_file, data):
    
    work_bam = dd.get_work_bam(data)
    base, ext = os.path.splitext(work_bam)
    out_file = base + ".transcriptome" + ext
    if utils.file_exists(out_file):
        data = dd.set_transcriptome_bam(data, out_file)
        return data
    # bwa mem needs phred+33 quality, so convert if it is Illumina
    if dd.get_quality_format(data).lower() == "illumina":
        logger.info("bwa mem does not support the phred+64 quality format, "
                    "converting %s and %s to phred+33.")
        fastq_file = fastq.groom(fastq_file, data, in_qual="fastq-illumina")
        if pair_file:
            pair_file = fastq.groom(pair_file, data, in_qual="fastq-illumina")
    bwa = config_utils.get_program("bwa", data["config"])
    gtf_file = dd.get_gtf_file(data)
    gtf_fasta = index_transcriptome(gtf_file, ref_file, data)
    args = " ".join(_bwa_args_from_config(data["config"]))
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    samtools = config_utils.get_program("samtools", data["config"])
    cmd = ("{bwa} mem {args} -a -t {num_cores} {gtf_fasta} {fastq_file} "
           "{pair_file} ")
    with file_transaction(data, out_file) as tx_out_file:
        message = "Aligning %s and %s to the transcriptome." % (fastq_file, pair_file)
        cmd += "| " + postalign.sam_to_sortbam_cl(data, tx_out_file, name_sort=True)
        do.run(cmd.format(**locals()), message)
    data = dd.set_transcriptome_bam(data, out_file)
    return data
360726
def _dict_seq_locus(list_c, loci_obj, seq_obj):
    
    seqs = defaultdict(set)
    # n = len(list_c.keys())
    for c in list_c.values():
        for l in c.loci2seq:
            [seqs[s].add(c.id) for s in c.loci2seq[l]]

    common = [s for s in seqs if len(seqs[s]) > 1]
    seqs_in_c = defaultdict(float)
    for c in list_c.values():
        for l in c.loci2seq:
            # total = sum([v for v in loci_obj[l].coverage.values()])
            for s in c.loci2seq[l]:
                if s in common:
                    pos = seq_obj[s].pos[l]
                    # cov = 1.0 * loci_obj[l].coverage[pos] / total
                    cov = 1.0 * loci_obj[l].coverage[pos]
                    if seqs_in_c[(s, c.id)] < cov:
                        seqs_in_c[(s, c.id)] = cov
    seqs_in_c = _transform(seqs_in_c)
    return seqs_in_c
833790
def get_refs(genome_build, aligner, galaxy_base, data):
    
    out = {}
    name_remap = {"samtools": "fasta"}
    if genome_build:
        galaxy_config = _get_galaxy_tool_info(galaxy_base)
        for name in [x for x in ("samtools", aligner) if x]:
            galaxy_dt = _get_galaxy_data_table(name, galaxy_config["tool_data_table_config_path"])
            loc_file, need_remap = _get_galaxy_loc_file(name, galaxy_dt, galaxy_config["tool_data_path"],
                                                        galaxy_base)
            cur_ref = _get_ref_from_galaxy_loc(name, genome_build, loc_file, galaxy_dt, need_remap,
                                               galaxy_config, data)
            base = os.path.normpath(utils.add_full_path(cur_ref, galaxy_config["tool_data_path"]))
            # Expand directories unless we are an aligner like minimap2 that uses the seq directory
            if os.path.isdir(base) and not (need_remap and os.path.basename(base) == "seq"):
                indexes = sorted(glob.glob(os.path.join(base, "*")))
            elif name != "samtools":
                indexes = sorted(glob.glob("%s*" % utils.splitext_plus(base)[0]))
            else:
                indexes = []
            name = name_remap.get(name, name)
            out[name] = {}
            if os.path.exists(base) and os.path.isfile(base):
                out[name]["base"] = base
            if indexes:
                out[name]["indexes"] = indexes
            # For references, add compressed inputs and indexes if they exist
            if name == "fasta" and "base" in out[name] and os.path.exists(out[name]["base"] + ".gz"):
                indexes = [out[name]["base"] + ".gz.fai", out[name]["base"] + ".gz.gzi",
                           utils.splitext_plus(out[name]["base"])[0] + ".dict"]
                out[name + "gz"] = {"base": out[name]["base"] + ".gz",
                                    "indexes": [x for x in indexes if os.path.exists(x)]}
        # add additional indices relative to the base
        if tz.get_in(["fasta", "base"], out):
            ref_dir, ref_filebase = os.path.split(out["fasta"]["base"])
            rtg_dir = os.path.normpath(os.path.join(ref_dir, os.path.pardir, "rtg",
                                                    "%s.sdf" % (os.path.splitext(ref_filebase)[0])))
            out["rtg"] = {"base": os.path.join(rtg_dir, "mainIndex"),
                          "indexes": [x for x in glob.glob(os.path.join(rtg_dir, "*"))
                                      if not x.endswith("/mainIndex")]}
            twobit = os.path.normpath(os.path.join(ref_dir, os.path.pardir, "ucsc",
                                                   "%s.2bit" % (os.path.splitext(ref_filebase)[0])))
            if os.path.exists(twobit):
                out["twobit"] = twobit
    return out
17183
def sample_annotation(data):
    
    names = data["rgnames"]['sample']
    tools = dd.get_expression_caller(data)
    work_dir = os.path.join(dd.get_work_dir(data), "mirbase")
    out_dir = os.path.join(work_dir, names)
    utils.safe_makedir(out_dir)
    out_file = op.join(out_dir, names)
    if dd.get_mirbase_hairpin(data):
        mirbase = op.abspath(op.dirname(dd.get_mirbase_hairpin(data)))
        if utils.file_exists(data["collapse"]):
            data['transcriptome_bam'] = _align(data["collapse"],
                                               dd.get_mirbase_hairpin(data),
                                               out_file,
                                               data)
            data['seqbuster'] = _miraligner(data["collapse"], out_file,
                                            dd.get_species(data),
                                            mirbase,
                                            data['config'])
            data["mirtop"] = _mirtop(data['seqbuster'],
                                     dd.get_species(data),
                                     mirbase,
                                     out_dir,
                                     data['config'])
        else:
            logger.debug("Trimmed collapsed file is empty for %s." % names)
    else:
        logger.debug("No annotation file from miRBase.")

    sps = dd.get_species(data) if dd.get_species(data) else "None"
    logger.debug("Looking for mirdeep2 database for %s" % names)
    if file_exists(op.join(dd.get_work_dir(data), "mirdeep2", "novel", "hairpin.fa")):
        data['seqbuster_novel'] = _miraligner(data["collapse"], "%s_novel" % out_file, sps,
                                              op.join(dd.get_work_dir(data),
                                                      "mirdeep2", "novel"),
                                              data['config'])

    if "trna" in tools:
        data['trna'] = _mint_trna_annotation(data)

    data = spikein.counts_spikein(data)
    return [[data]]
991145
def bamtag(sam):
    ''' 
    '''
    from pysam import AlignmentFile

    start_time = time.time()

    sam_file = open_bamfile(sam)
    out_file = AlignmentFile("-", "wh", template=sam_file)
    track = sam_file.fetch(until_eof=True)

    # peek at first alignment to determine the annotations
    if is_python3():
        queryalignment = next(track)
    else:
        queryalignment = track.next()
    annotations = detect_alignment_annotations(queryalignment)
    track = itertools.chain([queryalignment], track)

    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)

    for count, aln in enumerate(track, start=1):
        if count and not count % 1000000:
            logger.info("Processed %d alignments." % count)

        match = parser_re.match(aln.qname)
        tags = aln.tags

        if "cellular" in annotations:
            aln.tags += [('XC', match.group('CB'))]
        if "molecular" in annotations:
            aln.tags += [('RX', match.group('MB'))]
        if "sample" in annotations:
            aln.tags += [('XS', match.group('SB'))]

        out_file.write(aln)

    total_time = time.time() - start_time
    logger.info('BAM tag conversion done - {:.3}s, {:,} alns/min'.format(total_time, int(60. * count / total_time)))
    logger.info("Processed %d alignments." % count)
930151
def idsKEGG(organism):
    
    ORG=urlopen("http://rest.kegg.jp/list/"+organism).read()
    ORG=ORG.split("\n")
    final=[]
    for k in ORG:
        final.append(k.split("\t"))
    df=pd.DataFrame(final[0:len(final)-1])[[0,1]]
    df.columns=['KEGGid','description']
    field = pd.DataFrame(df['description'].str.split(';',1).tolist())[0]
    field = pd.DataFrame(field)
    df = pd.concat([df[['KEGGid']],field],axis=1)
    df.columns=['KEGGid','gene_name']
    df=df[['gene_name','KEGGid']]
    return df
1008613
def _gids_to_genes(gids, ssm_locs, cnv_ssms, data):
    
    locs = collections.defaultdict(set)
    for gid in gids:
        cur_locs = []
        try:
            cur_locs.append(ssm_locs[gid])
        except KeyError:
            for ssm_loc in cnv_ssms.get(gid, []):
                cur_locs.append(ssm_locs[ssm_loc])
        for chrom, pos in cur_locs:
            locs[chrom].add(pos)
    genes = set([])
    with tx_tmpdir(data) as tmpdir:
        chrom_prefix = "chr" if next(ref.file_contigs(dd.get_ref_file(data))).name.startswith("chr") else ""
        loc_file = os.path.join(tmpdir, "battenberg_find_genes.bed")
        with open(loc_file, "w") as out_handle:
            for chrom in sorted(locs.keys()):
                for loc in sorted(list(locs[chrom])):
                    out_handle.write("%s%s\t%s\t%s\n" % (chrom_prefix, chrom, loc - 1, loc))
        ann_file = annotate.add_genes(loc_file, data, max_distance=10000)
        for r in pybedtools.BedTool(ann_file):
            for gene in r.name.split(","):
                if gene != ".":
                    genes.add(gene)
    return sorted(list(genes))
561650
def matrix(records, scafstats_path):
    '''
    
    '''
    ncbi = ete3.NCBITaxa()
    
    # # Single file version 
    # contigs_lineages = {r.id: r.description.strip().partition(' ')[2].split('|')[2] for r in records}
    # df = pd.read_csv(scafstats_path, sep='\t').set_index('#name')
    # contigs_counts = pd.Series(df.assignedReads.values, index=df.index).to_dict()
    # lineages_counts = defaultdict(int)
    
    # for contig, lineage in contigs_lineages.items():
    #     lineages_counts[lineage] += contigs_counts.get(contig, 0)

    # lineages_counts_sorted = dict(reversed(sorted(lineages_counts.items(), key=lambda x: x[1])))
    # print(pd.DataFrame(lineages_counts_sorted, index=[scafstats_path]).transpose().to_csv())

    scafstats_paths = {fn.replace('.scafstats', ''): f'{scafstats_path}/{fn}'
                       for fn in os.listdir(scafstats_path)
                       if fn.endswith('.scafstats')}
    contigs_lineages = {r.id: r.description.strip().partition(' ')[2].split('|')[3] for r in records}
    samples = []
    
    for scafstat, path in scafstats_paths.items():
        raw_df = pd.read_csv(path, sep='\t').set_index('#name')
        contigs_counts = pd.Series(raw_df.assignedReads.values, index=raw_df.index).to_dict()
        lineages_counts = defaultdict(int)
        for contig, lineage in contigs_lineages.items():
            lineages_counts[lineage] += contigs_counts.get(contig, 0)
        counts_df = pd.DataFrame(lineages_counts, index=[scafstat]).transpose()
        samples.append(counts_df)

    return pd.concat(samples, axis=1, sort=False).rename_axis('lineage')
356599
def demultiplex_cells(fastq, out_dir, readnumber, prefix, cb_histogram, 
                      cb_cutoff):
    ''' 
    '''
    annotations = detect_fastq_annotations(fastq)
    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)
    readstring = "" if not readnumber else "_R{}".format(readnumber)
    filestring = "{prefix}{sample}{readstring}.fq"
    cb_set = set()
    if cb_histogram:
        cb_set = get_cb_depth_set(cb_histogram, cb_cutoff)
    sample_set = set()
    batch = collections.defaultdict(list)
    parsed = 0
    safe_makedir(out_dir)
    for read in read_fastq(fastq):
        parsed += 1
        match = parser_re.search(read).groupdict()
        sample = match['CB']
        if cb_set and sample not in cb_set:
            continue
        sample_set.add(sample)
        batch[sample].append(read)
        # write in batches to avoid opening up file handles repeatedly
        if not parsed % 10000000:
            for sample, reads in batch.items():
                out_file = os.path.join(out_dir, filestring.format(**locals()))
                with open(out_file, "a") as out_handle:
                    for read in reads:
                        out_handle.write(read)
            batch = collections.defaultdict(list)

    for sample, reads in batch.items():
        out_file = os.path.join(out_dir, filestring.format(**locals()))
        with open(out_file, "a") as out_handle:
            for read in reads:
                out_handle.write(read)