# [-fmt {bed,gtf}] [-o OUTPUT_FILE]

# Imports
import circRNA as circ
from collections import defaultdict
import sys, os, re, argparse, natsort
import numpy as np
import csv
import pybedtools as pybed

# Utility functions
def eprint(*args, **kwargs):
    print(*args,  file=sys.stderr, **kwargs)
    dsp_control = circ.DisplayControl.instance()
    if dsp_control.verbose:
        print(*args,  file=sys.stderr, **kwargs)


def exonic_annotations(circ_rnas, exons):
    """
    Annotate the circrnas according to the exon annotations
    circ_rnas: an array of cirRNAs object
    exons: an array of circRNAs objects
    returns the circ_rnas array with the circRNAs annotated

    We annotate each circRNA's splice junction with the corresponding
    exons spice junctions using a dict of exons splice junctions
    This means that splice junction correspondances must be exact
    """
    eprint("\nCompute junction dictionary...\n")
    annotation_junctions = compute_junction_dict(exons)
    eprint("\nCircRNA annotation...\n")
    set_annotation(circ_rnas, annotation_junctions)

    return circ_rnas


def junction_key(chrom, pos):
    return "%s:%s" %(chrom, pos)


def compute_junction_dict(annot_array):
    """
    Compute a dictionnary of splic junctions:
    keys: chrom:pos
    values: list of tuples
    tuple: (exon, end_type)
    end_type: 3 (resp 5) for donor (resp acceptor)
    """
    junction_dict = defaultdict(list)
    seen_exons = defaultdict()
    for annot in annot_array:
        if annot.get_att("exon_id") in seen_exons:
            continue
        seen_exons[annot.get_att("exon_id")] = 1
        key_start = junction_key(annot.chrom, annot.start)
        key_end = junction_key(annot.chrom, annot.end)
        if annot.strand == "+":
            start_end_type = "5"
            end_end_type = "3"
        else:
            start_end_type = "3"
            end_end_type = "5"
        # annotating start
        junction_dict[key_start].append((annot, start_end_type))
        # annotating end
        junction_dict[key_end].append((annot, end_end_type))
    return junction_dict


def set_annotation(circ_rnas, annotation_junctions):
    """
    An annotation string would be very much appreciated
    """
    for circ_rna in circ_rnas:
        junction_start = junction_key(circ_rna.chrom, circ_rna.start)
        junction_end = junction_key(circ_rna.chrom, circ_rna.end)
        if junction_start in annotation_junctions:
            circ_rna.add_start_annotation(annotation_junctions[junction_start])
            circ_rna.add_start_transcript_id(annotation_junctions[junction_start])
            circ_rna.add_start_gene_id(annotation_junctions[junction_start])
        if junction_end in annotation_junctions:
            circ_rna.add_end_annotation(annotation_junctions[junction_end])
            circ_rna.add_end_transcript_id(annotation_junctions[junction_end])
            circ_rna.add_end_gene_id(annotation_junctions[junction_end])


def intronic_annotations(circ_rnas, exons):
    """
    Annotate the circrnas according to the introns annotations
    circ_rnas: an array of cirRNAs object
    exons: an array of circRNAs objects
    returns the circ_rnas array with the circRNAs intrno-annotated

    Because at least on circRNA junction is away from existing
    splice junctions, the annotation is solely based on overlap between
    the circRNA and an intron :
      - [start, end] interval must be within the intron
      - this interval must cover at least 90% of the intron
    """

    #eprint("\nCompute introns from exons...\n")
    transcripts = get_exons_per_transcript(exons)
    introns = circ.compute_intronic_positions(transcripts)

    #eprint("\nCompute intersection with circRNAs...\n")
    circ_rnas = annotate_intron_intersection(circ_rnas, introns)

    return circ_rnas


def annotate_intron_intersection(circ_rnas, introns):
    pybed_introns = annot_to_pybed(introns)
    pybed_circrnas = annot_to_pybed(circ_rnas)
    intersections = pybed_circrnas.intersect(pybed_introns, f=0.95, wo=True) 

    circrna_d = {circ.name: circ for circ in circ_rnas}
    introns_d = {intron.name: intron for intron in introns}
    for i in intersections:
        intron = introns_d[i[9]]
        circ = circrna_d[i[3]]
        overlap = i[-1]
        p_overlap = float(overlap)/intron.length # taille circ sur taille intron
        circ.annotate_intron(intron, p_overlap)

    return circ_rnas


def annot_to_pybed(annots):
    pybed_annot = []
    for a in annots:
        pybed_annot.append(a.topybed())
    return pybed.BedTool(pybed_annot).sort()


def infra_exonic_annotations(circ_rnas, exons):
    """
    """
    pybed_exons = annot_to_pybed_ife(exons)
    pybed_circrnas = annot_to_pybed_ife(circ_rnas)  
    intersections = pybed_circrnas.intersect(pybed_exons, f=0.95, wo=True)

    circrna_d = {circ.name: circ for circ in circ_rnas}
    exons_d = {exon.name: exon for exon in exons}
    for i in intersections:
        exon = exons_d[i[9]]
        circ = circrna_d[i[3]]
        circ.annotate_ife(exon)

    return circ_rnas


def annot_to_pybed_ife(annots):
    pybed_annot = []
    for a in annots:
        pybed_annot.append(a.topybed_ife())
    return pybed.BedTool(pybed_annot).sort()


def get_exons_per_transcript(annots):
    """
    returns a dictionnary
    key: ...
    """
    d_transcript = defaultdict(list)
    for exon in annots:
        transcript_id = exon.get_att("transcript_id")
        d_transcript[transcript_id].append(exon)
    return d_transcript


def write_annotated_circrnas(circ_annotated, outputfile):
    circ_rnas = natsort.natsorted(circ_annotated, key=lambda e: (e.chrom, e.start))
    header = "\t".join(["chrom", "start", "end", "strand", "nb_ccr", "circ_rna_name", "exons_id_start", "exons_id_end", 
                        "transcript_id_start", "transcript_id_end", "gene_id_start", "gene_id_end", 
                        "intron_name", "start_i", "end_i", "gene_id_i", "exon_id_i", 
                        'exon_id_ife', "exons_start_end_ife", "gene_id_ife"])
    with open(outputfile, "w") as fout:
        fout.write(header+"\n")
        for circ_rna in circ_rnas:
            intron_annot = circ_rna.intron_annotation 
            if len(intron_annot) == 0:
                start_i = ""
                end_i = ""
                name_i = ""
            else:
                for annot_i in intron_annot:
                    intron = annot_i[0]
                    start_i = intron.start
                    end_i = intron.end
                    name_i = intron.name
            att_d = dict(re.split(" |=", item) for item in circ_rna.attributes.split("; "))
            nb_ccr = att_d["nb_ccr"].replace('\"', '')
            s = [circ_rna.chrom, circ_rna.start, circ_rna.end, 
                circ_rna.strand, nb_ccr, circ_rna.name,
                circ_rna.get_start_annotation_str(), 
                circ_rna.get_end_annotation_str(),
                circ_rna.get_start_transcript_id_str(), 
                circ_rna.get_end_transcript_id_str(), 
                circ_rna.get_start_gene_id_str(), 
                circ_rna.get_end_gene_id_str(),
                name_i, start_i, end_i, 
                circ_rna.get_intron_annot_gene_str('gene_id'),
                circ_rna.get_intron_annot_str('exon_id'), 
                circ_rna.get_infra_exonic_annot_str("exon_id"),
                circ_rna.get_infra_exonic_exons_start_end(),
                circ_rna.get_infra_exonic_gene_annot_str("gene_id")]
            fout.write( "%s\n" % "\t".join(map(str, s)))


def main():

    # Exonic:
    eprint("\nIdentification of EXONIC circRNAs:\n")
    eprint("\nReading circRNA file...\n")
    circ_rnas = circ.read_annotation(args.circ_rna_file, fmt=args.annot_format)
    eprint("\nReading annotation file...\n")
    annots = circ.read_annotation(args.annotation_file, fmt=args.annot_format)

    eprint("\nExonic annotations..\n")
    circ_rnas = exonic_annotations(circ_rnas, annots)

    eprint("\nIntronic annotations..\n")
    circ_rnas = intronic_annotations(circ_rnas, annots)

    #####
    eprint("\nInfra-exonic annotations..\n")
    circ_rnas = infra_exonic_annotations(circ_rnas, annots)
    #####

    eprint("\nWrite annotations..\n")
    write_annotated_circrnas(circ_rnas, args.output)


def parse_arguments():
    parser = argparse.ArgumentParser(description='Annotation')
    parser.add_argument('-circ', '--circ_rna_file',
                        required=True, help='Circular RNA file path')
    parser.add_argument('-fmt', '--annot_format',
                        required=False, default="bed", choices=["bed","gtf"],
                        help='Annotation file format')
    parser.add_argument('-annot', '--annotation_file',
                        required=True, help='Annotataion file path')
    parser.add_argument('-oe', '--exonic_output_file',
                        required=False, default='exonic_circ_rnas_annot.out',
                        help='Exonic output file path')
    parser.add_argument('-oi', '--intronic_output_file',
                        required=False, default='annotation_introns.bed',
                        help='Intronic output file path')
    parser.add_argument('-o', '--output',
                        required=True,
                        help='Annotation output file path')
    parser.add_argument('--verbose', help='Print more info', action='store_true')

    args = parser.parse_args()

    dsp_control = circ.DisplayControl.instance()
    dsp_control.verbose = args.verbose

    return args


if __name__ == '__main__':
    args = parse_arguments()
    main()

                                    # [-min_cr CR_THRESHOLD] [-tol TOLERANCE]
                                    # [-fmt {bed,gtf}] [-o OUTPUT_FILE]

# Imports
import sys, os, re, argparse, natsort, time
if not sys.warnoptions:
    import warnings
    warnings.simplefilter("ignore")
from collections import defaultdict
import pandas as pd
import numpy as np
from networkx import (draw, DiGraph, Graph)
import networkx as netwx
import statistics
from statistics import mode
from time import sleep
from tqdm import tqdm
import circRNA as circ


# Utility functions
def eprint(*args, **kwargs):
    print(*args,  file=sys.stderr, **kwargs)
    dsp_control = circ.DisplayControl.instance()
    if dsp_control.verbose:
        print(*args,  file=sys.stderr, **kwargs)


def merge_junctions(df_a):
    """ Merge two dataframes
    :df_r1: dataframe of R1 reads
    :df_r2: dataframe of R2 reads
    :return the merged dataframe
    """
    df_merged = pd.concat(df_a,ignore_index=True)
    return df_merged


def read_junction_file(file, read_type=None):
    """ Read the data, adds the header of the table
        and adds a column with the corresponding read type.
    :file: file to read
    :read_type: read type (R1 or R2)
    :return the dataframe of reads
    """
    names = ["chr_donor", "pos_donor","strand_donor","chr_acceptor",
            "pos_acceptor","strand_acceptor","junction_type",
            "repeat_left","repeat_right","read_name","pos_s1",
            "CIGAR_s1","pos_s2","CIGAR_s2"] # Columns names of the input data
    # Read data as table with pandas:
    df = pd.read_table(file, sep = '\t', names=names, dtype=str,
                       converters = {'pos_acceptor':int, 'pos_donor':int})
    df['read_type'] = read_type # Add column "read" to specify the read type
    return df


def get_valid_circjunctions(df):
    """ Selects valid circular chimeric reads (ccr)
        according to the conditions of the "valid_ccr" function.
    :df: datatframe to filter
    return: the filtered dataframe
    """
    # Parsing dataframe to select valid rows
    # Boolean array with true if the read is a valid ccr:
    valid_ccrs = []
    #for index, row in tqdm(df.iterrows(), total=df.shape[0]):
    for index, row in circ.tqdm_wrapper(df):
        valid_ccrs.append(valid_ccr(row))
    df_circ_cr = df[valid_ccrs]
    ccr_a = [] # Array of circular junctions
    eprint(" "+"\n"+"Generating array of selected ccr..."+"\n")
    #for index, row in tqdm(df_circ_cr.iterrows(), total=df_circ_cr.shape[0]):
    for index, row in circ.tqdm_wrapper(df_circ_cr):
        ccr =  circ.CCR(row) # Create CCR object
        ccr_a.append(ccr) # Add chimeric reads into array
    eprint("\n")
    return ccr_a


def valid_ccr(row):
    """ Select ccr according to several conditions.
    :row: row of the dataframe
    :return boolean "False" for not valid ccr
            or boolean "True" for valid ccr
    """
    # Regular expression of the left dangling CIGAR:
    dangling_left = re.compile('^\d+S\d+M$')
    # Regular expression of the rigth dangling CIGAR:
    dangling_right = re.compile('^\d+M\d+S$')
    if (row['chr_donor']!=row['chr_acceptor']
            or row['strand_donor']!=row['strand_acceptor']):
        return False
    elif abs(row['pos_donor'] - row['pos_acceptor']) > 10000000:
        return False
    else:
        if row['strand_donor'] == '-' :
            if (row['pos_donor']>row['pos_acceptor']
                    or row['pos_s1']>row['pos_s2']):
                return False
            else:
                if (not dangling_left.match(row['CIGAR_s1'])
                        or not dangling_right.match(row['CIGAR_s2'])):
                    return False
                else:
                    return True
        elif row['strand_donor'] == '+' :
            if (row['pos_donor']<row['pos_acceptor']
                    or row['pos_s1']<row['pos_s2']):
                return False
            else:
                if (not dangling_right.match(row['CIGAR_s1'])
                        or not dangling_left.match(row['CIGAR_s2'])):
                    return False
                else:
                    return True


def get_exact_circrnas(circ_junctions, cr_threshold):
    """
    Detects circ RNAs as CR sharing exactly the same circ signature
    chrom start end strand
    """
    eprint("Exact circular RNA detection")
    circ_rna_dict = defaultdict(list)
    for circ_junc in circ_junctions:
        circ_rna_dict[circ_junc.key].append(circ_junc)

    circ_array = []
    for key, cr_array in circ_rna_dict.items():
        if len(cr_array) >= cr_threshold:
            circ_array.append(circ.CircRNA(args.output_file_format, ccr_array=cr_array))
    return circ_array


def same_circ_signature(cr_a, cr_b, tol):
    """ Two-by-two comparison of chimeric reads
    :cr_a: chimeric read a
    :cr_b: chimeric read b
    :tol: number of different start-end positions tolerated
    :return boolean "False" for not valid conditions
    and "True" for valid conditions
    """
    if cr_a.chrom != cr_b.chrom:
        return False
    if abs(cr_a.start - cr_b.start) > tol:
        return False
    if abs(cr_a.end - cr_b.end) > tol:
        return False
    if cr_a.read_type == cr_b.read_type:
        if cr_a.strand != cr_b.strand:
            return False
    else:
        if cr_a.strand == cr_b.strand:
            return False
    return True


def merge_circ_rnas(component):
    if len(component) == 1:
        return component[0]
    elements = []
    for circ_rna in component:
        elements.extend(circ_rna.ccr_array)
    return circ.CircRNA(args.output_file_format, ccr_array=elements)


def connected_components(circ_a, neighbours, cr_threshold):
    """ Generate an undirected graph and
        cluters ccr into connected components
    :ccr_a: array of ccr object
    :neighbours: dictionary of neighbours chimeric reads
    :return an array of circular RNAs
    """
    undirected = netwx.Graph()
    undirected.add_nodes_from(range(len(circ_a)))
    undirected.add_edges_from(neighbours)
    circ_rnas = []  # Circular RNAs array
    # Parsing of connected components composed of chimeric reads:
    for comp in netwx.connected_components(undirected):
        circ_rna = merge_circ_rnas([circ_a[i] for i in comp])
        if circ_rna.nb_ccr >= cr_threshold:
            circ_rnas.append(circ_rna)
    return circ_rnas


def cumulative_cr(circ_a):
    return sum([c.nb_ccr for c in circ_a])


def compute_fuzzy_circ_rnas(ccr_a, cr_threshold, tolerance):
    """ Get the circular RNA corresponding to each
        connected component of chimeric reads
    :ccr_a: array of ccr object
    :return array of circular RNAs
    """
    if cumulative_cr(ccr_a) < cr_threshold:
        return []
    n = len(ccr_a)
    # Neighbour cr share sign the same circular RNA
    neighbours = []
    circ_rnas = []
    for i in range(n):
        for j in range(i+1, n):
            if ccr_a[j].start > ccr_a[i].end:
                break
            # Two-by-two comparison of chimeric reads:
            if same_circ_signature(ccr_a[i], ccr_a[j], tolerance):
                neighbours.append((i, j))
    circ_rnas = connected_components(ccr_a, neighbours, cr_threshold)
    return circ_rnas


def get_independent_intervals(circ_rnas):
    # Get independant intervals
    independant_intervals = []
    circ_interval = []
    chrom = None
    for circ_rna in circ_rnas:
        if not chrom:
            chrom = circ_rna.chrom
            max_end = circ_rna.end
        if circ_rna.chrom != chrom:
            independant_intervals.append(circ_interval)
            circ_interval = []
            chrom = circ_rna.chrom
            max_end = circ_rna.start
        if circ_rna.chrom != chrom or circ_rna.start > max_end:
            independant_intervals.append(circ_interval)
            circ_interval = []
            chrom = circ_rna.chrom
            max_end = circ_rna.start
        max_end = max(max_end, circ_rna.end)
        circ_interval.append(circ_rna)
    independant_intervals.append(circ_interval)
    return independant_intervals


def get_fuzzy_circrnas(circ_rnas, cr_threshold, tolerance):
    eprint("Fuzzy circular RNA detection")
    circ_rnas = natsort.natsorted(circ_rnas, key=lambda e: (e.chrom,
                                                            e.start))
    intervals = get_independent_intervals(circ_rnas)
    circ_rnas = []
    for interval in intervals:
        fuzzy_circ_rnas = compute_fuzzy_circ_rnas(interval, cr_threshold,
                                                  tolerance)
        if fuzzy_circ_rnas:
            circ_rnas.extend(fuzzy_circ_rnas)
    return circ_rnas


def circrna_detection(circ_juntcions, cr_threshold, tolerance):
    """ Detecte circular RNAs
    :df_ccr: dataframe of circular chimeric reads
    :return the array of circular RNAs detected
    """
    # First pass, detect exact circular RNAs (no tolerance)
    circ_rnas = get_exact_circrnas(circ_juntcions, cr_threshold)
    eprint("nb of tol %d circRNAS %d" % (tolerance, len(circ_rnas)))#########
    # if tolerance > 0 compute fuzzy circrnas from previous list
    if tolerance > 0:
        circ_rnas = get_fuzzy_circrnas(circ_rnas, cr_threshold, tolerance)
    else:
        circ_rnas = [c for c in circ_rnas if c.nb_ccr >= cr_threshold]
    return circ_rnas


def write_tab_circ_rnas(circ_rnas, outputfile, fmt):
    """ Write the table of circ RNA retained
    :circ_rnas: array of circular RNAs object
    :outputfile: the output file path
    :return the table of circular RNAs retained
    :start Chimeric.out.junction file: last intron base in
    the reference genome order in the circular transcript
    :end Chimeric.out.junction file: first intron base in the
    reference genome order in the circular transcript
    :start GTF: first exon position
    :end GTF: last exon position
    """
    nb_line = 0 # name bed file
    circ_rnas = natsort.natsorted(circ_rnas, key=lambda e: (e.chrom, e.start))
    eprint(" "+"\n"+"Writing the table of circular RNAs..."+"\n")
    with open(outputfile, "w") as fout:
        eprint("\nConvertion the Chimeric.out.junction format to the %s format...\n" % fmt)
        for circ_rna in circ_rnas:
            # Convert Chimeric.out.junction format to out format:
            circ_rna.start += 1
            circ_rna.end -= 1
            # Convert out format to the chosen format:
            circ_rna.write_annot(nb_line, fout, fmt)
            nb_line += 1


def stats(min_cr, tol, nb_ccr_tot, nb_ccr, nb_circ_rna_detected, circ_rnas):
          nb_tot_cr_circ = sum([c.nb_ccr for c in circ_rnas])
          return "%d\t%d\t%d\t%d\t%d\t%d" % (min_cr, tol, nb_ccr_tot, nb_ccr,
                                             nb_circ_rna_detected, nb_tot_cr_circ)


def stats2(df):
    return df.groupby('nb_ccr')['nb_ccr'] \
           .count().reset_index(name='nb_circ_rnas')


def main(r1_input_file, r2_input_file, cr_threshold, tolerance,
         output_file):

    eprint("\nReading junction files...\n")
    df_junctions_a = []
    if r1_input_file:
        df_junctions_a.append(read_junction_file(r1_input_file, read_type="R1"))
    if r2_input_file:
        df_junctions_a.append(read_junction_file(r2_input_file, read_type="R2"))

    eprint("Merging junction files...\n")
    df = merge_junctions(df_junctions_a)
    eprint("%d chimeric reads \n" % len(df.index))

    # Filtering data:
    eprint("Selecting valid ccr...\n")
    circ_junctions = get_valid_circjunctions(df)
    eprint("%d circular chimeric reads" % len(circ_junctions))

    # Circular RNAs detection:
    circ_rnas = circrna_detection(circ_junctions, cr_threshold, tolerance) #### 'circRNA' object has no attribute 'ccr_array'
    eprint("%d circRNAs detected" % len(circ_rnas))

    #for cir in circ_rnas:
    #    print(cir.var_info_str())

    # Write the table of circular RNAs retained:
    circ.write_annotation(circ_rnas, output_file, fmt=args.output_file_format,
                          attributes=["nb_ccr", "genomic_size", "left",
                                      "right", "complete", "nb_distinct"])


def parse_arguments():
    parser = argparse.ArgumentParser(description='Identification of chimeric reads and circular RNAs')
    parser.add_argument('-r1', '--r1_input_file',
                        required=False, help='R1 input file name')
    parser.add_argument('-r2', '--r2_input_file',
                        required=False, help='R2 input file name')
    parser.add_argument('-min_cr', '--cr_threshold', type=int,
                        required=False, default=5,
                        help='Minimum number of chimeric reads')
    parser.add_argument('-tol', '--tolerance', type=int, required=False,
                        default=3,
                        help='Tolerance of different positions to start and end')
    parser.add_argument('-fmt', '--output_file_format',
                        required=False, default='gtf', choices=["bed","gtf"], # bed ou gtf
                        help='Output file format : 0-based (bed) or 1-based (gtf)')
    parser.add_argument('-o', '--output_file', required=False, default='circ_rnas_detected.out',
                        help='Output file path') # gft par défaut # attribut : nb ccr, min_cr_distinct
    parser.add_argument('--verbose', help='Print more info', action='store_true')

    args = parser.parse_args()

    dsp_control = circ.DisplayControl.instance()
    dsp_control.verbose = args.verbose

    if not (args.r1_input_file or args.r2_input_file):
        parser.error('At least one input file is requested')
    return args


if __name__ == '__main__':
    args = parse_arguments()
    main(args.r1_input_file, args.r2_input_file,
         args.cr_threshold, args.tolerance,
         args.output_file)

from collections import defaultdict
import operator, pandas, re
import pybedtools as pybed
import natsort
from tqdm import tqdm
import pandas as pd


from Singleton import Singleton

@Singleton
class DisplayControl():
     def __init__(self):
         self.verbose = False


def tqdm_wrapper(df):
    display_control = DisplayControl.instance()
    if display_control.verbose:
        return tqdm(df.iterrows(), total=df.shape[0])
    else:
        return df.iterrows()


class CircJunction():
    def __init__(self, chrom, start, end, transcription_strand):
        """Initialize parameters of the CCR object"""
        self.chrom = chrom
        self.start = start
        self.end = end
        self.transcript_strand = transcription_strand

    @property
    def key(self):
        return "-".join(map(str,[self.chrom, self.start, self.end,
                                 self.transcript_strand]))


class CCR(CircJunction):
    def __init__(self, record):
        """
        Initialize parameters of the CCR object from junction record of the
        data frame
        """
        chrom = record['chr_donor']
        start, end, transcript_strand = self._setjunction(record)
        super(CCR, self).__init__(chrom, start, end, transcript_strand)
        self.record = record

    @property
    def CIGAR_s1(self):
        return self.record['CIGAR_s1']

    @property
    def CIGAR_s2(self):
        return self.record['CIGAR_s2']

    @property
    def read_type(self):
        return self.record['read_type']

    def _setjunction(self, record):
        strand = record['strand_donor']  # == record['strand_acptor']
        if strand == "+":
            start = record['pos_acceptor']
            end = record['pos_donor']
        else:
            start = record['pos_donor']
            end = record['pos_acceptor']
        if record['read_type'] == "R2":
            transcript_strand = strand
        else:
            transcript_strand = complement(strand)
        return start, end, transcript_strand


class Annotation():

    # initialise class variable
    counter = 0

    def __init__(self, chrom, start, end, strand, feature, attributes, source,
                 name=None):
        Annotation.counter += 1
        self.chrom = chrom
        self.start = int(start)
        self.end = int(end)
        self.strand = strand
        self.feature = feature
        self.source = source
        self.attributes = attributes
        if name:
            self.name = name
        else:
            self.name = "circRNA_%d" % Annotation.counter
        if isinstance(attributes, dict):
            if "name" not in attributes:
                attributes["name"] = name
            self.attributes_d = attributes
            self.attribues = "; ".join(["%s = %s " %(key, value) for key, value in attributes.items()])
        else:
            self.attributes += "; name=%s" % str(self.name)
            self.attributes_d = dict(re.split(" |=", item) for item in self.attributes.split("; "))
            self.clean_quotes()


    def _to_gtf_record(self, attributes=["name"]):
        s = [
            self.chrom, self.source, self.feature,
            self.start, self.end, self._get_score(),
            self.strand, self._get_frame(),
            self._get_attributes(attributes)
        ]
        return "\t".join(map(str, s))


    def _to_bed_record(self, attributes=["name"]):
        start_b, end_b = (self.start - 1, self.end)
        s = [
            self.chrom, start_b, end_b,
            self.name, self._get_score(),
            self.strand, self.source,
            self.feature, self._get_phase(),
            self._get_attributes(attributes)
        ]
        return "\t".join(map(str, s))

    def topybed(self):
        start_b, end_b = (self.start - 1, self.end)
        record = [self.chrom, start_b, end_b, self.name, ".", self.strand]
        return pybed.create_interval_from_list(record)

    def topybed_ife(self):
        record = [self.chrom, self.start, self.end, self.name, ".", self.strand]
        return pybed.create_interval_from_list(record)

    def _get_attributes(self, attributes):
        attribute_str = "; ".join([ "%s=%s" % (key, self.get_att(key)) for key in attributes]) + ";"
        return attribute_str


    @property
    def biotype(self):
        long_biotype = self.get_att("transcript_biotype")
        biotype_keys = ["protein_coding", "pseudogene", "snRNA", "snoRNA", "miRNA", "lncRNA"]
        biotype_values = ["c", "pseudo", "snc", "sno", "mi", "lnc"]
        succint_biotype = dict(zip(biotype_keys, biotype_values))
        if long_biotype in succint_biotype:
            return succint_biotype[long_biotype]
        else:
            # ribozyme, Y_RNA, scaRNA
            return long_biotype # "misc" = miscRNA: Miscellaneous RNA. A non-coding RNA that cannot be classified.

    @property
    def length(self):
        return self.end - self.start + 1

    def _get_score(self):
        return "."

    def _get_frame(self):
        return "."

    def _get_phase(self):
        return "."

    def get_att(self, key):
        if key not in self.attributes_d:
            print("Warning %s is not a valid attribute" % key)
            exit(1)
        return self.attributes_d[key]

    def clean_quotes(self):
        d = self.attributes_d
        for key in d:
            d[key] = d[key].replace('\"', '')

    def get_annotation_tag(self, end_type):
        return "_".join([self.get_att("exon_id"), end_type, self.biotype, self.strand])

    def annot_str(self, fmt="gtf", attributes=["name"]):
        if fmt=="gtf":
            return self._to_gtf_record(attributes)
        elif fmt=="bed":
            return self._to_bed_record(attributes)
        else:
            return "ERROR: Unknown format..."

    def __str__(self):
        return self.annot_str()


class CircRNA(Annotation):
    def __init__(self, fmt, *args, **kwargs):
        """ Initialize parameters of the circRNA object
        :ccr_array: array of circular chimeric reads
        """
        for key, value in kwargs.items():###
            if "ccr_array" in kwargs:
                self._init_from_ccr(kwargs.get("ccr_array"))
            if len(args) > 0:
                self._init_default(fmt, args[0], args[1], args[2], args[3], args[4], args[5], name=value) ###
            self.start_annotation = []
            self.end_annotation = []
            self.start_transcript_id = []
            self.end_transcript_id = []
            self.start_gene_id = []
            self.end_gene_id = []
            self.intron_annotation = []
            self.infra_exonic_annotation = []


    def _init_from_ccr(self, ccr_array):
        # Number of distinct CIGARs:
        self.ccr_array = ccr_array  # Array of cicular chimeric reads
        #rep_ccr = self._consensus_start_end(ccr_array)
        rep_ccr = ccr_array[0]
        chrom = rep_ccr.chrom
        start, end, start_var, end_var = self._compute_consensus_start_end(ccr_array)
        start = rep_ccr.start + 1 ###
        end = rep_ccr.end - 1 ###
        strand = rep_ccr.transcript_strand
        feature = "circRNA"
        source = "circRNA_detection"
        attributes = "nb_ccr=%d" % self._compute_nb_cr(ccr_array)
        super(CircRNA, self).__init__(chrom, start, end, strand, feature,
                                      attributes, source) # name=None
        self.source = source
        self.feature = feature


        self.CIGAR_s1 = rep_ccr.CIGAR_s1  # CIGAR of the first segment
        self.CIGAR_s2 = rep_ccr.CIGAR_s2  # CIGAR of the second segment
        self.read_type = rep_ccr.read_type  # Read type (R1 or R2)

        self.attributes_d['nb_ccr'] = self._compute_nb_cr(ccr_array)
        self.attributes_d['nb_distinct'] = compute_distinct_cr(ccr_array)
        self.attributes_d['genomic_size'] = self._compute_genomic_size()
        self.attributes_d['left'] = self._compute_left()
        self.attributes_d['right'] = self._compute_right()
        self.attributes_d['complete'] = self._compute_complete()

        self.start_var = start_var
        self.end_var = end_var

    @property
    def nb_ccr(self):
        return self.get_att("nb_ccr")

    @property
    def nccr(self):
        return int(self.get_att("nb_ccr"))

    @property
    def complete(self):
        return self.get_att("complete")

    @property
    def left(self):
        return self.get_att("left")

    @property
    def right(self):
        return self.get_att("right")

    @property
    def nb_distinct(self):
        return self.get_att("nb_distinct")

    def _get_score(self):
        return self.nb_ccr


    def _init_default(self, chrom, start, end, strand, feature, attributes,
                      source, name):
        super(CircRNA, self).__init__(chrom, start, end, strand,
                                      feature, attributes, source, name)

    @property
    def read_names(self):
        """ Return an array of chimeric reads names """
        reads_str = []
        for cr in self.ccr_array:
            reads_str.append(cr.name)
        return reads_str


    def _compute_nb_cr(self, ccr_array):
        return len(self.ccr_array)

    def _compute_genomic_size(self):
        """ Compute the genomic size """
        return self.end - self.start + 1

    def _compute_left(self):
        """ Compute the coordinates of genomic regions to
            left [chr:start strand] format
        """
        return str(self.chrom)+":"+str(self.start)+self.strand

    def _compute_right(self):
        """ Compute the coordinates of genomic regions to
            right [chr:end strand] format
        """
        return str(self.chrom)+":"+str(self.end)+self.strand

    def _compute_complete(self):
        """ Compute the coordinates of genomic regions to
            complete [chr:start-end strand] format
        """
        return str(self.chrom)+":"+str(self.start)+":"+str(self.end)+self.strand


    def _compute_consensus_start_end(self, ccr_array):
        """
           Some documentation would be welcome
        """
        d_start = defaultdict(list)
        d_end = defaultdict(list)
        for ccr in ccr_array:
            d_start[ccr.start].append(ccr)
            d_end[ccr.end].append(ccr)
        dd_start = defaultdict(list)
        dd_end = defaultdict(list)
        for key_start, value_start in d_start.items():
            dd_start[key_start] = compute_distinct_cr(value_start)
        for key_end, value_end in d_end.items():
            dd_end[key_end] = compute_distinct_cr(value_end)
        cons_start = max(dd_start.items(), key=operator.itemgetter(1))[0]
        cons_end = max(dd_end.items(), key=operator.itemgetter(1))[0]
        # Now variability arround consensus start and end
        start_var = compute_variability(dd_start, cons_start)
        end_var = compute_variability(dd_end, cons_end)
        # for cr in ccr_array: I don't understand the reason for this loop
        #     cr.start = max_start
        #     cr.end = max_end
        #     return cr
        return cons_start, cons_end, start_var, end_var

    def var_info_str(self):
        starts = ",".join(map(str,[cr.start for cr in self.ccr_array]))
        ends =  ",".join(map(str,[cr.end for cr in self.ccr_array]))
        return "\t".join(map(str, [self.chrom, self.start, self.end, self.nb_ccr,
                                   self.start_var,
                                   self.end_var]))

    def str4diffexp(self):
        key = "%s:%d-%d:%s" % (self.chrom, self.start, self.end, self.strand)
        return "%s\t%d" %(key, self.nccr)

    def add_start_annotation(self, annot):
        self.start_annotation.extend(annot)

    def add_end_annotation(self, annot):
        self.end_annotation.extend(annot)

    def add_start_transcript_id(self, annot):
        self.start_transcript_id.extend(annot)

    def add_end_transcript_id(self, annot):
        self.end_transcript_id.extend(annot)

    def add_start_gene_id(self, annot):
        self.start_gene_id.extend(annot)

    def add_end_gene_id(self, annot):
        self.end_gene_id.extend(annot)

    def get_start_annotation_str(self):
        annotations = []
        for exon_tuple in self.start_annotation:
            annotations.append(exon_tuple[0].get_annotation_tag(exon_tuple[1]))
        return ",".join(annotations)

    def get_end_annotation_str(self):
        annotations = []
        for exon_tuple in self.end_annotation:
            annotations.append(exon_tuple[0].get_annotation_tag(exon_tuple[1]))
        return ",".join(annotations)

    def get_start_transcript_id_str(self):
        annotations = []
        for transcript_tuple in self.start_transcript_id:
            annotations.append(transcript_tuple[0].attributes_d["transcript_id"])
        return ",".join(annotations)

    def get_end_transcript_id_str(self):
        annotations = []
        for transcript_tuple in self.end_transcript_id:
            annotations.append(transcript_tuple[0].attributes_d["transcript_id"])
        return ",".join(annotations)

    def get_start_gene_id_str(self):
        annotations = []
        for gene_tuple in self.start_gene_id:
            annotations.append(gene_tuple[0].attributes_d["gene_id"])
        return ",".join(annotations)

    def get_end_gene_id_str(self):
        annotations = []
        for gene_tuple in self.end_gene_id:
            annotations.append(gene_tuple[0].attributes_d["gene_id"])
        return ",".join(annotations)

    def get_intron_annot_str(self, att):
        annotations = []
        for intron, p_overlap in self.intron_annotation:
            upstream = intron.get_att("up_exon").get_att(att)
            downstream = intron.get_att("down_exon").get_att(att)
            annotations.append("%s-%s" % (upstream, downstream))
        return ",".join(annotations)


    def get_intron_annot_gene_str(self, att):
        annotations = []
        for intron, p_overlap in self.intron_annotation:
            upstream = intron.get_att("up_exon").get_att(att)
            strand_up = intron.get_att("up_exon").strand
            biotype_up = intron.get_att("up_exon").biotype
            downstream = intron.get_att("down_exon").get_att(att)
            strand_d = intron.get_att("down_exon").strand
            biotype_d = intron.get_att("down_exon").biotype
            annotations.append("%s_%s_%s-%s_%s_%s" % (upstream, strand_up, biotype_up,
                                                      downstream, strand_d, biotype_d))
        return ",".join(annotations)


    def annotate_intron(self, intron, p_overlap):
        self.intron_annotation.append((intron, p_overlap))


    def get_infra_exonic_annot_str(self, att):
        annotations = []
        for exon in self.infra_exonic_annotation:
            id_ife = exon.get_att(att)
            annotations.append("%s" % (id_ife))
        return ",".join(annotations)


    def get_infra_exonic_gene_annot_str(self, att):
        annotations = []
        for exon in self.infra_exonic_annotation:
            id_ife = exon.get_att(att)
            strand = exon.strand
            biotype = exon.biotype
            annotations.append("%s_%s_%s" % (id_ife, strand, biotype))
        return ",".join(annotations)


    def get_infra_exonic_exons_start_end(self):
        annotations = []
        for exon in self.infra_exonic_annotation:
            start_exon_ife = exon.start
            end_exon_ife = exon.end
            annotations.append("%s_%s" % (start_exon_ife, end_exon_ife))
        return ",".join(annotations)


    def annotate_ife(self, exon):
        self.infra_exonic_annotation.append(exon)


def complement(strand):
    """ molecular biology complement """
    return "-" if strand == "+" else "+"


def compute_distinct_cr(ccr_array):
    """ From a list of CR compute the number of distincts CR
    Read type with different alignments and different CIGAR
    :ccr_array: array of circular chimeric reads
    """
    cigar_dict = defaultdict(list)
    for cr in ccr_array:
        key = "-".join(map(str, [cr.read_type, cr.CIGAR_s1]))
        cigar_dict[key].append(cr)
    return len(cigar_dict.keys())

def compute_variability(d_values, center):
    deviations = []
    occurences = 0
    for pos, occ in d_values.items():
        deviations.append(occ * abs(pos - center))
        occurences += occ
    return sum(deviations)/occurences


def static_vars(**kwargs):
    def decorate(func):
        for k in kwargs:
            setattr(func, k, kwargs[k])
        return func
    return decorate


@static_vars(counter=0, introns_key_dict=defaultdict())
def get_intron_name(chrom, start, end):
    key = "-".join(map(str, [chrom, start, end]))
    if key not in get_intron_name.introns_key_dict:
        get_intron_name.counter += 1
        name = "INTRON_%d" % get_intron_name.counter
        get_intron_name.introns_key_dict[key] = name
    return get_intron_name.introns_key_dict[key]


@static_vars(introns_name_dict=defaultdict())
def new_intron(upstream_exon, donwstreem_exon):
    chrom = upstream_exon.chrom
    start_i = int(upstream_exon.end) + 1
    end_i = int(donwstreem_exon.start) - 1
    strand = upstream_exon.strand
    name = get_intron_name(chrom, start_i, end_i)
    if name in new_intron.introns_name_dict:
        return None
    new_intron.introns_name_dict[name] = 1
    source = "exon_annot"
    feature = "intron"
    attributes = {"up_exon": upstream_exon, "down_exon":  donwstreem_exon}
    intron = Annotation(chrom, start_i, end_i, strand,
                        feature, attributes, source, name=name)
    return intron


def compute_intronic_positions(d_transcripts):
    """
    returns an array of intron Annotation objects
    """
    introns_a = []
    for transcript_exons in d_transcripts.values():
        # sort by transcript
        sorted_exons = natsort.natsorted(transcript_exons,
                                         key=lambda e: (e.chrom, e.start))
        previous = None
        for exon in sorted_exons:
            if previous:
                intron = new_intron(previous, exon)
                if intron is not None:
                    introns_a.append(intron)
            previous = exon
    return introns_a


def write_annotation(annot_a, outputfile, fmt, attributes):
    with open(outputfile, "w") as fout:
        for record in annot_a:
            fout.write(record.annot_str(fmt, attributes)+"\n")


def read_annotation(file, fmt):
    if fmt=="gtf":
        header = ["chrom", "source", "feature", "start", "end",
                  "score", "strand", "phase", "attributes"]
    elif fmt=="bed":
        header = ["chrom", "start", "end", "gene_id", "score",
                  "strand", "source", "feature", "phase", "attributes"]
    else:
        print("Error: Unknown format...")
    return read_annotation_from_file(file, header)


def read_annotation_from_file(file, header):
    annot_df = pd.read_table(file, sep = '\t', names=header,
                             dtype={"chrom":str, "start":str, "end":str})
    annot_a = []
    for index, row in tqdm_wrapper(annot_df):
        chrom = row.chrom
        start = row.start
        end = row.end
        strand = row.strand
        feature = row.feature
        source = row.source if "source" in header else "."
        attributes = row.attributes
        name = row["gene_id"]
        if feature == "circRNA":
            annot = CircRNA(chrom, start, end, strand, feature, attributes,
                            source, name=name)
        else:
            annot = Annotation(chrom, start, end, strand, feature, attributes,
                               source)
        annot_a.append(annot)
    return annot_a