<commit_before>#!/usr/bin/env python
# -*- coding: UTF-8 -*-

"""
Scripts for the pineapple genome manuscript (unpublished).
"""

import sys

from jcvi.graphics.base import plt, savefig, panel_labels
from jcvi.graphics.karyotype import Karyotype
from jcvi.graphics.synteny import Synteny, draw_gene_legend
from jcvi.graphics.glyph import TextCircle
from jcvi.apps.base import OptionParser, ActionDispatcher


def main():

    actions = (
        ('ploidy', 'plot pineapple macro-synteny (requires data)'),
            )
    p = ActionDispatcher(actions)
    p.dispatch(globals())


def ploidy(args):
    """
    %prog cotton seqids karyotype.layout mcscan.out all.bed synteny.layout

    Build a figure that calls graphics.karyotype to illustrate the high ploidy
    of WGD history of pineapple genome. The script calls both graphics.karyotype
    and graphic.synteny.
    """
    p = OptionParser(ploidy.__doc__)
    p.add_option("--switch", help="Rename the seqid with two-column file")
    opts, args, iopts = p.set_image_options(args, figsize="9x7")

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

    seqidsfile, klayout, datafile, bedfile, slayout = args

    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    Karyotype(fig, root, seqidsfile, klayout)
    Synteny(fig, root, datafile, bedfile, slayout, switch=opts.switch)

    # legend showing the orientation of the genes
    draw_gene_legend(root, .27, .37, .52)

    # annotate the WGD events
    fc = 'lightslategrey'
    x = .09
    radius = .012
    TextCircle(root, x, .825, r'$\tau$', radius=radius, fc=fc)
    TextCircle(root, x, .8, r'$\sigma$', radius=radius, fc=fc)
    TextCircle(root, x, .72, r'$\rho$', radius=radius, fc=fc)
    for ypos in (.825, .8, .72):
        root.text(.12, ypos, r"$\times2$", color=fc, ha="center", va="center")
    root.plot([x, x], [.85, .775], ":", color=fc, lw=2)
    root.plot([x, x], [.75, .675], ":", color=fc, lw=2)

    labels = ((.04, .96, 'A'), (.04, .54, 'B'))
    panel_labels(root, labels)

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    pf = "pineapple-karyotype"
    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)


if __name__ == '__main__':
    main()<commit_msg>add scripts related to pseudomolecule in projects.pineapple<commit_after>#!/usr/bin/env python
# -*- coding: UTF-8 -*-

"""
Scripts for the pineapple genome manuscript (unpublished).
"""

import sys
import logging

from jcvi.formats.sizes import Sizes
from jcvi.graphics.base import plt, savefig, panel_labels
from jcvi.graphics.karyotype import Karyotype
from jcvi.graphics.synteny import Synteny, draw_gene_legend
from jcvi.graphics.glyph import TextCircle
from jcvi.apps.base import OptionParser, ActionDispatcher


def main():

    actions = (
        ('ploidy', 'plot pineapple macro-synteny (requires data)'),
        # build pseudomolecule
        ('agp', 'make agp file'),
        ('breakpoints', 'make breakpoints'),
        ('check', 'check agreement'),
            )
    p = ActionDispatcher(actions)
    p.dispatch(globals())


def ploidy(args):
    """
    %prog cotton seqids karyotype.layout mcscan.out all.bed synteny.layout

    Build a figure that calls graphics.karyotype to illustrate the high ploidy
    of WGD history of pineapple genome. The script calls both graphics.karyotype
    and graphic.synteny.
    """
    p = OptionParser(ploidy.__doc__)
    p.add_option("--switch", help="Rename the seqid with two-column file")
    opts, args, iopts = p.set_image_options(args, figsize="9x7")

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

    seqidsfile, klayout, datafile, bedfile, slayout = args

    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    Karyotype(fig, root, seqidsfile, klayout)
    Synteny(fig, root, datafile, bedfile, slayout, switch=opts.switch)

    # legend showing the orientation of the genes
    draw_gene_legend(root, .27, .37, .52)

    # annotate the WGD events
    fc = 'lightslategrey'
    x = .09
    radius = .012
    TextCircle(root, x, .825, r'$\tau$', radius=radius, fc=fc)
    TextCircle(root, x, .8, r'$\sigma$', radius=radius, fc=fc)
    TextCircle(root, x, .72, r'$\rho$', radius=radius, fc=fc)
    for ypos in (.825, .8, .72):
        root.text(.12, ypos, r"$\times2$", color=fc, ha="center", va="center")
    root.plot([x, x], [.85, .775], ":", color=fc, lw=2)
    root.plot([x, x], [.75, .675], ":", color=fc, lw=2)

    labels = ((.04, .96, 'A'), (.04, .54, 'B'))
    panel_labels(root, labels)

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    pf = "pineapple-karyotype"
    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)


scaffold = "scaffold_"


def check(args):
    fp = open("assembly-order.txt")
    fp.next()
    d = {}
    for row in fp:
        atoms = row.split()
        scaf, tag, linkage, no = atoms[:4]
        d[scaf] = tag

    fp = open("chimeric-scaffolds.txt")
    fp.next()
    for row in fp:
        old, new, tag, start, end = row.strip().split("\t")
        if new not in d:
            print new, "not in sheet1"
            continue
        if d[new] != tag:
            print "{0} => {1} in sheet1 but {2} in sheet2".format(new, d[new], tag)


def agp(args):
    fp = open("assembly-order.txt")
    fp.next()
    sizes = Sizes("SCAFFOLD-SPLIT.fasta").mapping
    for row in fp:
        atoms = row.split()
        assert len(atoms) in (4, 5)
        if len(atoms) == 4:
            atoms.append('?')
        scaf, tag, linkage, no, strand = atoms
        strand = strand.lower()
        strand = {'f': '+', 'r': '-', '?': '?'}[strand]
        scaf = "scaffold_" + scaf
        scaf_size = sizes[scaf]
        linkage = "LG{0:02d}".format(ord(linkage.lower()) - ord('a') + 1)
        print "\t".join(str(x) for x in \
                    (scaf, 0, scaf_size, linkage, 1000, strand))


def breakpoints(args):
    fp = open("chimeric-scaffolds.txt")
    fp.next()
    scaffolds = set()
    nbreaks = 0
    for row in fp:
        atoms = row.strip().split("\t")
        if len(atoms) == 3:
            continue
        old, new, tag, start, end = atoms
        old = scaffold + old
        start, end = int(start), int(end)
        if start >= end:
            logging.warn("{0} {1} >= {2}".format(old, start, end))
            start, end = end, start
        print "\t".join(str(x) for x in (old, start - 1, end))
        nbreaks += 1
        scaffolds.add(old)
    print >> sys.stderr, "{0} breakpoints in total, {1} scaffolds broken"\
                    .format(nbreaks, len(scaffolds))


if __name__ == '__main__':
    main()
{'fraction_non_alphanumeric': 0.15421185559352252, 'fraction_numerical': 0.021690684890803745, 'mean_word_length': 2.75041782729805, 'pattern_counts': {'":': 4, '<': 3, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2186761', 'n_tokens_mistral': 2564, 'n_tokens_neox': 2398, 'n_words': 633}