Parse pedigree line and remove excluded individuals from geno

        Translates alleles into numerical genotypes (0, 1, 2) counting
        number of minor alleles.

        Throws exceptions if an there are not 2 distinct alleles
482974
def process_genotypes(self, data):
        
        # Get a list of uniq entries in the data, except for missing
        alleles = list(set(data[4:]) - set(DataParser.missing_representation))
        if len(alleles) > 2:
            raise TooManyAlleles(chr=self.chr, rsid=self.rsid, alleles=alleles)

        # We don't have a way to know this in advance, so we want to just iterate onward
        # if we encounter one of these
        if len(alleles) == 1:
            raise TooFewAlleles(chr=self.chr, rsid=self.rsid, alleles=alleles)

        # Strip out any excluded individuals
        allelic_data = numpy.ma.MaskedArray(numpy.array(data[4:], dtype="S2"), self.ind_mask).compressed().reshape(-1, 2)

        maj_allele_count = numpy.sum(allelic_data==alleles[0])

        min_allele_count = numpy.sum(allelic_data==alleles[1])

        effect_allele_count = min_allele_count
        if min_allele_count > maj_allele_count:
            alleles = [alleles[1], alleles[0]]
            allele_count = maj_allele_count
            maj_allele_count = min_allele_count
            min_allele_count = allele_count

        #genotypes = []
        major_allele       = alleles[0]
        minor_allele       = alleles[1]

        # Genotypes represent the sum of minor alleles at each sample
        genotype_data = numpy.sum(allelic_data==minor_allele, axis=1)
        missing_alleles = allelic_data[:, 0]==DataParser.missing_representation
        genotype_data[missing_alleles] = DataParser.missing_storage
        hetero_count = numpy.sum(genotype_data==1)

        return (genotype_data,
                major_allele,
                minor_allele,
                hetero_count,
                maj_allele_count,
                min_allele_count,
                numpy.sum(missing_alleles),
                effect_allele_count)
923295
def fix_vcf_line(parts, ref_base):
    
    swap = {"1/1": "0/0", "0/1": "0/1", "0/0": "1/1", "./.": "./."}
    complements = {"G": "C", "A": "T", "C": "G", "T": "A", "N": "N"}
    varinfo, genotypes = fix_line_problems(parts)
    ref, var = varinfo[3:5]
    # non-reference regions or non-informative, can't do anything
    if ref_base in [None, "N"] or set(genotypes) == set(["./."]):
        varinfo = None
    # matching reference, all good
    elif ref_base == ref:
        assert ref_base == ref, (ref_base, parts)
    # swapped reference and alternate regions
    elif ref_base == var or ref in ["N", "0"]:
        varinfo[3] = var
        varinfo[4] = ref
        genotypes = [swap[x] for x in genotypes]
    # reference is on alternate strand
    elif ref_base != ref and complements.get(ref) == ref_base:
        varinfo[3] = complements[ref]
        varinfo[4] = ",".join([complements[v] for v in var.split(",")])
    # unspecified alternative base
    elif ref_base != ref and var in ["N", "0"]:
        varinfo[3] = ref_base
        varinfo[4] = ref
        genotypes = [swap[x] for x in genotypes]
    # swapped and on alternate strand
    elif ref_base != ref and complements.get(var) == ref_base:
        varinfo[3] = complements[var]
        varinfo[4] = ",".join([complements[v] for v in ref.split(",")])
        genotypes = [swap[x] for x in genotypes]
    else:
        print "Did not associate ref {0} with line: {1}".format(
            ref_base, varinfo)
    if varinfo is not None:
        return varinfo + genotypes
242229
def populate_iteration(self, iteration):
        """"""

        cur_idx = iteration.cur_idx
        genotypes = self.genotype_file.next().split()
        iteration.chr, iteration.rsid, junk, iteration.pos = genotypes[0:4]
        iteration.chr = int(iteration.chr)
        iteration.pos = int(iteration.pos)

        if DataParser.boundary.TestBoundary(iteration.chr, iteration.pos, iteration.rsid):
            try:
                [iteration.genotype_data,
                    iteration.major_allele,
                    iteration.minor_allele,
                    iteration.hetero_count,
                    iteration.maj_allele_count,
                    iteration.min_allele_count,
                    iteration.missing_allele_count,
                    iteration.allele_count2] = self.process_genotypes(genotypes)
                return iteration.maf >= DataParser.min_maf and iteration.maf <= DataParser.max_maf
            except TooFewAlleles:
                print "\n\n\nSkipping %s:%s %s %s" % (iteration.chr, iteration.pos, iteration.rsid, cur_idx)

        return False
474440
def normalize_alleles_left(ref, start, stop, alleles, bound, ref_step, shuffle=True):
    

    normalized_alleles = namedtuple('shuffled_alleles', 'start stop alleles')

    if len(alleles) < 2:
        return normalized_alleles(start, stop, alleles)

    # STEP 1: Trim common suffix
    trimmed, alleles = trim_common_suffixes(alleles)
    stop -= trimmed

    # STEP 2: Trim common prefix
    trimmed, alleles = trim_common_prefixes(alleles)
    start += trimmed

    # assert bound <= start,'start={:d}, left bound={:d}'.format(start, bound)

    # STEP 3: While a null allele exists, left shuffle by prepending alleles
    #         with reference and trimming common suffixes
    while shuffle and '' in alleles and start > bound:
        step = min(ref_step, start - bound)

        r = ref[start - step:start].upper()
        new_alleles = [r + a for a in alleles]

        trimmed, new_alleles = trim_common_suffixes(new_alleles)

        if not trimmed:
            break

        start -= trimmed
        stop -= trimmed

        if trimmed == step:
            alleles = new_alleles
        else:
            left = step - trimmed
            alleles = [a[left:] for a in new_alleles]
            break

    return normalized_alleles(start, stop, tuple(alleles))
958309
def iter_genotypes(self):
        
        # Iterating over all markers
        for i, (_, genotypes) in enumerate(self.bed.iter_geno()):
            info = self.bim.iloc[i, :]

            yield Genotypes(
                Variant(info.name, CHROM_INT_TO_STR[info.chrom],
                        info.pos, [info.a1, info.a2]),
                self._normalize_missing(genotypes),
                reference=info.a2,
                coded=info.a1,
                multiallelic=info.multiallelic
            )
747531
def discrepant_snps(self):
        
        df = self._discrepant_positions.append(self._discrepant_genotypes)
        if len(df) > 1:
            df = df.drop_duplicates()
        return df
1311204
def get_normal_sample(in_file):
    
    with utils.open_gzipsafe(in_file) as in_handle:
        for line in in_handle:
            if line.startswith("##PEDIGREE"):
                parts = line.strip().split("Original=")[1][:-1]
                return parts
377458
def parse_phlat_file(phlatfile, mhc_alleles):
    
    for line in phlatfile:
        if line.startswith('Locus'):
            continue
        line = line.strip().split()
        if line[0].startswith('HLA_D'):
            line[0] = line[0][:-1]  # strip the last character
        # Sometimes we get an error saying there was insufficient read
        # converage. We need to drop that line.
        # E.g. HLA_DQB1 no call due to insufficient reads at this locus
        if line[1] == 'no':
            continue
        if line[4] != 'NA':
            split_field = line[1].split(':')
            if len(split_field) >= 2 and not split_field[1] == 'xx':
                mhc_alleles[line[0]].append((line[1], line[4]))
        if line[5] != 'NA':
            split_field = line[2].split(':')
            if len(split_field) >= 2 and not split_field[1] == 'xx':
                mhc_alleles[line[0]].append((line[2], line[5]))
    return mhc_alleles
149324
private boolean parseVariant(SequenceLocation sequenceLocation, String line) {
        boolean validVariant;
        String[] fields = line.split("\t", -1);
        String mutationCds = fields[17];

        if (mutationCds.contains(">")) {
            validVariant = parseSnv(mutationCds, sequenceLocation);
            if (!validVariant) {
                invalidSubstitutionLines++;
            }
        } else if (mutationCds.contains("del")) {
            validVariant = parseDeletion(mutationCds, sequenceLocation);
            if (!validVariant) {
                invalidDeletionLines++;
            }
        } else if (mutationCds.contains("ins")) {
            validVariant = parseInsertion(mutationCds, sequenceLocation);
            if (!validVariant) {
                invalidInsertionLines++;
            }
        } else if (mutationCds.contains("dup")) {
            validVariant = parseDuplication(mutationCds);
            if (!validVariant) {
                invalidDuplicationLines++;
            }
        } else {
            validVariant = false;
            invalidMutationCDSOtherReason++;
        }

        return validVariant;
    }
136325
def _splitGenoGeneWindow(self,annotation_file=None,cis=1e4,funct='protein_coding',minSnps=1.,maxSnps=SP.inf):
        
         #1. load annotation
        assert annotation_file is not None, 'Splitter:: specify annotation file'
        try:
            f = h5py.File(annotation_file,'r')
            geneID = f['geneID'][:]
            gene_chrom = f['chrom'][:]
            gene_start = f['start'][:]
            gene_end = f['end'][:]
            gene_strand = f['strand'][:]
            gene_function = f['function'][:]
            f.close()
        except:
            print('Splitter:: format annotation file not valid')

        # if funct is not None, it has to be a list
        if funct is not None and funct!=list:   funct=[funct]

        windows = []
        nSnps   = []
        Igene   = []
        #2. calculates windows 
        for gene_i in range(geneID.shape[0]):
            if funct is not None:
                if gene_function[gene_i] not in funct:
                    Igene.append(False)
                    continue
            wnd = [gene_chrom[gene_i],gene_start[gene_i]-cis,gene_end[gene_i]+cis]
            Ir = (self.chrom==wnd[0])*(self.pos>=wnd[1])*(self.pos<=wnd[2])
            _nSnps = Ir.sum()
            if _nSnps>=minSnps and _nSnps<=maxSnps:
                windows.append(wnd)
                nSnps.append(_nSnps)
                Igene.append(True)
            else:
                Igene.append(False)
        Igene = SP.array(Igene)
        self.info['nSnps'] = SP.array(nSnps)
        self.info['geneID']        = geneID[Igene]
        self.info['gene_start']    = gene_start[Igene]
        self.info['gene_end']      = gene_end[Igene]
        self.info['gene_strand']   = gene_strand[Igene]
        self.info['gene_function'] = gene_function[Igene]
        return SP.array(windows)
353074
def parse_args():
    """"""
    parser = argparse.ArgumentParser(
        prog="geneparse-indexer",
        description="Genotype file indexer."
    )

    # IMPUTE2 files
    group = parser.add_argument_group("IMPUTE2 index")
    group.add_argument(
        "--impute2", metavar="IMPUTE2", type=str, nargs="+",
        help="Index an IMPUTE2 genotype file format. The file can be plain "
             "text or bgzipped.",
    )

    # BGEN files
    group = parser.add_argument_group("BGEN index")
    group.add_argument(
        "--bgen", metavar="BGEN", type=str, nargs="+",
        help="Index a BGEN genotype file. This requires 'bgenix' to be in the "
             "PATH.",
    )
    group.add_argument(
        "--legacy", action="store_true",
        help="Index the file using the '-with-rowid' option. This flag "
             "enables compatibility with SQLITE prior to version 3.8.2. See "
             "https://bitbucket.org/gavinband/bgen/wiki/bgenix for more "
             "information.",
    )

    return parser.parse_args()
800909
def process_mgi_relationship_transgene_genes(self, limit=None):
        
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        LOG.info("getting transgene genes")
        raw = '/'.join((self.rawdir, 'mgi_relationship_transgene_genes'))
        geno = Genotype(graph)
        col = [
            'rel_key', 'allele_key', 'allele_id', 'allele_label', 'category_key',
            'category_name', 'property_key', 'property_name', 'gene_num'
        ]
        with open(raw, 'r', encoding="utf8") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            header = next(filereader)
            if header != col:
                LOG.error('expected columns:  %s\n\tBut got:\n%s', col, header)
            for row in filereader:
                # rel_key,
                allele_key = int(row[col.index('allele_key')])
                allele_id = row[col.index('allele_id')]
                # allele_label,
                # category_key,
                # category_name,
                # property_key,
                # property_name,
                gene_num = int(row[col.index('gene_num')])

                if self.test_mode and allele_key not in self.test_keys.get('allele')\
                        and gene_num not in self.test_ids:
                    continue

                gene_id = 'NCBIGene:' + str(gene_num)

                # geno.addParts(gene_id, allele_id, self.globaltt['has_variant_part'])
                seqalt_id = self.idhash['seqalt'].get(allele_key)
                if seqalt_id is None:
                    seqalt_id = allele_id
                geno.addSequenceDerivesFrom(seqalt_id, gene_id)

                if not self.test_mode and limit is not None and \
                        filereader.line_num > limit:
                    break

        return
153005
public int geneCount() {
		int count = 0;
		for (int i = 0, n = _chromosomes.length(); i < n; ++i) {
			count += _chromosomes.get(i).length();
		}
		return count;
	}
281734
public void doMutation(PermutationSolution<T> solution) {
    int permutationLength ;
    permutationLength = solution.getNumberOfVariables() ;

    if ((permutationLength != 0) && (permutationLength != 1)) {
      if (mutationRandomGenerator.getRandomValue() < mutationProbability) {
        int pos1 = positionRandomGenerator.getRandomValue(0, permutationLength - 1);
        int pos2 = positionRandomGenerator.getRandomValue(0, permutationLength - 1);

        while (pos1 == pos2) {
          if (pos1 == (permutationLength - 1))
            pos2 = positionRandomGenerator.getRandomValue(0, permutationLength - 2);
          else
            pos2 = positionRandomGenerator.getRandomValue(pos1, permutationLength - 1);
        }

        T temp = solution.getVariableValue(pos1);
        solution.setVariableValue(pos1, solution.getVariableValue(pos2));
        solution.setVariableValue(pos2, temp);
      }
    }
  }
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
1230699
def g_to_n(self, g_interval):
        """"""

        # This code is extremely convoluted. To begin with, it
        # confuses interbase intervals for a *single* hgvs position
        # with intervals of *two* hgvs positions.  This is the origin
        # of bug #205 (now closed). Although the bug is itself fixed,
        # the mappers really need an overhaul. See #208.

        def _hgvs_offset(g_position, grs, gre, strand):
            """ Calculates the HGVS coordinate offset from a given genomic position"""
            if g_position == grs or g_position == gre:
                return 0
            mid = (grs + gre) / 2
            if g_position < mid or (g_position == mid and strand == 1):
                offset = g_position - grs
            else:
                offset = g_position - gre
            return offset

        def map_g_to_n_pos(pos):
            g_ci = _hgvs_coord_to_ci(pos, pos)
            # frs, fre = (f)orward (r)na (s)tart & (e)nd; forward w.r.t. genome
            frs, fre = self.im.map_ref_to_tgt(
                g_ci[0] - self.gc_offset, g_ci[1] - self.gc_offset, max_extent=False)
            if self.strand == -1:
                frs, fre = self.tgt_len - fre, self.tgt_len - frs
            # get BaseOffsetPosition information for n.
            if self.strand == 1:
                grs, gre = self.im.map_tgt_to_ref(frs, fre, max_extent=False)
            elif self.strand == -1:
                grs, gre = self.im.map_tgt_to_ref((self.tgt_len - fre), (self.tgt_len - frs),
                                                  max_extent=False)
            grs, gre = grs + self.gc_offset, gre + self.gc_offset
            # 0-width interval indicates an intron.  Need to calculate offsets but we"re are in ci coordinates
            # requires adding 1 strategically to get the HGVS position
            # (shift coordinates to 3' end of the ref nucleotide)

            if frs == fre:
                start_offset = _hgvs_offset(g_ci[0] + 1, grs, gre + 1, self.strand)
                end_offset = _hgvs_offset(g_ci[1], grs, gre + 1, self.strand)
            else:
                start_offset = _hgvs_offset(g_ci[0], grs, gre, self.strand)
                end_offset = _hgvs_offset(g_ci[1], grs, gre, self.strand)
            if self.strand == -1:
                start_offset, end_offset = self.strand * end_offset, self.strand * start_offset
            if start_offset > 0:
                frs -= 1
            if end_offset < 0:
                fre += 1
            start_base = _ci_to_hgvs_coord(frs, fre)[0]
            end_base = _ci_to_hgvs_coord(frs, fre)[1]
            return start_base, start_offset, end_base, end_offset

        start_bo = map_g_to_n_pos(g_interval.start.base)[0:2]
        end_bo = start_bo if g_interval.start.base == g_interval.end.base else map_g_to_n_pos(
            g_interval.end.base)[2:4]
        if self.strand == -1:
            start_bo, end_bo = end_bo, start_bo

        return hgvs.location.BaseOffsetInterval(
            start=hgvs.location.BaseOffsetPosition(
                base=start_bo[0], offset=start_bo[1], datum=Datum.SEQ_START),
            end=hgvs.location.BaseOffsetPosition(
                base=end_bo[0], offset=end_bo[1], datum=Datum.SEQ_START),
            uncertain=g_interval.uncertain)
1150509
def parse_dereplicated_uc(dereplicated_uc_lines):
    

    dereplicated_clusters = {}

    seed_hit_ix = 0
    seq_id_ix = 8
    seed_id_ix = 9

    for line in dereplicated_uc_lines:
        if line.startswith("#") or len(line.strip()) == 0:
            continue
        curr_line = line.strip().split('\t')
        if curr_line[seed_hit_ix] == "S":
            dereplicated_clusters[curr_line[seq_id_ix]] = []
        if curr_line[seed_hit_ix] == "H":
            curr_seq_id = curr_line[seq_id_ix]
            dereplicated_clusters[curr_line[seed_id_ix]].append(curr_seq_id)

    return dereplicated_clusters