Use intersectBed to find the fused gene => split genes mappings.
490619
def get_splits(split_bed, gff_file, stype, key):
    
    bed_file = get_bed_file(gff_file, stype, key)
    cmd = "intersectBed -a {0} -b {1} -wao".format(split_bed, bed_file)
    cmd += " | cut -f4,10"
    p = popen(cmd)
    splits = defaultdict(set)
    for row in p:
        a, b = row.split()
        splits[a].add(b)

    return splits
241545
def subset_to_genome(in_file, out_file, data):
    
    if not utils.file_uptodate(out_file, in_file):
        contigs = set([x.name for x in ref.file_contigs(dd.get_ref_file(data))])
        with utils.open_gzipsafe(in_file) as in_handle:
            with file_transaction(data, out_file) as tx_out_file:
                with open(tx_out_file, "w") as out_handle:
                    for line in in_handle:
                        parts = line.split()
                        if parts and parts[0] in contigs:
                            out_handle.write(line)
    return out_file
929681
private int addBAISplits(List<InputSplit> splits,
                                 int i,
                                 List<InputSplit> newSplits,
                                 Configuration conf) throws IOException {
                int splitsEnd = i;

                final Path path = ((FileSplit)splits.get(i)).getPath();
                final Path baiPath = getBAIPath(path);
                final FileSystem fs = path.getFileSystem(conf);
                final Path sinPath;
                if (fs.exists(baiPath)) {
                    sinPath = baiPath;
                } else {
                    sinPath = new Path(path.toString().replaceFirst("\\.bam$", BAMIndex.BAMIndexSuffix));
                }
                try (final FSDataInputStream in = fs.open(path);
                     final SeekableStream guesserSin = WrapSeekable.openPath(fs, path);
                     final SeekableStream sin = WrapSeekable.openPath(fs, sinPath)) {

                        SAMFileHeader header = SAMHeaderReader.readSAMHeaderFrom(in, conf);
                        SAMSequenceDictionary dict = header.getSequenceDictionary();
                        
                        final BAMSplitGuesser guesser = new BAMSplitGuesser(guesserSin, conf);

                        final LinearBAMIndex idx = new LinearBAMIndex(sin, dict);

			// searches for the first contig that contains linear bins
			// a contig will have no linear bins if there are no reads mapped to that
			// contig (e.g., reads were aligned to a whole genome, and then reads from
			// only a single contig were selected)
                        int ctgIdx = -1;
                        int bin = 0;
                        LinearIndex linIdx;
                        int ctgBins;
                        long lastStart = 0;
                        do {
                                ctgIdx++;
                                linIdx = idx.getLinearIndex(ctgIdx);
                                ctgBins = linIdx.size();
                        } while(ctgBins == 0);
                        long nextStart = linIdx.get(bin);
                        
                        FileVirtualSplit newSplit = null;
                        boolean lastWasGuessed = false;

			// loop and process all of the splits that share a single .bai
                        while(splitsEnd < splits.size() &&
                              ((FileSplit)(splits.get(splitsEnd))).getPath() == path) {
                                FileSplit fSplit = (FileSplit)splits.get(splitsEnd);
                                splitsEnd++;
                                
                                if (splitsEnd >= splits.size()) {
                                        break;
                                }

                                long fSplitEnd = (fSplit.getStart() + fSplit.getLength()) << 16;
                                lastStart = nextStart;

				// we need to advance and find the first linear index bin
				// that starts after the current split ends.
				// this is the end of our split.
				while(nextStart < fSplitEnd && ctgIdx < dict.size()) {

					// are we going off of the end of this contig?
					// if so, advance to the next contig with a linear bin
					if (bin + 1 >= ctgBins) {
						do {
                                                        ctgIdx += 1;
                                                        bin = 0;
                                                        if (ctgIdx >= dict.size()) {
                                                                break;
                                                        }
                                                        linIdx = idx.getLinearIndex(ctgIdx);
                                                        ctgBins = linIdx.size();
                                                } while (ctgBins == 0);
                                        }
                                        if (ctgIdx < dict.size() && linIdx.size() > bin) {
                                                nextStart = linIdx.get(bin);
                                                bin++;
                                        }
                                }

				// is this the first split?
				// if so, split ranges from where the reads start until the identified end
                                if (fSplit.getStart() == 0) {
                                        try (final SeekableStream inFile = WrapSeekable.openPath(path.getFileSystem(conf), path)) {
                                            SamReader open = SamReaderFactory.makeDefault().setUseAsyncIo(false)
                                                .open(SamInputResource.of(inFile));
                                            SAMFileSpan span = open.indexing().getFilePointerSpanningReads();
                                            long bamStart = ((BAMFileSpan) span).getFirstOffset();
                                            newSplit = new FileVirtualSplit(fSplit.getPath(),
                                                                        bamStart,
                                                                        nextStart - 1,
                                                                        fSplit.getLocations());
                                            newSplits.add(newSplit);
                                        }
                                } else {

					// did we find any blocks that started in the last split?
					// if yes, then we're fine
					// if no, then we need to guess a split start (in the else clause)
					if (lastStart != nextStart) {
                                                if (lastWasGuessed) {
                                                        newSplit.setEndVirtualOffset(lastStart - 1);
                                                        lastWasGuessed = false;
                                                }
                                                newSplit = new FileVirtualSplit(fSplit.getPath(),
                                                                                lastStart,
                                                                                nextStart - 1,
                                                                                fSplit.getLocations());
                                                newSplits.add(newSplit);
                                        } else {
						// guess the start
                                                long alignedBeg = guesser.guessNextBAMRecordStart(fSplit.getStart(),
                                                                                                  fSplit.getStart() + fSplit.getLength());
                                                newSplit.setEndVirtualOffset(alignedBeg - 1);
                                                lastStart = alignedBeg;
                                                nextStart = alignedBeg;
                                                newSplit = new FileVirtualSplit(fSplit.getPath(),
                                                                                alignedBeg,
                                                                                alignedBeg + 1,
                                                                                fSplit.getLocations());
                                                lastWasGuessed = true;
                                                newSplits.add(newSplit);
                                        }
                                }
                                lastStart = nextStart;
                        }
			// clean up the last split
                        if (splitsEnd == splits.size()) {
                                if (lastWasGuessed) {
                                        newSplit.setEndVirtualOffset(lastStart - 1);
                                        lastWasGuessed = false;
                                }
                                FileSplit fSplit = (FileSplit)splits.get(splitsEnd - 1);
                                long fSplitEnd = (fSplit.getStart() + fSplit.getLength()) << 16;
                                newSplit = new FileVirtualSplit(fSplit.getPath(),
                                                                lastStart,
                                                                fSplitEnd,
                                                                fSplit.getLocations());
                                newSplits.add(newSplit);
                        }
                }
                return splitsEnd + 1;
        }
691158
def get_splicejunction_file(align_dir, data):
    
    samplename = dd.get_sample_name(data)
    align_dir = os.path.dirname(dd.get_work_bam(data))
    knownfile = get_known_splicesites_file(align_dir, data)
    novelfile = os.path.join(align_dir, "%s-novelsplicesites.bed" % samplename)
    bed_files = [x for x in [knownfile, novelfile] if file_exists(x)]
    splicejunction = bed.concat(bed_files)
    splicejunctionfile = os.path.join(align_dir,
                                      "%s-splicejunctions.bed" % samplename)
    if splicejunction:
        splicejunction.saveas(splicejunctionfile)
        return splicejunctionfile
    else:
        return None
149479
def standardize_genes(self, inplace=False):
        """"""
        matrix = self.center_genes(inplace=inplace)
        matrix.X[:,:] = matrix.X / \
            np.tile(np.std(matrix.X, axis=1, ddof=1), (matrix.n, 1)).T
        return matrix
1231139
def add_genes(in_file, data, max_distance=10000, work_dir=None):
    
    gene_file = regions.get_sv_bed(data, "exons", out_dir=os.path.dirname(in_file))
    if gene_file and utils.file_exists(in_file):
        out_file = "%s-annotated.bed" % utils.splitext_plus(in_file)[0]
        if work_dir:
            out_file = os.path.join(work_dir, os.path.basename(out_file))
        if not utils.file_uptodate(out_file, in_file):
            fai_file = ref.fasta_idx(dd.get_ref_file(data))
            with file_transaction(data, out_file) as tx_out_file:
                _add_genes_to_bed(in_file, gene_file, fai_file, tx_out_file, data, max_distance)
        return out_file
    else:
        return in_file
1109604
def _whole_gene_covered_by_nucmer_hits(nucmer_hits, ref_seq, percent_threshold, max_nt_extend):
        ''''''
        coords = AssemblyCompare.nucmer_hits_to_ref_coords(nucmer_hits)
        covered = []
        for coords_list in coords.values():
            covered.extend(coords_list)
        pyfastaq.intervals.merge_overlapping_in_list(covered)
        return (2 * max_nt_extend + pyfastaq.intervals.length_sum_from_list(covered)) / len(ref_seq) >= percent_threshold
1318662
def match_subfeats(f1, f2, dbx1, dbx2, featuretype=None, slop=False):
    
    f1c, f2c = list(dbx1.children(f1, featuretype=featuretype, order_by='start')), \
            list(dbx2.children(f2, featuretype=featuretype, order_by='start'))

    lf1c, lf2c = len(f1c), len(f2c)
    if match_nchildren(f1c, f2c):
        if lf1c > 0 and lf2c > 0:
            exclN = set()
            if featuretype.endswith('UTR') or featuretype == 'exon':
                N = []
                if featuretype.startswith('five_prime'):
                    N = [1] if f1.strand == "+" else [lf1c]
                elif featuretype.startswith('three_prime'):
                    N = [lf1c] if f1.strand == "+" else [1]
                else:   # infer UTR from exon collection
                    N = [1] if 1 == lf1c else [1, lf1c]

                for n in N:
                    if match_Nth_child(f1c, f2c, N=n, slop=slop):
                        exclN.add(n-1)
                    else:
                        return False

            for i, (cf1, cf2) in enumerate(zip(f1c, f2c)):
                if i in exclN: continue
                if not match_span(cf1, cf2):
                    return False
    else:
        if (lf1c, lf2c) in [(0, 1), (1, 0)] and slop \
                and featuretype.endswith('UTR'):
            return True

        return False

    return True
112257
def fit_cluster(data, assignments, k, P_init, R_init, means):
    
    for c in range(k):
        if data[:,assignments==c].shape[1] == 0:
            _, assignments = kmeans_pp(data, k)
    genes, cells = data.shape
    nb_gene_indices = np.array([True for i in range(genes)])
    for c in range(k):
        c_data = data[:,assignments==c]
        nb_gene_indices = nb_gene_indices & find_nb_genes(c_data)
    for c in range(k):
        c_data = data[:,assignments==c]
        nb_genes = c_data[nb_gene_indices,:]
        poisson_genes = c_data[~nb_gene_indices, :]
        P_init[nb_gene_indices, c], R_init[nb_gene_indices, c] = nb_fit(nb_genes)
        means[~nb_gene_indices, c] = poisson_genes.mean(1)
    return nb_gene_indices
908881
def merge(args):
    
    from jcvi.formats.base import DictFile

    p = OptionParser(merge.__doc__)
    opts, args = p.parse_args(args)

    if len(args) != 3:
        sys.exit(not p.print_help())

    quartets, registry, lost = args
    qq = DictFile(registry, keypos=1, valuepos=3)
    lost = DictFile(lost, keypos=1, valuepos=0, delimiter='|')
    qq.update(lost)
    fp = open(quartets)
    cases = {
        "AN,CN": 4,
        "BO,AN,CN": 8,
        "BO,CN": 2,
        "BR,AN": 1,
        "BR,AN,CN": 6,
        "BR,BO": 3,
        "BR,BO,AN": 5,
        "BR,BO,AN,CN": 9,
        "BR,BO,CN": 7,
    }
    ip = {
        "syntenic_model": "Syntenic_model_excluded_by_OMG",
        "complete": "Predictable",
        "partial": "Truncated",
        "pseudogene": "Pseudogene",
        "random": "Match_random",
        "real_ns": "Transposed",
        "gmap_fail": "GMAP_fail",
        "AN LOST": "AN_LOST",
        "CN LOST": "CN_LOST",
        "BR LOST": "BR_LOST",
        "BO LOST": "BO_LOST",
        "outside": "Outside_synteny_blocks",
        "[NF]": "Not_found",
    }
    for row in fp:
        atoms = row.strip().split("\t")
        genes = atoms[:4]
        tag = atoms[4]
        a, b, c, d = [qq.get(x, ".").rsplit("-", 1)[-1] for x in genes]
        qqs = [c, d, a, b]
        for i, q in enumerate(qqs):
            if atoms[i] != '.':
                qqs[i] = "syntenic_model"
        # Make comment
        comment = "Case{0}".format(cases[tag])
        dots = sum([1 for x in genes if x == '.'])
        if dots == 1:
            idx = genes.index(".")
            status = qqs[idx]
            status = ip[status]
            comment += "-" + status
        print(row.strip() + "\t" + "\t".join(qqs + [comment]))
439319
def center_genes(self, use_median=False, inplace=False):
        """"""
        if use_median:
            X = self.X - \
                np.tile(np.median(self.X, axis=1), (self.n, 1)).T
        else:
            X = self.X - \
                np.tile(np.mean(self.X, axis=1), (self.n, 1)).T

        if inplace:
            self.X[:,:] = X
            matrix = self
        else:
            matrix = ExpMatrix(genes=self.genes, samples=self.samples,
                               X=X)
        return matrix
479847
def get_gene2gos(go2nt):
        """"""
        gene2gos = cx.defaultdict(set)
        nt0 = next(iter(go2nt.values()))
        b_str = isinstance(nt0.study_items, str)
        # print("NNNNTTTT", nt0)
        for goid, ntgo in go2nt.items():
            study_items = ntgo.study_items.split(', ') if b_str else ntgo.study_items
            for geneid in study_items:
                gene2gos[geneid].add(goid)
        if b_str:
          b_set = set(isinstance(g.isdigit(), int) for g in nt0.study_items.split(', '))
          if b_set == set([True]):
            return {int(g):gos for g, gos in gene2gos.items()}
        return {g:gos for g, gos in gene2gos.items()}
26359
def calculate_new_length(gene_split, gene_results, hit):
   ''' 
   '''
   # Looping over splitted hits and calculate new length
   first = 1
   for split in gene_split[hit['sbjct_header']]:
      new_start = int(gene_results[split]['sbjct_start'])
      new_end = int(gene_results[split]['sbjct_end'])
      
      # Get the frist HSP
      if first == 1:
         new_length = int(gene_results[split]['HSP_length'])
         old_start = new_start
         old_end = new_end
         first = 0
         continue
      if new_start < old_start:
         new_length = new_length + (old_start - new_start)
         old_start = new_start
      
      if new_end > old_end:
         new_length = new_length + (new_end - old_end)
         old_end = new_end
   
   return(new_length)
415475
def ngram_intersection(args, parser):
    """"""
    store = utils.get_data_store(args)
    corpus = utils.get_corpus(args)
    catalogue = utils.get_catalogue(args)
    store.validate(corpus, catalogue)
    store.intersection(catalogue, sys.stdout)
357778
def parse(self,tolerance=0,downsample=None,evidence=2,use_gene_names=False):
      """"""
      g = Graph()
      nodes = [Node(x) for x in self._transcripts]
      for n in nodes: g.add_node(n)
      for i in range(0,len(nodes)):
         for j in range(0,len(nodes)):
            if i == j: continue
            jov = nodes[i].payload.junction_overlap(nodes[j].payload,tolerance)
            sub = jov.is_subset()
            if not sub: continue
            if sub == 1:
               g.add_edge(Edge(nodes[i],nodes[j]))
               g.add_edge(Edge(nodes[j],nodes[i]))
            if sub == 2:
               g.add_edge(Edge(nodes[i],nodes[j]))
      g.merge_cycles()
      roots = g.roots
      groups = []
      for r in roots:
         g2 = g.get_root_graph(r)
         c = CompatibleGraph(g2,tolerance,downsample,evidence,use_gene_names=use_gene_names)
         groups.append(c)
      return groups
884803
def fuse(args):
    
    from jcvi.algorithms.graph import BiGraph

    p = OptionParser(fuse.__doc__)
    opts, args = p.parse_args(args)

    if len(args) < 1:
        sys.exit(not p.print_help())

    bedfiles = [x for x in args if x.endswith(".bed")]
    anchorfiles = [x for x in args if x.endswith(".anchors")]

    # TODO: Use Markov clustering to sparsify the edges
    families = Grouper()
    for anchorfile in anchorfiles:
        af = AnchorFile(anchorfile)
        for a, b, block_id in af.iter_pairs():
            families.join(a, b)

    allowed = set(families.keys())
    logging.debug("Total families: {}, Gene members: {}"
                  .format(len(families), len(allowed)))

    # TODO: Use C++ implementation of BiGraph() when available
    # For now just serialize this to the disk
    G = BiGraph()
    for bedfile in bedfiles:
        bed = Bed(bedfile, include=allowed)
        #add_bed_to_graph(G, bed, families)
        print_edges(G, bed, families)
954657
def _get_gene2pubmed(self, limit):
        
        src_key = 'gene2pubmed'
        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)
        assoc_counter = 0
        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:
                line_counter += 1
                # skip comments
                row = line.decode().strip().split('\t')
                if row[0][0] == '#':
                    continue
                # (tax_num, gene_num, pubmed_num) = line.split('\t')

                # ## set id_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()
                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

                pubmed_num = row[col.index('PubMed_ID')].strip()
                if gene_num == '-' or pubmed_num == '-':
                    continue

                gene_id = ':'.join(('NCBIGene', gene_num))
                pubmed_id = ':'.join(('PMID', pubmed_num))

                if self.class_or_indiv.get(gene_id) == 'C':
                    model.addClassToGraph(gene_id, None)
                else:
                    model.addIndividualToGraph(gene_id, None)
                # add the publication as a NamedIndividual
                # add type publication
                model.addIndividualToGraph(pubmed_id, None, None)
                reference = Reference(
                    graph, pubmed_id, self.globaltt['journal article'])
                reference.addRefToGraph()
                graph.addTriple(
                    pubmed_id, self.globaltt['is_about'], gene_id)
                assoc_counter += 1
                if not self.test_mode and limit is not None and line_counter > limit:
                    break

        LOG.info(
            "Processed %d pub-gene associations", assoc_counter)

        return