Add a genomic range to the locus

    :param grange:
    :type grange: GenomicRange
518153
def add_member(self,grange):
    
    if self.use_direction and not grange.direction:
      sys.stderr.write("ERROR if using direction then direction of input members must be set\n")
      sys.exit()
    # Get range set properly
    if not self.range:
      self.range = GenomicRange(grange.chr,grange.start,grange.end)
      if self.use_direction:
        self.range.set_direction(grange.get_direction())
    elif self.range.chr != grange.chr:
      sys.stderr.write("WARNING cannot add member with chromosomes are not equal\n")
      return False
    elif self.use_direction and self.range.direction != grange.direction:
      sys.stderr.write("WARNING cannot add member with different directions\n")
      return False
    else:
      if grange.start < self.range.start:  self.range.start = grange.start
      if grange.end > self.range.end: self.range.end = grange.end
    self.members.append(grange)
1237336
def mutate(self, chromosomes, p_mutate):
        
        assert 0 <= p_mutate <= 1
        
        for chromosome in chromosomes:
            chromosome.mutate(p_mutate)
253649
def allele_counts_dataframe(variant_and_allele_reads_generator):
    
    df_builder = DataFrameBuilder(
        AlleleCount,
        extra_column_fns={
            "gene": lambda variant, _: ";".join(variant.gene_names),
        })
    for variant, allele_reads in variant_and_allele_reads_generator:
        counts = count_alleles_at_variant_locus(variant, allele_reads)
        df_builder.add(variant, counts)
    return df_builder.to_dataframe()
1335988
def add_range(self, start, stop, step=1, pad=0):
        
        assert start < stop, "please provide ordered node index ranges"
        assert step > 0
        assert pad >= 0
        assert stop - start < 1e9, "range too large"

        if pad > 0 and self.padding is None:
            self.padding = pad
        set.update(self, range(start, stop, step))
1242986
def _add_genes_to_bed(in_file, gene_file, fai_file, out_file, data, max_distance=10000):
    
    try:
        input_rec = next(iter(pybedtools.BedTool(in_file)))
    except StopIteration:  # empty file
        utils.copy_plus(in_file, out_file)
        return
    # keep everything after standard chrom/start/end, 1-based
    extra_fields = list(range(4, len(input_rec.fields) + 1))
    # keep the new gene annotation
    gene_index = len(input_rec.fields) + 4
    extra_fields.append(gene_index)
    columns = ",".join([str(x) for x in extra_fields])
    max_column = max(extra_fields) + 1
    ops = ",".join(["distinct"] * len(extra_fields))
    # swap over gene name to '.' if beyond maximum distance
    # cut removes the last distance column which can cause issues
    # with bedtools merge: 'ERROR: illegal character '.' found in integer conversion of string'
    distance_filter = (r"""awk -F$'\t' -v OFS='\t' '{if ($NF > %s || $NF < -%s) $%s = "."} {print}'""" %
                       (max_distance, max_distance, gene_index))
    sort_cmd = bedutils.get_sort_cmd(os.path.dirname(out_file))
    cat_cmd = "zcat" if in_file.endswith(".gz") else "cat"
    # Ensure gene transcripts match reference genome
    ready_gene_file = os.path.join(os.path.dirname(out_file), "%s-genomeonly.bed" %
                                   (utils.splitext_plus(os.path.basename(gene_file))[0]))
    ready_gene_file = bedutils.subset_to_genome(gene_file, ready_gene_file, data)
    exports = "export TMPDIR=%s && %s" % (os.path.dirname(out_file), utils.local_path_export())
    bcbio_py = sys.executable
    gsort = config_utils.get_program("gsort", data)
    cmd = ("{exports}{cat_cmd} {in_file} | grep -v ^track | grep -v ^browser | grep -v ^# | "
           "{bcbio_py} -c 'from bcbio.variation import bedutils; bedutils.remove_bad()' | "
           "{gsort} - {fai_file} | "
            "bedtools closest -g {fai_file} "
            "-D ref -t first -a - -b <({gsort} {ready_gene_file} {fai_file}) | "
            "{distance_filter} | cut -f 1-{max_column} | "
            "bedtools merge -i - -c {columns} -o {ops} -delim ',' -d -10 > {out_file}")
    do.run(cmd.format(**locals()), "Annotate BED file with gene info")
248808
def get_piles(allgenes):
    
    from jcvi.utils.range import Range, range_piles

    ranges = [Range(a.seqid, a.start, a.end, 0, i) \
                    for i, a in enumerate(allgenes)]

    for pile in range_piles(ranges):
        yield [allgenes[x] for x in pile]
629258
def copy(self):
        '''

        '''
        return type(self)(self.name, self.start, self.stop, self.feature_type,
                          gene=self.gene, locus_tag=self.locus_tag,
                          qualifiers=self.qualifiers, strand=self.strand)
921720
def _apply_ranges(self, pat, hosts):
         

        (loose_pattern, limits) = self._enumeration_info(pat)
        if not limits:
            return hosts

        (left, right) = limits
        enumerated = enumerate(hosts)
        if left == '':
            left = 0
        if right == '':
            right = 0
        left=int(left)
        right=int(right)
        enumerated = [ h for (i,h) in enumerated if i>=left and i<=right ]
        return enumerated
676387
def g_to_n(self, g_interval):
        """"""

        grs, gre = g_interval.start.base - 1 - self.gc_offset, g_interval.end.base - 1 - self.gc_offset
        # frs, fre = (f)orward (r)na (s)tart & (e)nd; forward w.r.t. genome
        frs, frs_offset, frs_cigar = self._map(
            from_pos=self.ref_pos, to_pos=self.tgt_pos, pos=grs, base="start")
        fre, fre_offset, fre_cigar = self._map(
            from_pos=self.ref_pos, to_pos=self.tgt_pos, pos=gre, base="end")

        if self.strand == -1:
            frs, fre = self.tgt_len - fre - 1, self.tgt_len - frs - 1
            frs_offset, fre_offset = -fre_offset, -frs_offset

        # The returned interval would be uncertain when locating at alignment gaps
        return hgvs.location.BaseOffsetInterval(
            start=hgvs.location.BaseOffsetPosition(
                base=frs + 1, offset=frs_offset, datum=Datum.SEQ_START),
            end=hgvs.location.BaseOffsetPosition(
                base=fre + 1, offset=fre_offset, datum=Datum.SEQ_START),
            uncertain=frs_cigar in 'DI' or fre_cigar in 'DI')
626158
def add(self, ranges):
        """
        if len(ranges) and isinstance(ranges[0], int_types):
            ranges = [ranges]
        iset = self._iset
        self._maxlen = max(self._maxlen, max(r[1] - r[0] + 1 for r in ranges))

        if len(ranges) > 30 or len(iset) < len(ranges):
            iset.extend(ranges)
            iset.sort()
        else:
            for o in ranges:
                iset.insert(binsearch_left_start(iset, o[0], 0, len(iset)), o)
942219
def _add_genetic_models(self, variant_obj, info_dict):
        
        genetic_models_entry = info_dict.get('GeneticModels')
        if genetic_models_entry:
            genetic_models = []
            for family_annotation in genetic_models_entry.split(','):
                for genetic_model in family_annotation.split(':')[-1].split('|'):
                    genetic_models.append(genetic_model)
            logger.debug("Updating genetic models to: {0}".format(
                ', '.join(genetic_models)))
                
            variant_obj.genetic_models = genetic_models
1019750
def apply(self, read_tuple_name, read_tuple_id=None, synchronize_widths=True):
        
        parts = read_tuple_name.split("__")
        parts[0] = self._fill_right(parts[0], "-", self.prefix_width)
        if read_tuple_id is not None:
            parts[1] = "{:x}".format(read_tuple_id)
        parts[1] = self._fill_left(parts[1], "0", self.read_tuple_id_width)

        if synchronize_widths:
            new_segments = []
            segments = parts[2][1:-1].split("),(")
            for segment in segments:
                values = segment.split(",")
                values[0] = values[0].zfill(self.genome_id_width)
                values[1] = values[1].zfill(self.chr_id_width)
                values[3] = values[3].zfill(self.coor_width)
                values[4] = values[4].zfill(self.coor_width)
                new_segments.append("(" + ",".join(values) + ")")
            parts[2] = ",".join(new_segments)

        return "__".join(parts)
1125416
def mstmap(args):
    
    from jcvi.assembly.geneticmap import MSTMatrix

    p = OptionParser(mstmap.__doc__)
    p.add_option("--population_type", default="RIL6",
                 help="Type of population, possible values are DH and RILd")
    p.add_option("--missing_threshold", default=.5,
                 help="Missing threshold, .25 excludes any marker with >25% missing")
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    lmd, = args
    fp = open(lmd)
    next(fp)  # Header
    table = {"0": "-", "1": "A", "2": "B", "3": "X"}
    mh = ["locus_name"] + fp.next().split()[4:]
    genotypes = []
    for row in fp:
        atoms = row.split()
        chr, pos, ref, alt = atoms[:4]
        locus_name = ".".join((chr, pos))
        codes = [table[x] for x in atoms[4:]]
        genotypes.append([locus_name] + codes)

    mm = MSTMatrix(genotypes, mh, opts.population_type, opts.missing_threshold)
    mm.write(opts.outfile, header=True)
1180970
public void visit(BioPAXElement domain, Object range, Model targetModel, PropertyEditor editor)
	{
		BioPAXElement targetDomain = targetModel.getByID(domain.getUri());
		
		if (range instanceof BioPAXElement)
		{
			BioPAXElement bpe = (BioPAXElement) range;
			BioPAXElement existing = targetModel.getByID(bpe.getUri());
			//set the property value if the value is already present in the target 
			if (existing != null) {
				editor.setValueToBean(existing, targetDomain);
			}
		}
		else
		{
			editor.setValueToBean(range, targetDomain);
		}
	}
137107
def _add_gene_disease(self, row):  # ::List  getting syntax error here
        
        col = self.files['developmental_disorders']['columns']
        if len(row) != len(col):
            raise ValueError("Unexpected number of fields for row {}".format(row))

        variant_label = "variant of {}".format(row[col.index('gene_symbol')])
        disease_omim_id = row[col.index('disease_omim_id')]
        if disease_omim_id == 'No disease mim':
            # check if we've manually curated
            disease_label = row[col.index('disease_label')]
            if disease_label in self.mondo_map:
                disease_id = self.mondo_map[disease_label]
            else:
                return  # sorry for this
        else:
            disease_id = 'OMIM:' + disease_omim_id

        hgnc_curie = 'HGNC:' + row[col.index('hgnc_id')]

        relation_curie = self.resolve(row[col.index('g2p_relation_label')])
        mutation_consequence = row[col.index('mutation_consequence')]
        if mutation_consequence not in ('uncertain', ''):
            consequence_relation = self.resolve(
                self._get_consequence_predicate(mutation_consequence))
            consequence_curie = self.resolve(mutation_consequence)
            variant_label = "{} {}".format(mutation_consequence, variant_label)
        else:
            consequence_relation = None
            consequence_curie = None

        allelic_requirement = row[col.index('allelic_requirement')]
        if allelic_requirement != '':
            requirement_curie = self.resolve(allelic_requirement)
        else:
            requirement_curie = None

        pmids = row[col.index('pmids')]
        if pmids != '':
            pmid_list = ['PMID:' + pmid for pmid in pmids.split(';')]
        else:
            pmid_list = []

        # build the model
        # Should we build a reusable object and/or tuple that
        # could be passed to a more general model builder for
        # this and orphanet (and maybe clinvar)
        self._build_gene_disease_model(
            hgnc_curie,
            relation_curie,
            disease_id,
            variant_label,
            consequence_relation,
            consequence_curie,
            requirement_curie,
            pmid_list
        )
549012
def _process_trait_mappings(self, raw, limit=None):
        
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        line_counter = 0
        model = Model(graph)

        with open(raw, 'r') as csvfile:
            filereader = csv.reader(csvfile, delimiter=',', quotechar='\"')
            next(filereader, None)  # skip header line
            for row in filereader:
                line_counter += 1
                # need to skip the last line
                if len(row) < 8:
                    LOG.info("skipping line %d: %s", line_counter, '\t'.join(row))
                    continue
                (vto_id, pto_id, cmo_id, ato_column, species, trait_class,
                 trait_type, qtl_count) = row

                ato_id = re.sub(
                    r'ATO #', 'AQTLTrait:', re.sub(
                        r'\].*', '', re.sub(r'\[', '', ato_column)))
                ato_id = ato_id.strip()

                ato_label = re.sub(r'.*\]\s*', '', ato_column)

                model.addClassToGraph(ato_id, ato_label.strip())

                if re.match(r'VT:.*', vto_id):
                    model.addClassToGraph(vto_id, None)
                    model.addEquivalentClass(ato_id, vto_id)
                if re.match(r'LPT:.*', pto_id):
                    model.addClassToGraph(pto_id, None)
                    model.addXref(ato_id, pto_id)
                if re.match(r'CMO:.*', cmo_id):
                    model.addClassToGraph(cmo_id, None)
                    model.addXref(ato_id, cmo_id)

        LOG.info("Done with trait mappings")
        return
439169
public void addTransientRange(String id, RowRange range) {
    String start = DatatypeConverter.printHexBinary(range.getStart().toArray());
    String end = DatatypeConverter.printHexBinary(range.getEnd().toArray());

    appConfig.setProperty(PREFIX + id, start + ":" + end);
  }