Retrieve larger chromosomes, avoiding the smaller ones for plotting.
545640
def _get_larger_chroms(ref_file):
    
    from scipy.cluster.vq import kmeans, vq
    all_sizes = []
    for c in ref.file_contigs(ref_file):
        all_sizes.append(float(c.size))
    all_sizes.sort()
    if len(all_sizes) > 5:
        # separate out smaller chromosomes and haplotypes with kmeans
        centroids, _ = kmeans(np.array(all_sizes), 2)
        idx, _ = vq(np.array(all_sizes), centroids)
        little_sizes = tz.first(tz.partitionby(lambda xs: xs[0], zip(idx, all_sizes)))
        little_sizes = [x[1] for x in little_sizes]
        # create one more cluster with the smaller, removing the haplotypes
        centroids2, _ = kmeans(np.array(little_sizes), 2)
        idx2, _ = vq(np.array(little_sizes), centroids2)
        little_sizes2 = tz.first(tz.partitionby(lambda xs: xs[0], zip(idx2, little_sizes)))
        little_sizes2 = [x[1] for x in little_sizes2]
        # get any chromosomes not in haplotype/random bin
        thresh = max(little_sizes2)
    else:
        thresh = 0
    larger_chroms = []
    for c in ref.file_contigs(ref_file):
        if c.size > thresh:
            larger_chroms.append(c.name)
    return larger_chroms
1237336
def mutate(self, chromosomes, p_mutate):
        
        assert 0 <= p_mutate <= 1
        
        for chromosome in chromosomes:
            chromosome.mutate(p_mutate)
95639
def effective_genome_size(fasta, read_length, nb_cores, tmpdir="/tmp"):
    # type: (str, int, int, str) -> None
    """"""

    idx = Fasta(fasta)

    genome_length = sum([len(c) for c in idx])

    logging.info("Temporary directory: " + tmpdir)
    logging.info("File analyzed: " + fasta)
    logging.info("Genome length: " + str(genome_length))
    print("File analyzed: ", fasta)
    print("Genome length: ", genome_length)

    chromosomes = ", ".join([c.name for c in idx])

    if "_" in chromosomes:
        print("Warning. The following chromosomes are part of your genome:\n",
              chromosomes.replace(">", "") + "\n",
              file=sys.stderr)
        print(
            "You probably want to remove all chromosomes in your fasta containing '_' for the effective genome size computation to be accurate.",
            file=sys.stderr)

    # if tmpdir is None:
    #     try:
    #         tmpdir = os.environ['TMPDIR']
    #     except KeyError:
    #         tmpdir = '/tmp'
    output_file = os.path.join(tmpdir, '{1}.jf'.format(read_length,
                                                       basename(fasta)))
    atexit.register(
        lambda: call("rm {output_file}".format(output_file=output_file), shell=True))

    call(
        "jellyfish count -t {nb_cores} -m {read_length} -s {genome_length} -L 1 -U 1 --out-counter-len 1 --counter-len 1 {fasta} -o {output_file}".format(
            **vars()),
        shell=True)

    stats = check_output("jellyfish stats {output_file}".format(
        output_file=output_file),
                         shell=True)

    unique_kmers = int(stats.split()[1])

    effective_genome_size = unique_kmers / genome_length

    logging.info("Number unique {read_length}-mers: ".format(
        read_length=read_length) + str(unique_kmers))
    logging.info("Effective genome size: " + str(effective_genome_size))
    print("Number unique {read_length}-mers: ".format(read_length=read_length),
          unique_kmers)
    print("Effective genome size: ", effective_genome_size)
    assert effective_genome_size < 1, "Something wrong happened, effective genome size over 1!"
841006
def get_chromosome_lengths(fasta_file, fancy_sort=True):
    """"""
    chromlen = []
    with gzip.open(fasta_file, 'rt', encoding='ascii') as fh:
        fasta = SeqIO.parse(fh, 'fasta')
        for i, f in enumerate(fasta):
            chromlen.append((f.id, len(f.seq)))
            _LOGGER.info('Processed chromosome "%s"...', f.id)
            #print(dir(f))
            #if i == 1: break
            
    # convert to pandas Series
    chromlen = pd.Series(OrderedDict(chromlen))
    chromlen.index.name = 'Chromosome'
    chromlen.name = 'Length'

    if fancy_sort:
        # sort using fancy ordering
        chrom_for_sorting = chromlen.index.to_series().apply(_transform_chrom)
        a = chrom_for_sorting.argsort(kind='mergesort')
        chromlen = chromlen.iloc[a]

    return chromlen
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)
1257589
def get_gene2section2gos(gene2gos, sec2gos):
        """"""
        gene2section2gos = {}
        for geneid, gos_gene in gene2gos.items():
            section2gos = {}
            for section_name, gos_sec in sec2gos.items():
                gos_secgene = gos_gene.intersection(gos_sec)
                if gos_secgene:
                    section2gos[section_name] = gos_secgene
            gene2section2gos[geneid] = section2gos
        return gene2section2gos
799136
def sample_chromosomes(job, genome_fai_file):
    
    work_dir = os.getcwd()
    genome_fai = untargz(job.fileStore.readGlobalFile(genome_fai_file), work_dir)
    return chromosomes_from_fai(genome_fai)
954107
def _check_bam_contigs(in_bam, ref_file, config):
    
    # GATK allows chromosome M to be in multiple locations, skip checking it
    allowed_outoforder = ["chrM", "MT"]
    ref_contigs = [c.name for c in ref.file_contigs(ref_file, config)]
    with pysam.Samfile(in_bam, "rb") as bamfile:
        bam_contigs = [c["SN"] for c in bamfile.header["SQ"]]
    extra_bcs = [x for x in bam_contigs if x not in ref_contigs]
    extra_rcs = [x for x in ref_contigs if x not in bam_contigs]
    problems = []
    warnings = []
    for bc, rc in zip_longest([x for x in bam_contigs if (x not in extra_bcs and
                                                                     x not in allowed_outoforder)],
                                         [x for x in ref_contigs if (x not in extra_rcs and
                                                                     x not in allowed_outoforder)]):
        if bc != rc:
            if bc and rc:
                problems.append("Reference mismatch. BAM: %s Reference: %s" % (bc, rc))
            elif bc:
                warnings.append("Extra BAM chromosomes: %s" % bc)
            elif rc:
                warnings.append("Extra reference chromosomes: %s" % rc)
    for bc in extra_bcs:
        warnings.append("Extra BAM chromosomes: %s" % bc)
    for rc in extra_rcs:
        warnings.append("Extra reference chromosomes: %s" % rc)
    if problems:
        raise ValueError("Unexpected order, name or contig mismatches between input BAM and reference file:\n%s\n"
                         "Setting `bam_clean: remove_extracontigs` in the configuration can often fix this issue."
                         % "\n".join(problems))
    if warnings:
        print("*** Potential problems in input BAM compared to reference:\n%s\n" %
              "\n".join(warnings))
599422
def normalize(self):
        """"""
        length = self.length()
        if length > 0:
            self.X /= length
            self.Y /= length
        else:
            print "Length 0, cannot normalize."
794753
def batch_scan(self, cells: np.ndarray = None, genes: np.ndarray = None, axis: int = 0, batch_size: int = 1000, layer: str = None) -> Iterable[Tuple[int, np.ndarray, np.ndarray]]:
		
		deprecated("'batch_scan' is deprecated. Use 'scan' instead")
		if cells is None:
			cells = np.fromiter(range(self.shape[1]), dtype='int')
		if genes is None:
			genes = np.fromiter(range(self.shape[0]), dtype='int')
		if layer is None:
			layer = ""
		if axis == 1:
			cols_per_chunk = batch_size
			ix = 0
			while ix < self.shape[1]:
				cols_per_chunk = min(self.shape[1] - ix, cols_per_chunk)

				selection = cells - ix
				# Pick out the cells that are in this batch
				selection = selection[np.where(np.logical_and(selection >= 0, selection < cols_per_chunk))[0]]
				if selection.shape[0] == 0:
					ix += cols_per_chunk
					continue

				# Load the whole chunk from the file, then extract genes and cells using fancy indexing
				vals = self.layers[layer][:, ix:ix + cols_per_chunk]
				vals = vals[genes, :]
				vals = vals[:, selection]

				yield (ix, ix + selection, vals)
				ix += cols_per_chunk

		if axis == 0:
			rows_per_chunk = batch_size
			ix = 0
			while ix < self.shape[0]:
				rows_per_chunk = min(self.shape[0] - ix, rows_per_chunk)

				selection = genes - ix
				# Pick out the genes that are in this batch
				selection = selection[np.where(np.logical_and(selection >= 0, selection < rows_per_chunk))[0]]
				if selection.shape[0] == 0:
					ix += rows_per_chunk
					continue

				# Load the whole chunk from the file, then extract genes and cells using fancy indexing
				vals = self.layers[layer][ix:ix + rows_per_chunk, :]
				vals = vals[selection, :]
				vals = vals[:, cells]
				yield (ix, ix + selection, vals)
				ix += rows_per_chunk
805533
def _get_sv_exclude_file(items):
    
    sv_bed = utils.get_in(items[0], ("genome_resources", "variation", "sv_repeat"))
    if sv_bed and os.path.exists(sv_bed):
        return sv_bed
644293
def include_block(self, x):
        
        last_pos = self._chr_last_blocks.get(x.chrom, 0)
        # Region excludes an entire chromosome, typically decoy/haplotypes
        if last_pos <= self._end_buffer and x.stop >= self._ref_sizes.get(x.chrom, 0) - self._end_buffer:
            return True
        # Do not split on smaller decoy and haplotype chromosomes
        elif self._ref_sizes.get(x.chrom, 0) <= self._target_size:
            return False
        elif (x.start - last_pos) > self._target_size:
            self._chr_last_blocks[x.chrom] = x.stop
            return True
        else:
            return False
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
134319
def _importSequence(chromosome, fastaFile, targetDir) :
    ""

    f = gzip.open(fastaFile)
    header = f.readline()
    strRes = f.read().upper().replace('\n', '').replace('\r', '')
    f.close()

    fn = '%s/chromosome%s.dat' % (targetDir, chromosome.number)
    f = open(fn, 'w')
    f.write(strRes)
    f.close()
    chromosome.dataFile = fn
    chromosome.header = header
    return len(strRes)
270330
def untrim(self, n='all'):
        
        if len(self._set_xdata)==0 or len(self._set_ydata)==0:
            self._error("No data. Please use set_data() and plot() prior to zooming.")
            return

        if   _s.fun.is_a_number(n): n = [n]
        elif isinstance(n,str):     n = list(range(len(self._set_xdata)))

        # loop over the specified plots
        for i in n:
            self['xmin'][i] = None
            self['xmax'][i] = None
            self['ymin'][i] = None
            self['ymax'][i] = None
            
        # now show the update.
        self.clear_results()
        if self['autoplot']: self.plot()
        return self
1327429
def get_variants_in_region(self, chrom, start, end):
        """"""
        if self.chrom is not None and chrom == self.chrom:
            # We are going to search for 'NA' since the chromosome was set
            chrom = "NA"

        iterator = self._bgen.iter_variants_in_region(
            CHROM_STR_DECODE.get(chrom, chrom), start, end,
        )
        for info, dosage in iterator:
            yield Genotypes(
                Variant(
                    info.name, CHROM_STR_ENCODE.get(info.chrom, info.chrom),
                    info.pos, [info.a1, info.a2],
                ),
                dosage,
                reference=info.a1,
                coded=info.a2,
                multiallelic=True,
            )
443155
def convert_from_sliced_object(data):
    """"""
    if data.base is not None and isinstance(data, np.ndarray) and isinstance(data.base, np.ndarray):
        if not data.flags.c_contiguous:
            warnings.warn("Usage of np.ndarray subset (sliced data) is not recommended "
                          "due to it will double the peak memory cost in LightGBM.")
            return np.copy(data)
    return data