@files( importAffymetrixAnnotation, "probeset2transcript.table")
def buildProbeset2Transcript( infile, outfile):
    '''build map relating a probeset to a transcript_id, gene_id, ...

    1. only take those that can be mapped to transcripts in table transcripts
    2. check consistency of location with transcripts in table transcripts
    3. remove those with cross hybridyzation signals

    '''

    dbhandle = sqlite3.connect( PARAMS["database"] )

    E.debug( "collecting exons")
    raise NotImplementError("check fastgtf and gzip.")
    with gzip.open(PARAMS["transcripts"]) as inf:
        exons = GTF.readAndIndex( GTF.iterator( inf ) )

    E.debug( "collecting transcripts")
    cc = dbhandle.cursor()
    statement = "SELECT transcript_id, gene_id, source, contig, min(start), max(end), strand FROM transcripts GROUP BY transcript_id, gene_id, source, contig"
    cc.execute( statement )
    map_transcript2gene = {}
    for transcript_id, gene_id, source, contig, start, end, strand in cc.fetchall():
        map_transcript2gene[transcript_id] = (gene_id, source, contig, start, end, strand )

    counter = E.Counter()

    cc = dbhandle.cursor()
    statement = "SELECT probeset, transcript_cluster_id, seqname, start, stop, mrna_assignment, crosshyb_type FROM affymetrix"
    cc.execute( statement )

    outf_notfound = open( outfile + ".notfound", "w")
    outf_notfound.write("probeset\tcause\tcontig\tstart\tend\tmrna_assignment\n" )

    outf = open( outfile, "w" )
    outf.write( "probeset\tcluster_id\tprobe_contig\tprobe_start\tprobe_end\ttranscript_id\tgene_id\tsource\tcontig\tstart\tend\tstrand\n" )

    genes_found, transcripts_found = set(), set()

    for probeset, transcript_cluster_id, contig, start, end, mrna_assignment, crosshyb_type in cc.fetchall():
        start -= 1 
        if mrna_assignment == "---":
            counter.probesets_unassigned += 1
            continue

        data = [ x.strip().split(" // ") for x in mrna_assignment.split("///") ]
        found, mapped, mismapped, nxhyb = 0, 0, 0, 0

        # collect transcripts by exon overlap
        try:
            overlaps = [ x[2] for x in exons.get( contig, start, end ) ]
        except KeyError:
            overlaps = []

        transcripts = set( [x.transcript_id for x in overlaps ] )

        for transcript_id, seqname, score, direct_probes, possible_probes, xhyb in data:

            if xhyb != "0":
                nxhyb += 1
                continue

            # collect transcripts from annotation file
            if transcript_id.startswith("ENST"):
                found += 1
                    
                try:
                    gene_id, other_source, other_contig, other_start, other_end, other_strand = map_transcript2gene[transcript_id]
                except KeyError:
                    continue

                mapped += 1

                if transcript_id not in transcripts:
                    mismapped += 1
                    continue
                
                if other_contig != contig or \
                        min(other_end,end) - max(other_start, start) < 0:
                    mismapped += 1
                    continue

                outf.write("%s\t%s\t%s\t%i\t%i\t%s\t%s\t%s\t%s\t%i\t%i\t%s\n" % \
                               (probeset, transcript_cluster_id, contig, start, end,
                                transcript_id, gene_id, other_source,
                                other_contig, other_start, other_end, other_strand))

                genes_found.add(gene_id)
                transcripts_found.add(transcript_id)

        if found == 0: 
            counter.probesets_notfound += 1
            code = "no ENSEMBL"
        elif mapped == 0:
            counter.probesets_unmapped += 1
            code = "not in transcripts"
        elif mapped > 0 and mapped == mismapped:
            counter.probesets_mismapped += 1
            code = "mismapped"
        elif mapped > 0 and mapped == nxhyb:
            counter.probesets_crosshyb += 1
            code = "cross-hybridized"
        elif mapped > 0 and mapped == nxhyb + mismapped:
            counter.probesets_errors += 1
            code = "errors"
        else:
            counter.probesets_found += 1
            code = None

        if code:
            outf_notfound.write("%s\t%s\t%s\t%i\t%i\t%s\n" % (probeset, code, contig, start, end, mrna_assignment) )

    outf_notfound.close()
    outf.close()

    E.info( "mapping results: %s" % str(counter))
    E.info( "genes found: %i(%5.2f%%): transcripts found %i(%5.2f%%)" %\
                (len(genes_found),
                 100.0 * len(genes_found) / len( set( [x[0] for x in map_transcript2gene.values()]) ),
                 len(transcripts_found),
                 100.0 * len(transcripts_found) / len( map_transcript2gene.keys() ) ) )

############################################################
############################################################
############################################################
##
############################################################
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_vitaminD.py/buildProbeset2Transcript/original