%prog statstable *.gff

    Print gene statistics table.
547969
def statstable(args):
    
    p = OptionParser(statstable.__doc__)
    opts, args = p.parse_args(args)

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

    gff_files = args
    for metric in metrics:
        logging.debug("Parsing files in `{0}`..".format(metric))

        table = {}
        for x in gff_files:
            pf = op.basename(x).split(".")[0]
            numberfile = op.join(metric, pf + ".txt")
            ar = [int(x.strip()) for x in open(numberfile)]
            sum = SummaryStats(ar).todict().items()
            keys, vals = zip(*sum)
            keys = [(pf, x) for x in keys]
            table.update(dict(zip(keys, vals)))

        print(tabulate(table), file=sys.stderr)
581824
def yeasttruth(args):
    
    p = OptionParser(yeasttruth.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    pillars = args[0]
    gffiles = args[1:]
    aliases = {}
    pivot = {}
    for gffile in gffiles:
        is_pivot = op.basename(gffile).startswith("Saccharomyces_cerevisiae")
        gff = Gff(gffile)
        for g in gff:
            if g.type != "gene":
                continue
            for a in g.attributes["Alias"]:
                aliases[a] = g.accn
                if is_pivot:
                    pivot[a] = g.accn
    logging.debug("Aliases imported: {0}".format(len(aliases)))
    logging.debug("Pivot imported: {0}".format(len(pivot)))
    fw = open("yeast.aliases", "w")
    for k, v in sorted(aliases.items()):
        print("\t".join((k, v)), file=fw)
    fw.close()

    fp = open(pillars)
    pairs = set()
    fw = must_open(opts.outfile, "w")
    for row in fp:
        atoms = [x for x in row.split() if x != "---"]
        pps = [pivot[x] for x in atoms if x in pivot]
        atoms = [aliases[x] for x in atoms if x in aliases]
        for p in pps:
            for a in atoms:
                if p == a:
                    continue
                pairs.add(tuple(sorted((p, a))))

    for a, b in sorted(pairs):
        print("\t".join((a, b)), file=fw)
    fw.close()
325134
def run(_, data, out_dir):
    
    out = []
    vcinfo = get_active_vcinfo(data)
    if vcinfo:
        if dd.get_phenotype(data) == "normal" and "germline" in vcinfo:
            out.append(_bcftools_stats(data, out_dir, "germline", germline=True))
        elif dd.get_phenotype(data) != "germline":
            out.append(_bcftools_stats(data, out_dir))
            if "germline" in vcinfo:
                out.append(_bcftools_stats(data, out_dir, "germline", germline=True))
        else:
            out.append(_bcftools_stats(data, out_dir, germline=True))
    out.append(_snpeff_stats(data, out_dir))

    out = [item for item in out if item]
    if out:
        return {"base": out[0], "secondary": out[1:]}
483873
def make_pizzly_gtf(gtf_file, out_file, data):
    
    if file_exists(out_file):
        return out_file
    db = gtf.get_gtf_db(gtf_file)
    with file_transaction(data, out_file) as tx_out_file:
        with open(tx_out_file, "w") as out_handle:
            for gene in db.features_of_type("gene"):
                children = [x for x in db.children(id=gene)]
                for child in children:
                    if child.attributes.get("gene_biotype", None):
                        gene_biotype = child.attributes.get("gene_biotype")
                gene.attributes['gene_biotype'] = gene_biotype
                gene.source = gene_biotype[0]
                print(gene, file=out_handle)
                for child in children:
                    child.source = gene_biotype[0]
                    # gffread produces a version-less FASTA file
                    child.attributes.pop("transcript_version", None)
                    print(child, file=out_handle)
    return out_file
635109
def stats(args):
    
    from jcvi.utils.cbook import percentage

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

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

    blocksfile, = args
    fp = open(blocksfile)
    counts = defaultdict(int)
    total = orthologous = 0
    for row in fp:
        atoms = row.rstrip().split("\t")
        hits = [x for x in atoms[1:] if x != '.']
        counts[len(hits)] += 1
        total += 1
        if atoms[1] != '.':
            orthologous += 1

    print("Total lines: {0}".format(total), file=sys.stderr)
    for i, n in sorted(counts.items()):
        print("Count {0}: {1}".format(i, percentage(n, total)), file=sys.stderr)

    print(file=sys.stderr)

    matches = sum(n for i, n in counts.items() if i != 0)
    print("Total lines with matches: {0}".\
                format(percentage(matches, total)), file=sys.stderr)
    for i, n in sorted(counts.items()):
        if i == 0:
            continue

        print("Count {0}: {1}".format(i, percentage(n, matches)), file=sys.stderr)

    print(file=sys.stderr)
    print("Orthologous matches: {0}".\
                format(percentage(orthologous, matches)), file=sys.stderr)
730306
def print_table(lines, separate_head=True):
	""""""
	#Count the column width
	widths = []
	for line in lines:
			for i,size in enumerate([len(x) for x in line]):
					while i >= len(widths):
							widths.append(0)
					if size > widths[i]:
							widths[i] = size
	   
	#Generate the format string to pad the columns
	print_string = ""
	for i,width in enumerate(widths):
			print_string += "{" + str(i) + ":" + str(width) + "} | "
	if (len(print_string) == 0):
			return
	print_string = print_string[:-3]
	   
	#Print the actual data
	for i,line in enumerate(lines):
			print(print_string.format(*line))
			if (i == 0 and separate_head):
					print("-"*(sum(widths)+3*(len(widths)-1)))
289160
def print_table(*args, **kwargs):
    
    t = format_table(*args, **kwargs)
    click.echo(t)
1230495
def print_result(self, f=sys.stdout, verbose=False):
        
        var_count = len(self.betas)/2
        if verbose:
            results = [str(x) for x in [
                self.chr,
                self.pos,
                self.rsid,
                self.ph_label,
                self.non_miss,
                self.maj_allele,
                self.min_allele,
                "%.4e" % self.eff_alcount,
                "%.4e" % self.p_mvtest,
                "%.4e" % self.lmpv
            ]]

            for i in range(var_count):
                results.append("%.3e" % self.betas[i])
                results.append(self.stringify(self.beta_stderr[i]))
                results.append(self.stringify(self.beta_pvalues[i]))
                results.append(self.stringify(self.betas[i+var_count]))
                results.append(self.stringify(self.beta_stderr[i+var_count]))
                results.append(self.stringify(self.beta_pvalues[i+var_count]))
            print >> f, "\t".join([str(x) for x in results])

        else:
            print >> f, "\t".join(str(x) for x in [
                self.chr,
                self.pos,
                self.rsid,
                self.ph_label,
                self.non_miss,
                self.maj_allele,
                self.min_allele,
                "%.3e" % self.eff_alcount,
                "%.3e" % self.p_mvtest,
                "%.3e" % self.betas[1],
                self.stringify(self.beta_stderr[1]),
                "%.3e" % self.beta_pvalues[1],
                "%.3e" % self.betas[1+var_count],
                self.stringify(self.beta_stderr[1+var_count]),
                "%.3e" % self.beta_pvalues[1+var_count],
            ])
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
616859
def read_gtf_line(cols, field="name"):
    """"""
    field = field.lower()
    try:
        group = cols[2]
        attrs = cols[8].split(";")
        name = [attr.strip().split(" ")[1] for attr in attrs if attr.strip().split(" ")[0].lower().endswith(field)]
        if not name:
            name = [attr.strip().split(" ")[1] for attr in attrs if attr.strip().split(" ")[0].lower().endswith("gene_id")]
        if not name:
            name = ["None"]
        biotype = [attr.strip().split(" ")[1] for attr in attrs if attr.strip().split(" ")[0].lower().endswith("biotype")]
        if biotype:
            group = biotype[0]
        c = cols[0]
        s = int(cols[3])
        e = int(cols[4])
        st = cols[6]
        return [c, s, e, st, group, name[0]]
    except(Exception, e):
        logger.error(cols)
        logger.error("File is not in correct format")
        logger.error("Expect chr source feature start end . strand attributes")
        logger.error("Attributes are 'gene_name SNCA; gene_id ENSG; '")
        logger.error("The 3rd column is used as type of small RNA (like miRNA)")
        logger.error("at least should contains '; *name NAME; '")
        logger.error(e)
        raise
522731
def summary(args):
    
    from jcvi.formats.base import SetFile
    from jcvi.formats.bed import BedSummary
    from jcvi.utils.table import tabulate

    p = OptionParser(summary.__doc__)
    p.add_option("--isoform", default=False, action="store_true",
                 help="Find longest isoform of each id")
    p.add_option("--ids", help="Only include features from certain IDs")
    opts, args = p.parse_args(args)

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

    gff_file, = args
    ids = opts.ids

    if ids:
        ids = SetFile(ids)
        logging.debug("Total ids loaded: {0}".format(len(ids)))

        if opts.isoform:
            pids = set()
            gff = Gff(gff_file)
            for g in gff:
                if g.type != "mRNA":
                    continue
                if g.parent not in ids:
                    continue
                if "longest" not in g.attributes:
                    pids = set(x + ".1" for x in ids)
                    break
                if g.attributes["longest"][0] == "0":
                    continue
                pids.add(g.id)
            ids = pids
            logging.debug("After checking longest: {0}".format(len(ids)))

        # Collects aliases
        gff = Gff(gff_file)
        for g in gff:
            if g.name in ids:
                ids.add(g.id)
        logging.debug("Total ids including aliases: {0}".format(len(ids)))

    gff = Gff(gff_file)
    beds = defaultdict(list)
    for g in gff:
        if ids and not (g.id in ids or g.name in ids or g.parent in ids):
            continue

        beds[g.type].append(g.bedline)

    table = {}
    for type, bb in sorted(beds.items()):
        bs = BedSummary(bb)
        table[(type, "Features")] = bs.nfeats
        table[(type, "Unique bases")] = bs.unique_bases
        table[(type, "Total bases")] = bs.total_bases

    print(tabulate(table), file=sys.stdout)
669003
def print_genome_matrix(hits, fastas, id2desc, file_name):
    
    out = open(file_name, 'w')
    fastas = sorted(fastas)
    print('## percent identity between genomes', file=out)
    print('# - \t %s' % ('\t'.join(fastas)), file=out)
    for fasta in fastas:
        line = [fasta]
        for other in fastas:
            if other == fasta:
                average = '-'
            else:
                average = numpy.average([hits[fasta][other][i][3] for i in hits[fasta][other]])
            line.append(str(average))
        print('\t'.join(line), file=out)
    print('', file=out)
    print('## percent of orfs that are orthologous between genomes', file=out)
    print('# - \t %s' % ('\t'.join(fastas)), file=out)
    for fasta in fastas:
        line = [fasta]
        for other in fastas:
            if other == fasta:
                percent = '-'
            else:
                orthologs = float(len(hits[fasta][other]))
                orfs = float(len([i for i in id2desc if id2desc[i][0] == fasta]))
                percent = float(orthologs / orfs) * 100
            line.append(str(percent))
        print('\t'.join(line), file=out)
578194
def bcftools_stats_genstats_headers(self):
        
        stats_headers = OrderedDict()
        stats_headers['number_of_records'] = {
            'title': 'Vars',
            'description': 'Variations total',
            'min': 0, 'format': '{:,.0f}',
        }
        stats_headers['variations_hom'] = {
            'title': 'Hom',
            'description': 'Variations homozygous',
            'min': 0, 'format': '{:,.0f}',
        }
        stats_headers['variations_het'] = {
            'title': 'Het',
            'description': 'Variations heterozygous',
            'min': 0, 'format': '{:,.0f}',
        }
        stats_headers['number_of_SNPs'] = {
            'title': 'SNP',
            'description': 'Variation SNPs',
            'min': 0, 'format': '{:,.0f}',
        }
        stats_headers['number_of_indels'] = {
            'title': 'Indel',
            'description': 'Variation Insertions/Deletions',
            'min': 0, 'format': '{:,.0f}',
        }
        stats_headers['tstv'] = {
            'title': 'Ts/Tv',
            'description': 'Variant SNP transition / transversion ratio',
            'min': 0, 'format': '{:,.2f}',
        }
        stats_headers['number_of_MNPs'] = {
            'title': 'MNP',
            'description': 'Variation multinucleotide polymorphisms',
            'min': 0, 'format': '{:,.0f}', "hidden": True,
        }
        return stats_headers
1057932
def print_table(seqs, id2name, name):
    
    itable = open('%s.itable' % (name.rsplit('.', 1)[0]), 'w')
    print('\t'.join(['#sequence', 'gene', 'model', 'insertion', 'gene position', 'model position', 'length', 'orf?', 'intron?', 'orf?intron?', 'insertion', 'orf', 'intron']), file=itable)
    for seq, info in list(seqs.items()):
        gene, model, insertions = info
        name = id2name[seq]
        for i, ins in enumerate(insertions, 1):
            gene_pos, model_pos, length, iseq, \
                    orfs, introns, orfs_b, introns_b, orf_annotations = ins
            # append annotation to orf header
            for orf in orfs:
                parts = orf[0].split()
                annotation = orf_annotations[parts[0].split('>')[1]]
                orf[0] = '%s %s %s' % (parts[0], annotation, ' '.join(parts[1:]))
            # get orf position
            gene_pos = '-'.join([str(j) for j in gene_pos])
            # check if orf, intron is present
            if orfs_b is True or introns_b is True:
                orfs_introns_b = True
            else:
                orfs_introns_b = False
            out = [name, gene, model, i, gene_pos, model_pos, length, orfs_b, introns_b, orfs_introns_b]
            out.append('|'.join(iseq))
            out.append('|'.join(['|'.join(orf) for orf in orfs]))
            out.append('|'.join(['|'.join(intron) for intron in introns]))
            print('\t'.join([str(i) for i in out]), file=itable)
    itable.close()
168045
def verifybamid_general_stats_table(self):
		

		# create a dictionary to hold the columns to add to the general stats table
		headers = OrderedDict()
		# available columns are:
		#SEQ_ID RG  CHIP_ID #SNPS   #READS  AVG_DP  FREEMIX FREELK1 FREELK0 FREE_RH FREE_RA CHIPMIX CHIPLK1 CHIPLK0 CHIP_RH CHIP_RA DPREF   RDPHET  RDPALT
		#see https://genome.sph.umich.edu/wiki/VerifyBamID#Interpreting_output_files

		# add the CHIPMIX column if we have data
		if not self.hide_chip_columns:
			headers['CHIPMIX'] = dict(self.col_config_defaults, **{
				'title': 'Contamination (S+A)',
				'description': 'VerifyBamID: CHIPMIX -   Sequence+array estimate of contamination (NA if the external genotype is unavailable)'
			})

		# add the FREEMIX column. set the title and description
		headers['FREEMIX'] = dict(self.col_config_defaults, **{
			'title': 'Contamination (S)',
			'description': 'VerifyBamID: FREEMIX -   Sequence-only estimate of contamination.',
		})

		# pass the data dictionary and header dictionary to function to add to table.
		self.general_stats_addcols(self.verifybamid_data, headers)
231663
def bed(args):
    '''
    
    '''
    from jcvi.utils.cbook import gene_name

    p = OptionParser(bed.__doc__)
    p.add_option("--type", dest="type", default="gene",
            help="Feature type to extract, use comma for multiple [default: %default]")
    p.add_option("--key", default="ID", help="Key in the attributes to extract")
    p.add_option("--source",
            help="Source to extract from, use comma for multiple [default: %default]")
    p.add_option("--span", default=False, action="store_true",
            help="Use feature span in the score column")
    p.add_option("--score_attrib", dest="score_attrib", default=False,
            help="Attribute whose value is to be used as score in `bedline` [default: %default]")
    p.add_option("--append_source", default=False, action="store_true",
            help="Append GFF source name to extracted key value")
    p.add_option("--append_ftype", default=False, action="store_true",
            help="Append GFF feature type to extracted key value")
    p.add_option("--nosort", default=False, action="store_true",
            help="Do not sort the output bed file [default: %default]")
    p.set_stripnames(default=False)
    p.set_outfile()

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

    gffile, = args
    key = opts.key or None
    type = opts.type or set()
    source = opts.source or set()
    strip_names = opts.strip_names
    span = opts.span

    if opts.type:
        type = set(x.strip() for x in opts.type.split(","))
    if opts.source:
        source = set(x.strip() for x in opts.source.split(","))

    gff = Gff(gffile, key=key, append_source=opts.append_source, \
        append_ftype=opts.append_ftype, score_attrib=opts.score_attrib)
    b = Bed()

    for g in gff:
        if type and g.type not in type:
            continue
        if source and g.source not in source:
            continue

        bl = g.bedline
        if strip_names:
            bl.accn = gene_name(bl.accn)
        if span:
            bl.score = bl.span
        b.append(bl)

    sorted = not opts.nosort
    b.print_to_file(opts.outfile, sorted=sorted)
    logging.debug("Extracted {0} features (type={1} id={2})".\
                    format(len(b), ",".join(type), key))
110770
def stats(self, ops=(min, max, np.median, sum)):
        """"""
        names = [op.__name__ for op in ops]
        ops = [_zero_on_type_error(op) for op in ops]
        columns = [[op(column) for op in ops] for column in self.columns]
        table = type(self)().with_columns(zip(self.labels, columns))
        stats = table._unused_label('statistic')
        table[stats] = names
        table.move_to_start(stats)
        return table