Remove junctions that do not use an annotated donor or acceptor according
    to external junction annotation. Add strand and gene information for
    junctions according to external annoation (STAR strand ignored).

    Parameters
    ----------
    sj_outP : pandas.Panel
        Panel where each dataframe(?) corresponds to an sj_out file filtered to
        remove low coverage junctions.

    annotDF : pandas.DataFrame
        Annotation information from STAR. Should have the following columns:
        'chrom', 'start', 'end', 'intron_motif',
        'annotated'.

    extDF : pandas.DataFrame
        Dataframe containing information about annotated splice sites. These can
        differ from those provided to STAR for alignment. The dataframe should
        have the following columns: 'gene', 'chrom', 'start', 'end', 'strand',
        'chrom:start', 'chrom:end', 'donor', 'acceptor', 'intron'

    Returns
    -------
    countDF :  pandas.DataFrame
        Number of unique junction spanning reads for each junction that passed
        filtering criteria.

    annotDF : pandas.DataFrame
        Annotation information for junctions that passed filtering criteria.
    
    stats : list of strings
        Human readable statistics about the combining process.

def _filter_jxns_donor_acceptor(sj_outP, annotDF, extDF):
    """Remove junctions that do not use an annotated donor or acceptor according
    to external junction annotation. Add strand and gene information for
    junctions according to external annoation (STAR strand ignored).

    Parameters
    ----------
    sj_outP : pandas.Panel
        Panel where each dataframe(?) corresponds to an sj_out file filtered to
        remove low coverage junctions.

    annotDF : pandas.DataFrame
        Annotation information from STAR. Should have the following columns:
        'chrom', 'start', 'end', 'intron_motif',
        'annotated'.

    extDF : pandas.DataFrame
        Dataframe containing information about annotated splice sites. These can
        differ from those provided to STAR for alignment. The dataframe should
        have the following columns: 'gene', 'chrom', 'start', 'end', 'strand',
        'chrom:start', 'chrom:end', 'donor', 'acceptor', 'intron'

    Returns
    -------
    countDF :  pandas.DataFrame
        Number of unique junction spanning reads for each junction that passed
        filtering criteria.

    annotDF : pandas.DataFrame
        Annotation information for junctions that passed filtering criteria.
    
    stats : list of strings
        Human readable statistics about the combining process.
    
    """
    import re

    sjRE = re.compile('(.*:.*-.*):(\+|-)')
    juncRE = re.compile('(.*):(\d*)-(\d*):') 
    
    # Add column showing whether junction is in external annotation.
    annotDF = copy.deepcopy(annotDF)
    annotDF['ext_annotated'] = False
    annotDF.ix[set(annotDF.index) & set(extDF.intron),
               'ext_annotated'] = True
    
    # Replace strand information from STAR with strand information from external
    # database.
    strandSE = pd.Series(extDF.strand.values,index=extDF.intron)
    strandSE = strandSE[set(strandSE.index) & set(annotDF.index)]
    annotDF['strand'] = '*'
    annotDF.ix[strandSE.index,'strand'] = strandSE.values

    # Add column for start and end location (chromosome plus position for
    # uniqueness).
    annotDF['chrom:start'] = (annotDF.chrom + ':' + 
                              annotDF.start.astype(int).astype(str))
    annotDF['chrom:end'] = (annotDF.chrom + ':' +
                            annotDF.end.astype(int).astype(str))

    ext_startS = set(extDF['chrom:start'].values)
    ext_endS = set(extDF['chrom:end'].values)

    # Remove junctions that don't have a start or end shared with external
    # database.
    tdf = pd.DataFrame(True, columns=['temp_start_in_ext'], index=ext_startS)
    annotDF = annotDF.merge(tdf, left_on='chrom:start', right_index=True,
                            how='left')
    tdf = pd.DataFrame(True, columns=['temp_end_in_ext'], index=ext_endS)
    annotDF = annotDF.merge(tdf, left_on='chrom:end', right_index=True,
                            how='left')
    junctions_to_remove = annotDF[(annotDF.ext_annotated == False) &
                                  (annotDF.temp_start_in_ext.isnull()) &
                                  (annotDF.temp_end_in_ext.isnull())].index

    #junctions_to_removeSE = annotDF[annotDF.ext_annotated == False].apply(
    #    lambda x: (x['chrom:start'] in ext_startS) + 
    #    (x['chrom:end'] in ext_endS) == 0,axis=1)
    if len(junctions_to_remove) > 0:
        annotDF = annotDF.drop(junctions_to_remove)

    # Add column indicating which gene the junctions belong to for annotated
    # jxn's.
    geneSE = pd.Series(dict(zip(extDF.intron.values,extDF.gene)))
    annotDF['gene_id'] = ''
    annotDF['gene_id'] = geneSE[annotDF.index]

    # Now we'll figure out the genes for the junctions that aren't in our
    # database. We can associate each start and end with a gene and use this.
    start_gene = dict(zip(extDF['chrom:start'], extDF.gene))
    end_gene = dict(zip(extDF['chrom:end'], extDF.gene))
    
    genes = []
    for ind in annotDF[annotDF.ext_annotated == False].index:
        cur_start = annotDF.ix[ind,'chrom:start'] 
        cur_end = annotDF.ix[ind,'chrom:end'] 
        gene = start_gene.get(cur_start, '')
        if gene == '':
            gene = end_gene.get(cur_end, '')
        genes.append(gene)
    annotDF.ix[annotDF.ext_annotated == False, 'gene_id'] = genes

    # We can use the genes to assign strand to the novel splice junctions.
    strandSE = pd.Series(dict(zip(extDF.gene,extDF.strand)))
    ind = annotDF[annotDF.ext_annotated == False].index
    annotDF.ix[ind, 'strand'] = strandSE[annotDF.ix[ind, 'gene_id']].values

    # And re-index with the strand info.
    annotDF.index = [ x + ':' + annotDF.ix[x,'strand'] for x in annotDF.index ]

    # Now we'll add donor and acceptor info.
    pos_ind = annotDF[annotDF.strand == '+'].index
    neg_ind = annotDF[annotDF.strand == '-'].index
    annotDF.ix[pos_ind, 'donor'] = (annotDF.ix[pos_ind, 'chrom'] + ':' + 
                                    annotDF.ix[pos_ind, 'start'].astype(str) + 
                                    ':' +  annotDF.ix[pos_ind, 'strand'])
    annotDF.ix[pos_ind, 'acceptor'] = (annotDF.ix[pos_ind, 'chrom'] + ':' + 
                                       annotDF.ix[pos_ind, 'end'].astype(str) 
                                       + ':' + annotDF.ix[pos_ind, 'strand'])
    annotDF.ix[neg_ind, 'acceptor'] = (annotDF.ix[neg_ind, 'chrom'] + ':' + 
                                       annotDF.ix[neg_ind, 'start'].astype(str) 
                                       + ':' + annotDF.ix[neg_ind, 'strand'])
    annotDF.ix[neg_ind, 'donor'] = (annotDF.ix[neg_ind, 'chrom'] + ':' + 
                                    annotDF.ix[neg_ind, 'end'].astype(str) + 
                                    ':' + annotDF.ix[neg_ind, 'strand'])

    # And whether the donor or acceptor is in the external database or not.
    ext_donor = pd.DataFrame(False, index=list(set(extDF.donor)),
                             columns=['novel_donor'])
    ext_acceptor = pd.DataFrame(False, index=list(set(extDF.acceptor)),
                                columns=['novel_acceptor'])
    annotDF = annotDF.merge(ext_donor, left_on='donor', right_index=True,
                            how='left')
    annotDF.ix[annotDF.novel_donor.isnull(), 'novel_donor'] = True
    annotDF = annotDF.merge(ext_acceptor, left_on='acceptor', right_index=True,
                            how='left')
    annotDF.ix[annotDF.novel_acceptor.isnull(), 'novel_acceptor'] = True

    # ext_donorS = frozenset(extDF.donor)
    # ext_acceptorS = frozenset(extDF.acceptor)
    # annotDF['novel_donor'] = False
    # annotDF['novel_acceptor'] = False
    # ind = annotDF[annotDF.ext_annotated == False].index
    # if len(ind) > 0:
    #     novel_donor = []
    #     novel_acceptor = []
    #     for i in ind:
    #         novel_donor.append(annotDF.ix[i, 'donor'] not in ext_donorS)
    #         novel_acceptor.append(annotDF.ix[i, 'acceptor'] not in ext_acceptorS)
    #     annotDF.ix[ind, 'novel_donor'] = novel_donor
    #     annotDF.ix[ind, 'novel_acceptor'] = novel_acceptor

    # Sort by gene ID and start/end.
    annotDF = annotDF.sort_values(by=['gene_id', 'start', 'end'])

    # Make file with counts for the junctions we are interested in.
    L = [ juncRE.match(x).group().strip(':') for x in annotDF.index ]
    countDF = sj_outP.ix[:, L, 'unique_junction_reads']
    countDF.index = annotDF.index

    stats = []
    t = extDF.shape[0]
    stats.append('Junction filtering stats')
    stats.append(('Number of junctions in external '
                  'annotation\t{0:,}').format(t))

    t = annotDF.annotated.sum()
    stats.append(('Number observed junctions in STAR '
                  'annotation\t{0:,}').format(t))

    t = annotDF.ext_annotated.sum()
    stats.append(('Number observed junctions in external '
                  'annotation\t{0:,}').format(t))

    t = annotDF.shape[0] - annotDF.ext_annotated.sum()
    stats.append(('Number of observed junctions not in external '
                  'annotation\t{0:,}').format(t))

    t = annotDF.ix[annotDF.ext_annotated == False,'annotated'].sum()
    stats.append(('Number of observed junctions not in external annotation but '
                  'in STAR annotation\t{0:,}').format(t))

    t = sum(annotDF.ix[annotDF.ext_annotated == False,'annotated'].values == 0)
    stats.append(('Number of observed junctions not in external annotation and '
                  'not in STAR annotation\t{0:,}').format(t))
  
    t = len(set(annotDF.ix[annotDF.novel_donor, 'donor']))
    stats.append(('Number of novel donors\t{0:,}').format(t))

    t = annotDF.novel_donor.sum()
    stats.append(('Number of novel junctions with novel '
                  'donors\t{0:,}').format(t))

    t = len(set(annotDF.ix[annotDF.novel_acceptor, 'acceptor']))
    stats.append(('Number of novel acceptors\t{0:,}').format(t))

    t = annotDF.novel_acceptor.sum()
    stats.append(('Number of novel junctions with novel '
                  'acceptors\t{0:,}').format(t))

    t = (annotDF[annotDF.ext_annotated == False].shape[0] - 
         sum(annotDF.novel_donor) - 
         sum(annotDF.novel_acceptor))
    stats.append(('Number of novel junctions with new combination of donor and '
                  'acceptor\t{0:,}').format(t))
    return countDF, annotDF, stats