Summary of the chromosomes of SNPs.

    Parameters
    ----------
    snps : pandas.DataFrame

    Returns
    -------
    str
        human-readable listing of chromosomes (e.g., '1-3, MT'), empty str if no chromosomes
489997
def get_chromosomes_summary(snps):
    

    if isinstance(snps, pd.DataFrame):
        chroms = list(pd.unique(snps["chrom"]))

        int_chroms = [int(chrom) for chrom in chroms if chrom.isdigit()]
        str_chroms = [chrom for chrom in chroms if not chrom.isdigit()]

        # https://codereview.stackexchange.com/a/5202
        def as_range(iterable):
            l = list(iterable)
            if len(l) > 1:
                return "{0}-{1}".format(l[0], l[-1])
            else:
                return "{0}".format(l[0])

        # create str representations
        int_chroms = ", ".join(
            as_range(g)
            for _, g in groupby(int_chroms, key=lambda n, c=count(): n - next(c))
        )
        str_chroms = ", ".join(str_chroms)

        if int_chroms != "" and str_chroms != "":
            int_chroms += ", "

        return int_chroms + str_chroms
    else:
        return ""
1128962
def summary(args):
    
    from jcvi.utils.natsort import natsort_key

    p = OptionParser(summary.__doc__)
    p.add_option("--suffix", default="Mb",
            help="make the base pair counts human readable [default: %default]")
    p.add_option("--ids",
            help="write the ids that have >= 50% N's [default: %default]")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    idsfile = opts.ids
    header = "Seqid Real N's Total %_real".split()
    if idsfile:
        idsfile = open(idsfile, "w")
        nids = 0

    data = []
    for fastafile in args:
        for rec in SeqIO.parse(fastafile, "fasta"):
            seqlen = len(rec)
            nns = rec.seq.count('n') + rec.seq.count('N')
            reals = seqlen - nns
            pct = reals * 100. / seqlen
            pctreal = "{0:.1f}%".format(pct)
            if idsfile and pct < 50:
                nids += 1
                print(rec.id, file=idsfile)

            data.append((rec.id, reals, nns, seqlen, pctreal))

    data.sort(key=natsort_key)
    ids, reals, nns, seqlen, pctreal = zip(*data)
    reals = sum(reals)
    nns = sum(nns)
    seqlen = sum(seqlen)
    pctreal = "{0:.1f}%".format(reals * 100. / seqlen)
    data.append(("Total", reals, nns, seqlen, pctreal))

    write_csv(header, data, sep=" ", filename=opts.outfile, thousands=True)
    if idsfile:
        logging.debug("A total of {0} ids >= 50% N's written to {1}.".\
                      format(nids, idsfile.name))
        idsfile.close()

    return reals, nns, seqlen
683988
def gffatributes2ids(s):
    
    Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
    if '=' in s:
        d=dict([i.split('=') for i in s.split(';')])
        if 'Parent' in d:
            d[d['Parent'].split(':')[0]+'_id']=d['Parent'].split(':')[1]
        Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
        if 'Name' in d:    
            Name=d['Name']
        if 'gene_id' in d:    
            gene_id=d['gene_id']
        if 'transcript_id' in d:    
            transcript_id=d['transcript_id']
        if 'protein_id' in d:    
            protein_id=d['protein_id']
        if 'exon_id' in d:    
            exon_id=d['exon_id']
    return Name,gene_id,transcript_id,protein_id,exon_id
747531
def discrepant_snps(self):
        
        df = self._discrepant_positions.append(self._discrepant_genotypes)
        if len(df) > 1:
            df = df.drop_duplicates()
        return df
1114816
def parse_snpeff_log(self, f):
        

        keys = {
            '# Summary table': [
                'Genome', 'Number_of_variants_before_filter',
                'Number_of_known_variants', 'Number_of_effects',
                'Genome_total_length', 'Change_rate'
            ],
            '# Effects by impact': [ 'HIGH', 'LOW', 'MODERATE', 'MODIFIER' ],
            '# Effects by functional class': [ 'MISSENSE', 'NONSENSE', 'SILENT', 'Missense_Silent_ratio' ],
            '# Hom/Het table': ['Het', 'Hom', 'Missing'],
            '# Ts/Tv summary': [ 'Transitions', 'Transversions', 'Ts_Tv_ratio' ],
            '# Count by effects': 'all',
            '# Count by genomic region': 'all'
        }
        parsed_data = {}
        section = None
        for l in f['f']:
            l = l.strip()
            if l[:1] == '#':
                section = l
                self.snpeff_section_totals[section] = dict()
                continue
            s = l.split(',')

            # Quality values / counts
            if section == '# Quality':
                quals = OrderedDict()
                if l.startswith('Values'):
                    values = [int(c) for c in l.split(',')[1:] ]
                    counts = f['f'].readline()
                    counts = [int(c) for c in counts.split(',')[1:] ]
                    c = 0
                    total = sum(counts)
                    for i, v in enumerate(values):
                        if c < (total * 0.995):
                            quals[v] = counts[i]
                            c += counts[i]
                if len(quals) > 0:
                    self.snpeff_qualities[f['s_name']] = quals

            # Everything else
            elif section in keys:
                if keys[section] == 'all' or any([k in s[0].strip() for k in keys[section]]):
                    try:
                        parsed_data[ s[0].strip() ] = float(s[1].strip())
                    except ValueError:
                        parsed_data[ s[0].strip() ] = s[1].strip()
                    except IndexError:
                        pass
                    else:
                        # Parsing the number worked - add to totals
                        try:
                            self.snpeff_section_totals[section][s[0].strip()] += parsed_data[ s[0].strip() ]
                        except KeyError:
                            self.snpeff_section_totals[section][s[0].strip()] = parsed_data[ s[0].strip() ]
                    if len(s) > 2 and s[2][-1:] == '%':
                        parsed_data[ '{}_percent'.format(s[0].strip()) ] = float(s[2][:-1])

        if len(parsed_data) > 0:
            if f['s_name'] in self.snpeff_data:
                log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
            self.add_data_source(f)
            self.snpeff_data[f['s_name']] = parsed_data
737601
def toStr(self):
        
        # Note that the string representation contains much less
        # information than the HTML summary. E.g., it does not contain the
        # unique (de-duplicated, by id) read count, since that is only computed
        # when we are making combined FASTA files of reads matching a
        # pathogen.
        readCountGetter = itemgetter('readCount')
        result = []
        append = result.append

        append(self._title())
        append('')

        for pathogenName in sorted(self.pathogenNames):
            samples = self.pathogenNames[pathogenName]
            sampleCount = len(samples)
            append('%s (in %d sample%s)' %
                   (pathogenName,
                    sampleCount, '' if sampleCount == 1 else 's'))
            for sampleName in sorted(samples):
                proteins = samples[sampleName]['proteins']
                proteinCount = len(proteins)
                totalReads = sum(readCountGetter(p) for p in proteins.values())
                append('  %s (%d protein%s, %d read%s)' %
                       (sampleName,
                        proteinCount, '' if proteinCount == 1 else 's',
                        totalReads, '' if totalReads == 1 else 's'))
                for proteinName in sorted(proteins):
                    append(
                        '    %(coverage).2f\t%(medianScore).2f\t'
                        '%(bestScore).2f\t%(readCount)4d\t%(hspCount)4d\t'
                        '%(index)3d\t%(proteinName)s'
                        % proteins[proteinName])
            append('')

        return '\n'.join(result)
373018
def load_from_args_as_dataframe(args):
    '''
    
    '''
    if not args.variants and not args.single_variant:
        return None

    if args.variant_source_name:
        variant_source_names = util.expand(
            args.variant_source_name,
            'variant_source_name',
            'variant source',
            len(args.variants))
    else:
        variant_source_names = util.drop_prefix(args.variants)

    variant_to_sources = collections.defaultdict(list)

    dfs = []
    for i in range(len(args.variants)):
        name = variant_source_names[i]
        prefix = (
            'metadata:' if len(args.variants) == 1 else "metadata:%s:" % name)
        df = load_as_dataframe(
            args.variants[i],
            name=name,
            genome=args.genome,
            max_variants=args.max_variants_per_source,
            only_passing=not args.include_failing_variants,
            metadata_column_prefix=prefix)

        if df.shape[0] == 0:
            logging.warn("No variants loaded from: %s" % args.variants[i])
        else:
            for variant in df.variant:
                variant_to_sources[variant].append(name)
            dfs.append(df)

    if args.single_variant:
        variants = []
        extra_args = {}
        if args.genome:
            extra_args = {
                'ensembl': varcode.reference.infer_genome(args.genome)
            }
        for (locus_str, ref, alt) in args.single_variant:
            locus = Locus.parse(locus_str)
            variant = varcode.Variant(
                    locus.contig,
                    locus.inclusive_start,
                    ref,
                    alt,
                    **extra_args)
            variants.append(variant)
            variant_to_sources[variant].append("commandline")
        dfs.append(variants_to_dataframe(variants))

    df = dfs.pop(0)
    for other_df in dfs:
        df = pandas.merge(
            df,
            other_df,
            how='outer',
            on=["variant"] + STANDARD_DATAFRAME_COLUMNS)

    genomes = df["genome"].unique()
    if len(genomes) > 1:
        raise ValueError(
                "Mixing references is not supported. "
                "Reference genomes: %s" % (", ".join(genomes)))

    df["sources"] = [" ".join(variant_to_sources[v]) for v in df.variant]

    # Apply filters:
    if args.ref:
        df = df.ix[df.ref.isin(args.ref)]
    if args.alt:
        df = df.ix[df.alt.isin(args.alt)]
    loci = loci_util.load_from_args(
        util.remove_prefix_from_parsed_args(args, "variant"))
    if loci is not None:
        df = df.ix[[
            loci.intersects(pileup_collection.to_locus(v))
            for v in df.variant
        ]]
    return df
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)
194455
def df_representative_structures(self):
        """"""
        rep_struct_pre_df = []
        df_cols = ['gene', 'id', 'is_experimental', 'file_type', 'structure_file']

        for g in self.genes_with_a_representative_structure:
            repdict = g.protein.representative_structure.get_dict(df_format=True, only_attributes=df_cols)
            repdict['gene'] = g.id
            rep_struct_pre_df.append(repdict)

        df = pd.DataFrame.from_records(rep_struct_pre_df, columns=df_cols).set_index('gene')
        if df.empty:
            log.warning('Empty dataframe')
            return df
        else:
            return ssbio.utils.clean_df(df)
1238640
def dtypes_summary(df):
    
    output_df = pd.DataFrame([])
    row_count = df.shape[0]
    row_indexes = ['rows_numerical','rows_string','rows_date_time','category_count','largest_category','rows_na','rows_total']
    for colname in df:
        data = df[colname] # data is the pandas series associated with this column
        # number of numerical values in the column
        rows_numerical = pd.to_numeric(data,errors = 'coerce').count()
        # number of values that can't be coerced to a numerical
        rows_string = row_count - rows_numerical
        # number of values that can be coerced to a date-time object
        rows_date_time = pd.to_datetime(data,errors = 'coerce',infer_datetime_format = True).count()
        # categories in column
        value_counts = data.value_counts().reset_index()
        # number of different values in the dataframe
        categories = len(value_counts)
        # largest category
        largest_category = value_counts.iloc[0,1]
        # number of null/missing values
        rows_na = data.isnull().sum()
        # build the output list
        output_data = [rows_numerical, rows_string, rows_date_time, categories, 
                       largest_category,rows_na,row_count]
        # add to dataframe
        output_df.loc[:,colname] = pd.Series(output_data)

    # row names
    output_df.index = row_indexes
    return output_df
996778
def print_stats(genomes):
    
    header = ['#genome', 'contig', 'locus', 'position', 'strand', 'length', \
                'sample', 'coverage', \
                'ref. Pn/Ps', 'ref SNP density', \
                'consensus Pn/Ps', 'consensus SNP density', \
                'product']
    print('\t'.join(header))
    for genome, contigs in list(genomes.items()):
        for contig, samples in list(contigs.items()):
            for sample, stats in list(samples.items()):
                for locus, rates in list(stats['sub_rates'].items()):
                    length = rates['info']['length']
                    position, strand = rates['info']['position']
                    position = '%s-%s' % position
                    out = [genome, contig, locus, position, strand, length, \
                            sample, '%.2f' % (rates['cov']), \
                            rates['ref PnPs'], rates['ref SNP density'], \
                            rates['consensus PnPs'], rates['consensus SNP density'], \
                            rates['info']['product'][0]]
                    print('\t'.join([str(i) for i in out]))
1197641
def calc_frequencies(genomes, bp_table, min_cov, min_per):
    
    nucs = ['A', 'T', 'G', 'C', 'N']
    if bp_table is not False:
        bp_table = open(bp_table, 'w')
        header = ['#genome', 'contig', 'sample', 'position', \
                    'reference', 'ref. frequency', \
                    'consensus', 'con. frequency', \
                    'A', 'T', 'G', 'C', 'N', '# insertions', '# deletions']
        print('\t'.join(header), file=bp_table)
    for genome, contigs in list(genomes.items()):
        for contig, samples in list(contigs.items()):
            for sample, stats in list(samples.items()):
                for pos, ps in enumerate(stats['bp_stats'], 1):
                    coverage = sum([ps[nuc] for nuc in nucs])
                    for nuc in nucs:
                        # make sure support for base passes thresholds
                        nuc_cov = ps[nuc]
                        if coverage == 0:
                            nuc_per = 0
                        else:
                            nuc_per = (float(nuc_cov)/coverage)*100
                        if nuc_cov < min_cov or nuc_per < min_per:
                            del ps[nuc]
                    ps = find_consensus(ps)
                    genomes[genome][contig][sample][pos] = ps
                    if bp_table is not False:
                        out = [genome, contig, sample, pos]
                        for i in ['ref', 'reference frequency', \
                                  'consensus', 'consensus frequency', \
                                  'A', 'T', 'G', 'C', 'N', \
                                  'In', 'Del']:
                            try:
                                if i == 'consensus':
                                    out.append(ps[i][0])
                                elif i in ['In', 'Del']:
                                    out.append(len(ps[i]))
                                else:
                                    out.append(ps[i])
                            except:
                                out.append('n/a')
                        print('\t'.join([str(i) for i in out]), file=bp_table)
    return genomes
1331368
def get_sample_info(self, individual=-1):
        
        if isinstance(individual, str):
            individual = self.individuals[individual]
        extra = OrderedDict()
        for format_ in self.format:
            index = getattr(self, format_)
            if index != -1:
                if format_ == 'GT':
                    d = self.genotype
                elif format_ == 'GQ':
                    d = self.genome_quality
                elif format_ == 'DP':
                    d = self.depth
                if individual == -1:
                    if len(d) != len(self.samples):
                        [self.parse_sample(i) for i in six.moves.range(len(self.samples))]
                    extra[format_] = [d[i] for i in six.moves.range(len(d))]
                else:
                    if individual not in d:
                        self.parse_sample(individual)
                    extra[format_] = [d[individual]]
        return extra
125788
def _freqs_by_chromosome(in_file, params, somatic_info):
    
    freqs = []
    coords = []
    cur_chrom = None
    with pysam.VariantFile(in_file) as bcf_in:
        for rec in bcf_in:
            if _is_biallelic_snp(rec) and _passes_plus_germline(rec) and chromhacks.is_autosomal(rec.chrom):
                if cur_chrom is None or rec.chrom != cur_chrom:
                    if cur_chrom and len(freqs) > 0:
                        yield cur_chrom, freqs, coords
                    cur_chrom = rec.chrom
                    freqs = []
                    coords = []
                stats = _tumor_normal_stats(rec, somatic_info)
                if tz.get_in(["tumor", "depth"], stats, 0) > params["min_depth"]:
                    # not a ref only call
                    if len(rec.samples) == 0 or sum(rec.samples[somatic_info.tumor_name].allele_indices) > 0:
                        freqs.append(tz.get_in(["tumor", "freq"], stats))
                        coords.append(rec.start)
        if cur_chrom and len(freqs) > 0:
            yield cur_chrom, freqs, coords
387635
def get_chromosome_priority(chrom, chrom_dict={}):
    
    priority = 0
    
    chrom = str(chrom).lstrip('chr')
    
    if chrom_dict:
        priority = chrom_dict.get(chrom, 0)
    
    else:
        try:
            if int(chrom) < 23:
                priority = int(chrom)
        except ValueError:
            if chrom == 'X':
                priority = 23
            elif chrom == 'Y':
                priority = 24
            elif chrom == 'MT':
                priority = 25
            else:
                priority = 26
    
    return str(priority)
1034692
def to_dataframe(self):
        """"""
        def row_from_variant(variant):
            return OrderedDict([
                ("chr", variant.contig),
                ("start", variant.original_start),
                ("ref", variant.original_ref),
                ("alt", variant.original_alt),
                ("gene_name", ";".join(variant.gene_names)),
                ("gene_id", ";".join(variant.gene_ids))
            ])
        rows = [row_from_variant(v) for v in self]
        if len(rows) == 0:
            # TODO: return a DataFrame with the appropriate columns
            return pd.DataFrame()
        return pd.DataFrame.from_records(rows, columns=rows[0].keys())
492865
def exons(self, contig=None, strand=None):
        
        # DataFrame with single column called "exon_id"
        exon_ids = self.exon_ids(contig=contig, strand=strand)
        return [
            self.exon_by_id(exon_id)
            for exon_id in exon_ids
        ]