Read a 10X genomics compressed tarball containing expression data.
        
        Note: common prefix patterns:
        - "filtered_gene_bc_matrices/[annotations]/"
        - "filtered_matrices_mex/[annotations]/"

        TODO: docstring
161104
def read_10xgenomics(cls, tarball_fpath: str, prefix: str,
                         use_ensembl_ids: bool = False):
        """"""

        _LOGGER.info('Reading file: %s', tarball_fpath)

        with tarfile.open(tarball_fpath, mode='r:gz') as tf:
            ti = tf.getmember('%smatrix.mtx' % prefix)
            with tf.extractfile(ti) as fh:
                mtx = scipy.io.mmread(fh)

            ti = tf.getmember('%sgenes.tsv' % prefix)
            with tf.extractfile(ti) as fh:
                wrapper = io.TextIOWrapper(fh, encoding='ascii')
                i = 1
                if use_ensembl_ids:
                    i = 0
                gene_names = \
                        [row[i] for row in csv.reader(wrapper, delimiter='\t')]

            ti = tf.getmember('%sbarcodes.tsv' % prefix)
            with tf.extractfile(ti) as fh:
                wrapper = io.TextIOWrapper(fh, encoding='ascii')
                barcodes = \
                        [row[0] for row in csv.reader(wrapper, delimiter='\t')]


            assert mtx.shape[0] == len(gene_names)
            assert mtx.shape[1] == len(barcodes)
        
        _LOGGER.info('Matrix dimensions: %s', str(mtx.shape))
        X = mtx.todense()
        matrix = cls(X=X, genes=gene_names, cells=barcodes)
        
        return matrix
759328
def sb_filter(fastq, bc, cores, nedit):
    ''' 
    '''
    barcodes = set(sb.strip() for sb in bc)
    if nedit == 0:
        filter_sb = partial(exact_sample_filter2, barcodes=barcodes)
    else:
        barcodehash = MutationHash(barcodes, nedit)
        filter_sb = partial(correcting_sample_filter2, barcodehash=barcodehash)
    p = multiprocessing.Pool(cores)

    chunks = tz.partition_all(10000, read_fastq(fastq))
    bigchunks = tz.partition_all(cores, chunks)
    for bigchunk in bigchunks:
        for chunk in p.map(filter_sb, list(bigchunk)):
            for read in chunk:
                sys.stdout.write(read)
984494
def fasttagcount(sam, out, genemap, positional, minevidence, cb_histogram, 
                 cb_cutoff, subsample, parse_tags, gene_tags, umi_matrix):
    ''' 
    '''
    from pysam import AlignmentFile

    from io import StringIO
    import pandas as pd

    from utils import weigh_evidence

    if sam.endswith(".sam"):
        logger.error("To use the fasttagcount subcommand, the alignment file must be a "
                     "coordinate sorted, indexed BAM file.")
        sys.exit(1)

    logger.info('Reading optional files')

    gene_map = None
    if genemap:
        with open(genemap) as fh:
            try:
                gene_map = dict(p.strip().split() for p in fh)
            except ValueError:
                logger.error('Incorrectly formatted gene_map, need to be tsv.')
                sys.exit()

    if positional:
        tuple_template = '{0},{1},{2},{3}'
    else:
        tuple_template = '{0},{1},{3}'

    if not cb_cutoff:
        cb_cutoff = 0

    if cb_histogram and cb_cutoff == "auto":
        cb_cutoff = guess_depth_cutoff(cb_histogram)

    cb_cutoff = int(cb_cutoff)

    cb_hist = None
    filter_cb = False
    if cb_histogram:
        cb_hist = pd.read_csv(cb_histogram, index_col=0, header=-1, squeeze=True, sep="\t")
        total_num_cbs = cb_hist.shape[0]
        cb_hist = cb_hist[cb_hist > cb_cutoff]
        logger.info('Keeping {} out of {} cellular barcodes.'.format(cb_hist.shape[0], total_num_cbs))
        filter_cb = True

    parser_re = re.compile('.*:CELL_(?P<CB>.*):UMI_(?P<MB>.*)')

    if subsample:
        logger.info('Creating reservoir of subsampled reads ({} per cell)'.format(subsample))
        start_sampling  = time.time()

        reservoir = collections.defaultdict(list)
        cb_hist_sampled = 0 * cb_hist
        cb_obs = 0 * cb_hist

        track = stream_bamfile(sam)
        current_read = 'none_observed_yet'
        for i, aln in enumerate(track):
            if aln.qname == current_read:
                continue

            current_read = aln.qname

            if parse_tags:
                CB = aln.get_tag('CR')
            else:
                match = parser_re.match(aln.qname)
                CB = match.group('CB')

            if CB not in cb_hist.index:
                continue

            cb_obs[CB] += 1
            if len(reservoir[CB]) < subsample:
                reservoir[CB].append(i)
                cb_hist_sampled[CB] += 1
            else:
                s = pd.np.random.randint(0, cb_obs[CB])
                if s < subsample:
                    reservoir[CB][s] = i

        index_filter = set(itertools.chain.from_iterable(reservoir.values()))
        sam_file.close()
        sampling_time = time.time() - start_sampling
        logger.info('Sampling done - {:.3}s'.format(sampling_time))

    evidence = collections.defaultdict(lambda: collections.defaultdict(float))
    bare_evidence = collections.defaultdict(float)
    logger.info('Tallying evidence')
    start_tally = time.time()

    sam_mode = 'r' if sam.endswith(".sam") else 'rb'
    sam_file = AlignmentFile(sam, mode=sam_mode)
    transcript_map = collections.defaultdict(set)
    sam_transcripts = [x["SN"] for x in sam_file.header["SQ"]]
    if gene_map:
        for transcript, gene in gene_map.items():
            if transcript in sam_transcripts:
                transcript_map[gene].add(transcript)
    else:
        for transcript in sam_transcripts:
            transcript_map[transcript].add(transcript)
    missing_transcripts = set()
    alignments_processed = 0
    unmapped = 0
    kept = 0
    nomatchcb = 0
    current_read = 'none_observed_yet'
    current_transcript = None
    count_this_read = True
    transcripts_processed = 0
    genes_processed = 0
    cells = list(cb_hist.index)
    targets_seen = set()

    if umi_matrix:
        bare_evidence_handle = open(umi_matrix, "w")
        bare_evidence_handle.write(",".join(["gene"] + cells) + "\n")

    with open(out, "w") as out_handle:
        out_handle.write(",".join(["gene"] + cells) + "\n")
        for gene, transcripts in transcript_map.items():
            for transcript in transcripts:
                for aln in sam_file.fetch(transcript):
                    alignments_processed += 1

                    if aln.is_unmapped:
                        unmapped += 1
                        continue

                    if gene_tags and not aln.has_tag('GX'):
                        unmapped += 1
                        continue

                    if aln.qname != current_read:
                        current_read = aln.qname
                        if subsample and i not in index_filter:
                            count_this_read = False
                            continue
                        else:
                            count_this_read = True
                    else:
                        if not count_this_read:
                            continue

                    if parse_tags:
                        CB = aln.get_tag('CR')
                    else:
                        match = parser_re.match(aln.qname)
                        CB = match.group('CB')

                    if filter_cb:
                        if CB not in cb_hist.index:
                            nomatchcb += 1
                            continue

                    if parse_tags:
                        MB = aln.get_tag('UM')
                    else:
                        MB = match.group('MB')

                    if gene_tags:
                        target_name = aln.get_tag('GX').split(',')[0]
                    else:
                        txid = sam_file.getrname(aln.reference_id)
                        if gene_map:
                            if txid in gene_map:
                                target_name = gene_map[txid]
                            else:
                                missing_transcripts.add(txid)
                                continue
                        else:
                            target_name = txid
                    targets_seen.add(target_name)

                    # Scale evidence by number of hits
                    evidence[CB][MB] += weigh_evidence(aln.tags)
                    bare_evidence[CB] += weigh_evidence(aln.tags)
                    kept += 1
                transcripts_processed += 1
                if not transcripts_processed % 1000:
                    logger.info("%d genes processed." % genes_processed)
                    logger.info("%d transcripts processed." % transcripts_processed)
                    logger.info("%d alignments processed." % alignments_processed)

            earray = []
            for cell in cells:
                umis = [1 for _, v in evidence[cell].items() if v >= minevidence]
                earray.append(str(sum(umis)))
            out_handle.write(",".join([gene] + earray) + "\n")
            earray = []
            if umi_matrix:
                for cell in cells:
                    earray.append(str(int(bare_evidence[cell])))
                bare_evidence_handle.write(",".join([gene] + earray) + "\n")

            evidence = collections.defaultdict(lambda: collections.defaultdict(int))
            bare_evidence = collections.defaultdict(int)
            genes_processed += 1

    if umi_matrix:
        bare_evidence_handle.close()

    # fill dataframe with missing values, sort and output
    df = pd.read_csv(out, index_col=0, header=0)
    targets = pd.Series(index=set(transcript_map.keys()))
    targets = targets.sort_index()
    df = df.reindex(targets.index.values, fill_value=0)
    df = df.sort_index()
    df.to_csv(out)

    if umi_matrix:
        df = pd.read_csv(umi_matrix, index_col=0, header=0)
        df = df.reindex(targets.index.values, fill_value=0)
        df = df.sort_index()
        df.to_csv(umi_matrix)
1200675
def explore(args):
    
    logger.info("reading sequeces")
    data = load_data(args.json)
    logger.info("get sequences from json")
    #get_sequences_from_cluster()
    c1, c2 = args.names.split(",")
    seqs, names = get_sequences_from_cluster(c1, c2, data[0])
    loci = get_precursors_from_cluster(c1, c2, data[0])
    logger.info("map all sequences to all loci")
    print("%s" % (loci))
    map_to_precursors(seqs, names, loci, os.path.join(args.out, "map.tsv"), args)
    #map_sequences_w_bowtie(sequences, precursors)
    logger.info("plot sequences on loci")
    #get_matrix_position()
    #plot_sequences()
    logger.info("Done")
947891
def read_alignment(out_sam, loci, seqs, out_file):
    """"""
    hits = defaultdict(list)
    with open(out_file, "w") as out_handle:
        samfile = pysam.Samfile(out_sam, "r")
        for a in samfile.fetch():
            if not a.is_unmapped:
                nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                a = makeBED(a)
                if not a:
                    continue
                ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
        for hit in hits.values():
            nm = hit[0][0]
            for l in hit:
                if nm == l[0]:
                    print(l[1], file=out_handle)
    return out_file
178254
def seqs_from_fastacmd(acc_list, blast_db,is_protein=True):
    """"""
    fasta_cmd_res = fasta_cmd_get_seqs(acc_list, blast_db=blast_db, \
        is_protein=is_protein)
    recs = FastaCmdFinder(fasta_cmd_res['StdOut'])
    result = {}
    for rec in recs:
        try:
            result[rec[0][1:].strip()] = ''.join(map(strip, rec[1:]))
        except IndexError:  #maybe we didn't get a sequence?
            pass
    fasta_cmd_res.cleanUp()
    return result
314623
def prep_fastq_inputs(in_files, data):
    
    if len(in_files) == 1 and _is_bam_input(in_files):
        out = _bgzip_from_bam(in_files[0], data["dirs"], data)
    elif len(in_files) == 1 and _is_cram_input(in_files):
        out = _bgzip_from_cram(in_files[0], data["dirs"], data)
    elif len(in_files) in [1, 2] and _ready_gzip_fastq(in_files, data):
        out = _symlink_in_files(in_files, data)
    else:
        if len(in_files) > 2:
            fpairs = fastq.combine_pairs(in_files)
            pair_types = set([len(xs) for xs in fpairs])
            assert len(pair_types) == 1
            fpairs.sort(key=lambda x: os.path.basename(x[0]))
            organized = [[xs[0] for xs in fpairs]]
            if len(fpairs[0]) > 1:
                organized.append([xs[1] for xs in fpairs])
            in_files = organized
        parallel = {"type": "local", "num_jobs": len(in_files),
                    "cores_per_job": max(1, data["config"]["algorithm"]["num_cores"] // len(in_files))}
        inputs = [{"in_file": x, "read_num": i, "dirs": data["dirs"], "config": data["config"],
                   "is_cwl": "cwl_keys" in data,
                   "rgnames": data["rgnames"]}
                  for i, x in enumerate(in_files) if x]
        out = run_multicore(_bgzip_from_fastq_parallel, [[d] for d in inputs], data["config"], parallel)
    return out
78754
def build_dataset(filename, sequence_length, alphabet_size, max_lines=-1):
    """"""
    with open(filename) as f:
        content = f.readlines()
    content = content[:max_lines]
    content = [line for line in content if len(line) > 2]   # Remove blank lines
    seqs = np.zeros((sequence_length, len(content), alphabet_size))
    for ix, line in enumerate(content):
        padded_line = (line + " " * sequence_length)[:sequence_length]
        seqs[:, ix, :] = string_to_one_hot(padded_line, alphabet_size)
    return seqs
423223
def _wrpy_ncbi_gene_nts(fout_py, geneid2nt, log):
        """"""
        num_genes = len(geneid2nt)
        with open(fout_py, 'w') as ofstrm:
            docstr = "Data downloaded from NCBI Gene converted into Python namedtuples."
            ofstrm.write('"""{PYDOC}"""\n\n'.format(PYDOC=docstr))
            ofstrm.write("from collections import namedtuple\n\n")
            ofstrm.write('WRITTEN = "{DATE}"'.format(
                DATE=re.sub('-', '_', str(datetime.date.today()))))
            ofstrm.write(' # {N} items\n\n'.format(N=num_genes))
            ntd = next(iter(geneid2nt.values())) # Access one dictionary value in Python 2
            ofstrm.write("#pylint: disable=line-too-long,too-many-lines,invalid-name\n")
            ofstrm.write("{NtName} = namedtuple('{NtName}', '{FLDS}')\n\n".format(
                NtName=type(ntd).__name__, FLDS=' '.join(ntd._fields)))
            ofstrm.write("GENEID2NT = {{ # {N:,} items\n".format(N=num_genes))
            for geneid, ntd in sorted(geneid2nt.items(), key=lambda t: t[0]):
                ofstrm.write("    {GeneID} : {NT},\n".format(GeneID=geneid, NT=ntd))
            ofstrm.write("}\n")
            log.write("  {N:9} geneids WROTE: {PY}\n".format(N=num_genes, PY=fout_py))
1336550
def parse(self, limit=None):
        

        files_to_download, ftp = self._get_file_list(
            self.files['anat_entity']['path'],
            self.files['anat_entity']['pattern'])
        for dlname in files_to_download:
            localfile = '/'.join((self.rawdir, dlname))
            with gzip.open(localfile, 'rt', encoding='ISO-8859-1') as fh:
                LOG.info("Processing %s", localfile)
                self._parse_gene_anatomy(fh, limit)
        return
1081621
def parse_rsem_report (self, f):
        
        data = dict()
        multimapping_hist = dict()
        in_hist = False
        for l in f['f'].splitlines():
            s = l.split()
            if len(s) > 3:
                # Line: N0 N1 N2 N_tot
                # N0, number of unalignable reads;
                # N1, number of alignable reads; (nUnique + nMulti)
                # N2, number of filtered reads due to too many alignments
                # N_tot = N0 + N1 + N2
                data['Unalignable'] = int(s[0])
                data['Alignable'] = int(s[1])
                data['Filtered'] = int(s[2])
                data['Total'] = int(s[3])
                data['alignable_percent'] = (float(s[1]) / float(s[3])) * 100.0
            elif len(s) == 3:
                # Line: nUnique nMulti nUncertain
                # nUnique, number of reads aligned uniquely to a gene
                # nMulti, number of reads aligned to multiple genes; nUnique + nMulti = N1;
                # nUncertain, number of reads aligned to multiple locations in the given reference sequences, which include isoform-level multi-mapping reads
                data['Unique'] = int(s[0])
                data['Multi'] = int(s[1])
                data['Uncertain'] = int(s[2])
            elif len(s) == 2:
                # Number of gene alignments and read count
                if in_hist or int(s[0]) == 0:
                    in_hist = True
                    try:
                        multimapping_hist[int(s[0])] = int(s[1])
                    except ValueError:
                        pass
            else:
                break
        try:
            assert data['Unique'] + data['Multi'] == data['Alignable']
        except AssertionError:
            log.warn("Unique + Multimapping read counts != alignable reads! '{}'".format(f['fn']))
            return None
        except KeyError:
            log.warn("Error parsing RSEM counts file '{}'".format(f['fn']))
            return None

        # Save parsed data
        if len(data) > 0:
            if f['s_name'] in self.rsem_mapped_data:
                log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
            self.rsem_mapped_data[f['s_name']] = data
            self.add_data_source(f)
        if len(multimapping_hist) > 0:
            self.rsem_multimapping_data[f['s_name']] = multimapping_hist
616859
def read_gtf_line(cols, field="name"):
    """"""
    field = field.lower()
    try:
        group = cols[2]
        attrs = cols[8].split(";")
        name = [attr.strip().split(" ")[1] for attr in attrs if attr.strip().split(" ")[0].lower().endswith(field)]
        if not name:
            name = [attr.strip().split(" ")[1] for attr in attrs if attr.strip().split(" ")[0].lower().endswith("gene_id")]
        if not name:
            name = ["None"]
        biotype = [attr.strip().split(" ")[1] for attr in attrs if attr.strip().split(" ")[0].lower().endswith("biotype")]
        if biotype:
            group = biotype[0]
        c = cols[0]
        s = int(cols[3])
        e = int(cols[4])
        st = cols[6]
        return [c, s, e, st, group, name[0]]
    except(Exception, e):
        logger.error(cols)
        logger.error("File is not in correct format")
        logger.error("Expect chr source feature start end . strand attributes")
        logger.error("Attributes are 'gene_name SNCA; gene_id ENSG; '")
        logger.error("The 3rd column is used as type of small RNA (like miRNA)")
        logger.error("at least should contains '; *name NAME; '")
        logger.error(e)
        raise
291260
def _guess_seqtype_from_file(handle):
    ""
    if isinstance(handle, basestring):
        handle = StringIO(handle)

    for line in handle:
        if not line.strip():
            continue
        if line.lstrip().split()[0] in ('LOCUS', 'FEATURES', 'source', 'CDS',
                                        'gene'):
            return 'genbank'
        if len(line) > 2 and line[:3] in ('ID ', 'FT '):
            return 'embl'
        if line.startswith('>'):
            return 'fasta'
    handle.seek(0)
    import string
    from Bio.Data import IUPACData as iupac
    all_input_letters = set(handle.read().lower())
    all_valid = set(string.digits)
    all_valid.update(set(iupac.protein_letters.lower()))
    all_valid.update(set(iupac.unambiguous_dna_letters.lower()))
    all_valid.update(set('- \n'))
    if all_valid.issuperset(all_input_letters):
        return 'fasta'

    raise ValueError("Failed to guess format for input")
1088047
def get_assoc_ncbi_taxids(taxids, force_dnld=False, loading_bar=True, **kws):
    """"""
    fin = kws['gene2go'] if 'gene2go' in kws else os.path.join(os.getcwd(), "gene2go")
    dnld_ncbi_gene_file(fin, force_dnld, loading_bar=loading_bar)
    return read_ncbi_gene2go(fin, taxids, **kws)
1053360
def annotated_reads(txome,sorted_read_stream,output,args):
   if not os.path.exists(args.tempdir+'/'+output):
      os.makedirs(args.tempdir+'/'+output)
   ts = OrderedStream(iter(txome.transcript_stream()))
   rs = OrderedStream(GPDStream(sorted_read_stream))
   mls = MultiLocusStream([ts,rs])
   
   aof = gzip.open(args.tempdir+'/'+output+'/annotated.txt.gz','w')
   of1 = gzip.open(args.tempdir+'/'+output+'/unannotated_multiexon.gpd.gz','w')
   of2 = gzip.open(args.tempdir+'/'+output+'/partial_annotated_multiexon.gpd.gz','w')
   of3 = gzip.open(args.tempdir+'/'+output+'/unannotated_singleexon.gpd.gz','w')
   z = 0
   for ml in mls:
      refs, reads = ml.payload
      if len(refs) == 0: continue
      """Check and see if we have single exons annotated"""
      single_exon_refs = [x for x in refs if x.get_exon_count()==1]
      single_exon_reads = [x for x in reads if x.get_exon_count()==1]
      #print str(len(single_exon_refs))+" "+str(len(single_exon_reads))
      unannotated_single_exon_reads = []
      for seread in single_exon_reads:
         ovs  = [(x.exons[0].length,seread.exons[0].length,x.exons[0].overlap_size(seread.exons[0]),x) for x in single_exon_refs if x.exons[0].overlaps(seread.exons[0])]
         #check for minimum overlap
         ovs = [x for x in ovs if x[2] >= args.single_exon_minimum_overlap]
         ovs = [x for x in ovs if float(x[2])/float(max(x[0],x[1])) >= args.single_exon_mutual_overlap]
         ovs = sorted(ovs,key=lambda x: float(x[2])/float(max(x[0],x[1])),reverse=True)
         if len(ovs) == 0: 
            unannotated_single_exon_reads.append(seread)
            continue
         best_ref = ovs[0][3]
         aof.write(seread.name+"\t"+best_ref.name+"\tSE"+"\n")

      """May want an optional check for any better matches among exons"""
      #now we can look for best matches among multi-exon transcripts
      reads = [x for x in reads if x.get_exon_count() > 1]
      multiexon_refs = [x for x in refs if x.get_exon_count() > 1]
      unannotated_multi_exon_reads = []
      partial_annotated_multi_exon_reads = []
      for read in reads:
         # we dont' need to have multiple exons matched. one is enough to call
         ovs = [y for y in [(x,x.exon_overlap(read,
                                              multi_minover=args.multi_exon_minimum_overlap,
                                              multi_endfrac=args.multi_exon_end_frac,
                                              multi_midfrac=args.multi_exon_mid_frac,
                                              multi_consec=False)) for x in multiexon_refs] if y[1]]
         for o in ovs: o[1].analyze_overs()
         full =  sorted([x for x in ovs if x[1].is_subset()==1],
                        key = lambda y: float(y[1].overlap_size())/float(max(y[1].tx_obj1.length,y[1].tx_obj2.length)),
                        reverse=True
                       )
         if len(full) > 0:
            aof.write(read.name+"\t"+full[0][0].name+"\tfull"+"\n")
            continue
         subset =  sorted([x for x in ovs if x[1].is_subset()==2],
                          key = lambda y: (y[1].match_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         if len(subset) > 0:
            aof.write(read.name+"\t"+subset[0][0].name+"\tpartial"+"\n")
            continue
         #check for supersets
         superset = sorted([x for x in ovs if x[1].is_subset()==3],
                          key = lambda y: (y[1].match_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         if len(superset) > 0:
            partial_annotated_multi_exon_reads.append((read,superset[0][0]))
            #print read.name+"\t"+superset[0][0].name+"\tsuper"
            continue
         #check for noncompatible overlaps
         overset = sorted([x for x in ovs if x[1].match_exon_count > 0],
                          key = lambda y: (y[1].consecutive_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         #print [(x[1].consecutive_exon_count(), x[1].min_overlap_fraction()) for x in overset]
         if len(overset) > 0:
            partial_annotated_multi_exon_reads.append((read,overset[0][0]))
            #print read.name+"\t"+overset[0][0].name+"\tover"
            continue
         unannotated_multi_exon_reads.append(read)
      """Now we have partially annotated multi and unannotated multi and unannotated single"""
      if len(unannotated_multi_exon_reads) > 0:
         sys.stderr.write(str(z)+" "+str(len(unannotated_multi_exon_reads))+"   \r")
         d = Deconvolution(downsample(unannotated_multi_exon_reads,args.downsample_locus))
         groups = d.parse(tolerance=20,downsample=args.downsample)
         for tx in groups:
            z+=1
            of1.write(tx.get_gpd_line()+"\n")
      if len(partial_annotated_multi_exon_reads) > 0:
         sys.stderr.write(str(z)+" "+str(len(partial_annotated_multi_exon_reads))+"   \r")
         ### set the direction of the transcript
         for v in partial_annotated_multi_exon_reads:
            v[0].set_strand(v[1].direction)
         ### set the gene name of the transcript
         for v in partial_annotated_multi_exon_reads:
            v[0].set_gene_name(v[1].gene_name)         
         d = Deconvolution(downsample([x[0] for x in partial_annotated_multi_exon_reads],args.downsample_locus))
         groups = d.parse(tolerance=20,downsample=args.downsample,use_gene_names=True)
         for tx in groups:
            z += 1
            of2.write(tx.get_gpd_line()+"\n")
      """Do the unannotated single exon reads"""
      g = Graph()
      for r in unannotated_single_exon_reads:
         if len([x for x in partial_annotated_multi_exon_reads if x[0].overlaps(r)]) > 0: continue
         if len([x for x in unannotated_multi_exon_reads if x.overlaps(r)]) > 0: continue
         if len([x for x in txome.transcripts if x.overlaps(r)]) > 0: continue
         g.add_node(Node(r))
      for i in range(0,len(g.nodes)):
         for j in range(0,len(g.nodes)):
            if i == j: continue
            if g.nodes[i].payload.overlaps(g.nodes[j].payload):
               g.add_edge(Edge(g.nodes[i],g.nodes[j]))
      g.merge_cycles()
      for r in g.roots:
         se = []
         se += r.payload_list
         children = g.get_children(r)
         for child in children:  se += child.payload_list
         rng = GenomicRange(se[0].exons[0].chr,
                            min([x.exons[0].start for x in se]),
                            max([x.exons[0].end for x in se]))
         tx = Transcript([rng],Transcript.Options(direction='+'))
         of3.write(tx.get_gpd_line()+"\n")

   of1.close()
   of2.close()
   of3.close()
   sys.stderr.write("\n")      
   """"""
   sys.stderr.write("sort transcriptome\n")
   sort_gpd(args.tempdir+'/'+output+'/partial_annotated_multiexon.gpd.gz',args.tempdir+'/'+output+'/partial_annotated_multiexon.sorted.gpd.gz',args)
   sort_gpd(args.tempdir+'/'+output+'/unannotated_multiexon.gpd.gz',args.tempdir+'/'+output+'/unannotated_multiexon.sorted.gpd.gz',args)
   sort_gpd(args.tempdir+'/'+output+'/unannotated_singleexon.gpd.gz',args.tempdir+'/'+output+'/unannotated_singleexon.sorted.gpd.gz',args)
   """We still have the unannotated single exon reads to deal with"""
216013
def _get_gene_history(self, limit):
        
        src_key = 'gene_history'
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        model = Model(graph)
        LOG.info("Processing Gene records")
        line_counter = 0
        myfile = '/'.join((self.rawdir, self.files[src_key]['file']))
        LOG.info("FILE: %s", myfile)
        col = self.files[src_key]['columns']
        with gzip.open(myfile, 'rb') as tsv:
            row = tsv.readline().decode().strip().split('\t')
            row[0] = row[0][1:]  # strip comment
            if col != row:
                LOG.info(
                    '%s\nExpected Headers:\t%s\nRecived Headers:\t %s\n',
                    src_key, col, row)

            for line in tsv:
                # skip comments
                row = line.decode().strip().split('\t')
                if row[0][0] == '#':
                    continue
                # (tax_num, gene_num, discontinued_num, discontinued_symbol,
                # discontinued_date) = line.split('\t')

                # set filter=None in init if you don't want to have a filter
                # if self.id_filter is not None:
                #     if ((self.id_filter == 'taxids' and \
                #          (int(tax_num) not in self.tax_ids))
                #             or (self.id_filter == 'geneids' and \
                #                 (int(gene_num) not in self.gene_ids))):
                #         continue
                #  end filter

                gene_num = row[col.index('GeneID')].strip()
                discontinued_num = row[col.index('Discontinued_GeneID')].strip()

                if gene_num == '-' or discontinued_num == '-':
                    continue

                if self.test_mode and int(gene_num) not in self.gene_ids:
                    continue

                tax_num = row[col.index('tax_id')].strip()
                if not self.test_mode and tax_num not in self.tax_ids:
                    continue

                line_counter += 1
                gene_id = ':'.join(('NCBIGene', gene_num))
                discontinued_gene_id = ':'.join(('NCBIGene', discontinued_num))
                discontinued_symbol = row[col.index('Discontinued_Symbol')].strip()
                # add the two genes
                if self.class_or_indiv.get(gene_id) == 'C':
                    model.addClassToGraph(gene_id, None)
                    model.addClassToGraph(
                        discontinued_gene_id, discontinued_symbol)

                    # add the new gene id to replace the old gene id
                    model.addDeprecatedClass(discontinued_gene_id, [gene_id])
                else:
                    model.addIndividualToGraph(gene_id, None)
                    model.addIndividualToGraph(
                        discontinued_gene_id, discontinued_symbol)
                    model.addDeprecatedIndividual(discontinued_gene_id, [gene_id])

                # also add the old symbol as a synonym of the new gene
                model.addSynonym(gene_id, discontinued_symbol)

                if not self.test_mode and (limit is not None and line_counter > limit):
                    break

        return
85137
def load_just_in_time_genome_alignment(path, ref_spec, extensions=None,
                                       index_exts=None, fail_no_index=True,
                                       verbose=False):
    """"""
    if index_exts is None and fail_no_index:
      raise ValueError("Failure on no index specified for loading genome " +
                       "alignment, but no index extensions specified")

    extensions = __trim_extensions_dot(extensions)
    index_exts = __trim_extensions_dot(index_exts)

    partial_chrom_files = {}
    whole_chrom_files = {}
    for fn in os.listdir(path):
      pth = os.path.join(path, fn)
      if os.path.isfile(pth):
        base, ext = os.path.splitext(pth)
        ext = __trim_extension_dot(ext)
        if extensions is None or ext in extensions:
          idx_path = __find_index(pth, index_exts)
          if idx_path is None and fail_no_index:
            raise PyokitIOError("No index file for " + fn)
          chrom, start, end = __split_genomic_interval_filename(fn)
          assert((start is None and end is None) or
                 (start is not None and end is not None))
          if start is None:
            if chrom in whole_chrom_files:
              raise PyokitIOError("multiple files for chrom " + chrom)
            whole_chrom_files[chrom] = (pth, idx_path)
          else:
            k = (chrom, start, end)
            if k in partial_chrom_files:
              other_pth = partial_chrom_files[k][0]
              raise PyokitIOError("multiple files for " + str(k) + " --> " +
                                  pth + " and " + other_pth)
            partial_chrom_files[k] = (pth, idx_path)

    def factory(k):
      pth, idx = k
      return build_genome_alignment_from_file(pth, ref_spec, idx, verbose)

    return JustInTimeGenomeAlignment(whole_chrom_files, partial_chrom_files,
                                     factory)