def compareProteinFeatures(protein_ft,neg_coding_seq,pos_coding_seq):
    ###Parse out ft-information. Generate ft-fragment sequences for querying
    ###This is a modification of the original script from FeatureAlignment but simplified for exon analysis
    protein_ft_unique=[]; new_ft_list = []
    for ft_data in protein_ft:
        ft_name = ft_data.PrimaryAnnot(); domain_seq = ft_data.DomainSeq(); annotation = ft_data.SecondaryAnnot()
        protein_ft_unique.append((ft_name,annotation,domain_seq))
    ###Redundant entries that are class objects can't be eliminated, so save to a new list and eliminate redundant entries
    protein_ft_unique = unique.unique(protein_ft_unique)
    for (ft_name,annotation,domain_seq) in protein_ft_unique:
        ft_length = len(domain_seq)
        new_ft_data = 'null',domain_seq,ft_name,annotation
        new_ft_list.append(new_ft_data)
    new_ft_list = unique.unique(new_ft_list)
    pos_ft = []; neg_ft = []; all_fts = []
    for (pos,seq,ft_name,annot) in new_ft_list:
        if seq in pos_coding_seq:
            pos_ft.append([pos,seq,ft_name,annot]); all_fts.append([pos,seq,ft_name,annot])
        if seq in neg_coding_seq:
            neg_ft.append([pos,seq,ft_name,annot]); all_fts.append([pos,seq,ft_name,annot])
    all_fts = unique.unique(all_fts)
    pos_ft_missing=[]; neg_ft_missing=[]
    for entry in all_fts:
        if entry not in pos_ft: pos_ft_missing.append(entry)
        if entry not in neg_ft: neg_ft_missing.append(entry)
    pos_ft_missing2=[]; neg_ft_missing2=[]
    for entry in pos_ft_missing: entry[1] = ''; pos_ft_missing2.append(entry)
    for entry in neg_ft_missing: entry[1] = ''; neg_ft_missing2.append(entry)
        
    pos_ft_missing2 = unique.unique(pos_ft_missing2)
    neg_ft_missing2 = unique.unique(neg_ft_missing2)  
    return neg_ft_missing2,pos_ft_missing2

def grabRNAIdentifiers(mrna_assignment):
    ensembl_ids=[]; mRNA_ids=[]
    mRNA_entries = string.split(mrna_assignment,' /// ')
    for entry in mRNA_entries:
        mRNA_info = string.split(entry,' // '); mrna_ac = mRNA_info[0]
        if 'ENS' in mrna_ac: ensembl_ids.append(mrna_ac)
        else:
            try: int(mrna_ac[-3:]); mRNA_ids.append(mrna_ac)
            except ValueError: continue
    ensembl_ids = unique.unique(ensembl_ids)
    mRNA_ids = unique.unique(mRNA_ids)
    return ensembl_ids, mRNA_ids

def reformatHeatmapFile(input_file):
    import unique
    export_file=string.replace(input_file,'Clustering-','Input-')
    eo = export.ExportFile(export_file)
    first_row = True
    fn=filepath(input_file)
    for line in open(fn,'rU').xreadlines():
        data = cleanUpLine(line)
        t = string.split(data,'\t')
        if first_row == True:
            if 'column_clusters-flat' not in t:
                array_names = []
                for i in t[2:]:
                    array_names.append(string.replace(i,':','-'))
                    #array_names.append(i)
            elif 'column_clusters-flat' in t:
                array_clusters = t[2:]
                unique_clusters = unique.unique(array_clusters)
                ind=0; headers=[]
                for c in array_clusters:
                    headers.append(c+'::'+array_names[ind])
                    ind+=1
                headers = string.join(['uid']+headers,'\t')+'\n'
                eo.write(headers)
                first_row = False
        else:
            values = string.join([t[0]]+t[2:],'\t')+'\n'
            eo.write(values)
    return export_file, len(unique_clusters)

def reformatHeatmapFile(input_file):
    import unique

    export_file = string.replace(input_file, "Clustering-", "Input-")
    eo = export.ExportFile(export_file)
    first_row = True
    fn = filepath(input_file)
    for line in open(fn, "rU").xreadlines():
        data = cleanUpLine(line)
        t = string.split(data, "\t")
        if first_row == True:
            if "column_clusters-flat" not in t:
                array_names = []
                for i in t[2:]:
                    array_names.append(string.replace(i, ":", "-"))
                    # array_names.append(i)
            elif "column_clusters-flat" in t:
                array_clusters = t[2:]
                unique_clusters = unique.unique(array_clusters)
                ind = 0
                headers = []
                for c in array_clusters:
                    headers.append(c + "::" + array_names[ind])
                    ind += 1
                headers = string.join(["uid"] + headers, "\t") + "\n"
                eo.write(headers)
                first_row = False
        else:
            values = string.join([t[0]] + t[2:], "\t") + "\n"
            eo.write(values)
    return export_file, len(unique_clusters)

def getFeatureIsoformGenomePositions(species,protein_ft_db,mRNA_protein_seq_db,gene_transcript_db,coordinate_type):
    """ Adapted from compareProteinFeatures but for one isoform and returns genomic coordinates for each feature
    This function is designed to export all unique isoforms rather than just comparison isoforms """
    
    import export
    export_file = 'AltDatabase/ensembl/'+species+'/ProteinFeatureIsoform_complete.txt'                
    export_data = export.ExportFile(export_file)

    failed = 0
    worked = 0
    failed_ac=[]
    for gene in protein_ft_db:
        transcript_feature_db={}
        for ft in protein_ft_db[gene]:
            try:
                ft_name = ft.PrimaryAnnot(); annotation = ft.SecondaryAnnot()
                for (mRNA,type) in gene_transcript_db[gene]:
                    try:
                        protein,protein_seq = mRNA_protein_seq_db[mRNA]
                        error = False
                    except Exception:
                        failed_ac.append(mRNA)
                        error = True
                    if error == False:
                        if ft.DomainSeq() in protein_seq:
                            if coordinate_type == 'genomic':
                                pos1 = ft.GenomicStart(); pos2 = ft.GenomicStop()
                            else:
                                pos1 = str(ft.DomainStart()); pos2 = str(ft.DomainEnd())
    
                            ### There are often many features that overlap within a transcript, so consistently pick just one
                            if mRNA in transcript_feature_db:
                                db = transcript_feature_db[mRNA]
                                if (pos1,pos2) in db:
                                    db[pos1, pos2].append([protein,ft_name,annotation])
                                else:
                                    db[pos1, pos2]=[[protein,ft_name,annotation]]
                            else:
                                db={}
                                db[pos1, pos2]=[[protein,ft_name,annotation]]
                                transcript_feature_db[mRNA] = db
                                
                            #values = [mRNA, protein, pos1, pos2,ft_name,annotation]; unique_entries.append(values)
                            worked+=1
            except IOError:
                failed+=1

        for transcript in transcript_feature_db:
            db = transcript_feature_db[transcript]
            for (pos1,pos2) in db:
                db[pos1,pos2].sort() ### Pick the alphabetically listed first feature
                protein,ft_name,annotation = db[pos1,pos2][0]
                values = [transcript, protein, pos1, pos2,ft_name,annotation]
                export_data.write(string.join(values,'\t')+'\n')
                
    export_data.close()
    print failed,'features failed to have corresponding aligned genomic locations out of', worked+failed
    failed_ac = unique.unique(failed_ac)
    print len(failed_ac),'mRNAs without identified/in silico derived proteins'  ### Appear to be ncRNAs without ATGs
    print failed_ac[:20]

def eliminate_redundant_dict_values(database):
    db1={}
    for key in database:
        list = unique.unique(database[key])
        list.sort()
        db1[key] = list
    return db1

def eliminateRedundant(database):
    db1 = {}
    for key in database:
        list = unique.unique(database[key])
        list.sort()
        db1[key] = list
    return db1

def importSplicingEventsToVisualize(eventsToVisualizeFilename):
    splicing_events = []
    ### Import the splicing events to visualize from an external text file (multiple formats supported)
    type = None
    expandedSearch = False
    firstLine = True
    for line in open(eventsToVisualizeFilename, "rU").xreadlines():
        line = cleanUpLine(line)
        t = string.split(line, "\t")
        if firstLine:
            if "junctionID-1" in t:
                j1i = t.index("junctionID-1")
                j2i = t.index("junctionID-2")
                type = "ASPIRE"
                expandedSearch = True
            if "ANOVA" in t:
                type = "PSI"
            elif "independent confirmation" in t:
                type = "confirmed"
                expandedSearch = True
            elif "ANOVA" in eventsToVisualizeFilename:
                type = "ANOVA"
            firstLine = False
        if "|" in t[0]:
            type = "ANOVA"
        if " " in t[0] and ":" in t[0]:
            splicing_events.append(t[0])
        elif type == "ASPIRE":
            splicing_events.append(t[j1i] + " " + t[j2i])
            splicing_events.append(t[j2i] + " " + t[j1i])
        elif type == "ANOVA":
            try:
                a, b = string.split(t[0], "|")
                a = string.split(a, ":")
                a = string.join(a[1:], ":")
                splicing_events.append(a + " " + b)
                splicing_events.append(b + " " + a)
            except Exception:
                pass
        elif type == "PSI":
            try:
                j1, j2 = string.split(t[0], "|")
                a, b, c = string.split(j1, ":")
                j1 = b + ":" + c
                splicing_events.append(j1 + " " + j2)
                splicing_events.append(j2 + " " + j1)
            except Exception:
                # print traceback.format_exc();sys.exit()
                pass
        elif type == "confirmed":
            try:
                event_pair1 = string.split(t[1], "|")[0]
                a, b, c, d = string.split(event_pair1, "-")
                splicing_events.append(a + "-" + b + " " + c + "-" + d)
                splicing_events.append(c + "-" + d + " " + a + "-" + b)
            except Exception:
                pass
    splicing_events = unique.unique(splicing_events)
    return splicing_events, expandedSearch

def sangerImport(parse_sequences):
    """"Sanger center (miRBase) sequence was provided as a custom (requested) dump of their v5 target predictions
    (http://microrna.sanger.ac.uk/targets/v5/), containing Ensembl gene IDs, microRNA names, and putative target
    sequences, specific for either mouse or human. Mouse was requested in late 2005 whereas human in late 2007.
    These same annotation files, missing the actual target sequence but containing an ENS transcript and coordinate
    locations for that build (allowing seqeunce extraction with the appropriate Ensembl build) exist at:
    http://microrna.sanger.ac.uk/cgi-bin/targets/v5/download.pl"""
    
    if species == 'Hs': filename = 'AltDatabase/miRBS/'+species+'/'+'mirbase-v5_homo_sapiens.mirna.txt'; prefix = 'hsa-'
    if species == 'Rn': filename = 'AltDatabase/miRBS/'+species+'/'+'sanger_miR_target_predictions.txt'; prefix = 'rno-'
    if species == 'Mm': filename = 'AltDatabase/miRBS/'+species+'/'+'sanger_miR_target_predictions.txt'; prefix = 'mmu-'    

    print 'parsing', filename; count=0
    fn=filepath(filename); x=1; mir_sequences=[]
    verifyFile(filename,species) ### Makes sure file is local and if not downloads.
    for line in open(fn,'rU').xreadlines():         
        data = cleanUpLine(line)
        t = string.split(data,'\t')
        if x==0: x=1
        else:
            ensembl_geneids=[]
            if species == 'Hs':
                try:
                    mir = t[1]; ens_transcript = t[2]; ensembl_geneid = t[17]; mir_sequences = string.upper(t[14])
                    ensembl_geneids.append(ensembl_geneid)
                except IndexError: print line;kill
            elif species == 'Mm':
                ens_transcript,mir,mir_sequences = t
                if ens_transcript in ens_gene_to_transcript:
                    ensembl_geneids = ens_gene_to_transcript[ens_transcript]; ensembl_geneid = ensembl_geneids[0]
            elif species == 'Rn':
                ensembl_geneid,mir,mir_sequences = t
                mir_sequences = string.lower(mir_sequences); mir = string.replace(mir,'hsa','rno'); mir = string.replace(mir,'mmu','rno')
                ensembl_geneids=[ensembl_geneid]
            geneid_ls=[]
            #mir_sequences = string.replace(mir_sequences,'-',''); mir_sequences = string.replace(mir_sequences,'=','')
            #mir_sequences = string.upper(mir_sequences)
            #if 'GGCTCCTGTCACCTGGGTCCGT' in mir_sequences:
            #print ensembl_geneid, mir; sys.exit()
            for ensembl_geneid in ensembl_geneids:
                if ensembl_geneid in redundant_ensembl_by_build: ###Thus there are redundant geneids
                    geneid_ls += redundant_ensembl_by_build[ensembl_geneid]+[ensembl_geneid]
                else: geneid_ls += [ensembl_geneid]
                if species == 'Hs':
                    if ens_transcript in ens_gene_to_transcript: geneid_ls+= ens_gene_to_transcript[ens_transcript] 
            geneid_ls = unique.unique(geneid_ls)
            if len(geneid_ls) == 1 and geneid_ls[0]=='': null =[] ###not a valid gene
            elif prefix in mir:
                for ensembl_geneid in geneid_ls:
                    if parse_sequences == 'yes':
                        if (mir,ensembl_geneid) in combined_results:
                            mir_sequences = string.replace(mir_sequences,'-',''); mir_sequences = string.replace(mir_sequences,'=',''); count+=1
                            combined_results[(mir,ensembl_geneid)].append(string.upper(mir_sequences))
                    else:
                        if prefix in mir:
                            y = MicroRNATargetData(ensembl_geneid,'',mir,mir_sequences,'mirbase'); count+=1
                            try: microRNA_target_db[mir].append(y)
                            except KeyError: microRNA_target_db[mir] = [y]
    print count, 'miRNA-target relationships added for mirbase'

def importSplicingEventsToVisualize(eventsToVisualizeFilename):
    splicing_events=[]
    ### Import the splicing events to visualize from an external text file (multiple formats supported)
    type = None
    expandedSearch = False
    firstLine = True
    for line in open(eventsToVisualizeFilename,'rU').xreadlines():
        line = cleanUpLine(line)
        t = string.split(line,'\t')
        if firstLine:
            if 'junctionID-1' in t:
                j1i = t.index('junctionID-1')
                j2i = t.index('junctionID-2')
                type='ASPIRE'
                expandedSearch = True
            if 'ANOVA' in t:
                type='PSI'
            elif 'independent confirmation' in t:
                type='confirmed'
                expandedSearch = True
            elif 'ANOVA' in eventsToVisualizeFilename:
                type = 'ANOVA'
            firstLine=False
        if '|' in t[0]:
            type = 'ANOVA'
        if ' ' in t[0] and ':' in t[0]:
            splicing_events.append(t[0])
        elif type=='ASPIRE':
            splicing_events.append(t[j1i] +' '+ t[j2i])
            splicing_events.append(t[j2i] +' '+ t[j1i])
        elif type=='ANOVA':
            try:
                a,b = string.split(t[0],'|')
                a = string.split(a,':')
                a = string.join(a[1:],':')
                splicing_events.append(a +' '+ b)
                splicing_events.append(b +' '+ a)
            except Exception: pass
        elif type=='PSI':
            try:
                j1,j2 = string.split(t[0],'|')
                a,b,c = string.split(j1,':')
                j1 = b+':'+c
                splicing_events.append(j1 +' '+ j2)
                splicing_events.append(j2 +' '+ j1)
            except Exception:
                #print traceback.format_exc();sys.exit()
                pass
        elif type=='confirmed':
            try:
                event_pair1 = string.split(t[1],'|')[0]
                a,b,c,d = string.split(event_pair1,'-')
                splicing_events.append(a+'-'+b +' '+ c+'-'+d)
                splicing_events.append(c+'-'+d +' '+ a+'-'+b)
            except Exception: pass
	else:
	   splicing_events.append(t[0])
    splicing_events = unique.unique(splicing_events)
    return splicing_events,expandedSearch

def writeIsoformFile(isoform_junctions,o):
    for coord in isoform_junctions:
        isoform_junctions[coord] = unique.unique(isoform_junctions[coord])
        if '+' in coord:
            print coord, isoform_junctions[coord] 
    
    if '+' in coord:
        sys.exit()

def eliminateRedundant(database):
    for key in database:
        try:
            list = makeUnique(database[key])
            list.sort()
        except Exception: list = unique.unique(database[key])
        database[key] = list
    return database

 def Coordinates(self):
     x=0; coords=[]
     for i in self.start_set:
         coord = self.Chr()+':'+str(i)+'-'+str(self.end_set[x])
         coords.append(coord)
         x+=1
     coords = unique.unique(coords)
     coords = string.join(coords,'|') ###If multiple coordinates
     return coords

def grabNestedOntologyIDs():
    nested_ontology_tree={}
    for path in path_dictionary:
        parent_ontology_id = path_ontology_db[path]
        child_ontology_list=[]
        for child_path in path_dictionary[path]:
            child_ontology_id = path_ontology_db[child_path]; child_ontology_list.append(child_ontology_id)
        child_ontology_list = unique.unique(child_ontology_list)
        nested_ontology_tree[parent_ontology_id] = child_ontology_list
    return nested_ontology_tree

def findAvailableOntologies(species,mod_types):
    program_type,database_dir = unique.whatProgramIsThis()
    c = GrabFiles(); c.setdirectory('/'+database_dir+'/'+species+'/gene-go'); file_dirs=[]
    for mod in mod_types:
        file_dirs+= c.searchdirectory(mod+'-')
    avaialble_ontologies=[]
    for filedir in file_dirs:
        ontology_type = string.split(filedir,'-')[-1][:-4] ### remove the .txt
        avaialble_ontologies.append(ontology_type)
    avaialble_ontologies = unique.unique(avaialble_ontologies)
    return avaialble_ontologies

def mirandaImport(parse_sequences,force):
    """Miranda data is avaialble from two file types from different websites. The first is human-centric with multi-species
    target alignment information organized by Ensembl gene ID (http://cbio.mskcc.org/research/sander/data/miRNA2003/mammalian/index.html).
    A larger set of associations was also pulled from species specific files (http://www.microrna.org/microrna/getDownloads.do),
    where gene symbol was related to Ensembl gene. Both files provided target microRNA sequence."""
    
    ### Then download the latest annotations and sequences
    if parse_sequences == 'coordinates':
        export_object = export.ExportFile('miRanda/'+species+'/miRanda.txt')
        print "Exporting to:"+'miRanda/'+species+'/miRanda.txt'
    verify, filename = verifyExternalDownload('miRanda')
    if verify == 'no': filename = downloadFile('miRanda')
    print 'parsing', filename; count=0; null_count=[]
    fn=filepath(filename); x=1; mir_sequences=[]
    verifyFile(filename,species) ### Makes sure file is local and if not downloads.
    for line in open(fn,'rU').xreadlines():         
        data = cleanUpLine(line)
        t = string.split(data,'\t')
        if x==0: x=1
        else:
            symbol = string.upper(t[3]); mir = t[1]; entrez_gene = t[2]; mir_sequences = string.upper(t[8])
            mir_sequences = string.replace(mir_sequences,'-',''); mir_sequences = string.replace(mir_sequences,'=','')
            mir_sequences = string.replace(mir_sequences,'U','T')
            #if 'GGCTCCTGTCACCTGGGTCCGT' in mir_sequences:
            #print symbol, mir; sys.exit()
            ensembl_gene_ids = []
            if symbol in symbol_ensembl_current:
                ensembl_gene_ids = symbol_ensembl_current[symbol]
            else: ensembl_gene_ids=[]; null_count.append(symbol)

            if len(ensembl_gene_ids) > 0:
                for ensembl_geneid in ensembl_gene_ids:
                    if parse_sequences == 'yes':
                        if (mir,ensembl_geneid) in combined_results:
                            combined_results[(mir,ensembl_geneid)].append(string.upper(mir_sequences)); count+=1
                    else:
                        y = MicroRNATargetData(ensembl_geneid,'',mir,mir_sequences,'miRanda'); count+=1
                        try: microRNA_target_db[mir].append(y)
                        except KeyError: microRNA_target_db[mir] = [y]
                        if parse_sequences == 'coordinates':
                            """
                            genome_coord = string.split(t[13],':')[1:]; chr = 'chr'+ genome_coord[0]
                            strand = genome_coord[-1]; start, stop = string.split(genome_coord[1],'-')
                            """
                            genome_coord = t[13][1:-1]
                            align_score = t[15]
                            y.setScore(align_score); y.setCoordinates(genome_coord)
                            export_object.write(y.Output())
    print count, 'miRNA-target relationships added for miRanda'
    null_count = unique.unique(null_count)
    print len(null_count), 'missing symbols',null_count[:10]
    if parse_sequences == 'coordinates': export_object.close()

 def Ensembl(self):
     ens_list = unique.unique(self._ensembl)
     ens_str = string.join(ens_list,',')
     return ens_str

def analyzeCommonProteinClassesAndCompartments(sm,kw,ft_call,ft_string,rc,de,go):
    ### Used to assign "Common Protein Classes" annotations to Gene Expression summary file (ExpressionOutput folder)
    class_def=[]; annotation=[]; cellular_components = []
    if 'DNA-binding domain' in sm or 'Transcription' in go or 'Transcription regulation' in kw: class_def.append('transcription regulator')
    if 'protein kinase superfamily' in sm or 'Kinase' in go: class_def.append('kinase')
    if 'mRNA splicing' in kw or 'mRNA processing' in kw: class_def.append('splicing regulator')

    if 'G-protein coupled receptor' in sm or 'LU7TM' in sm:
        g_type = []
        if ('adenylate cyclase' in ft_call) or ('adenylyl cyclase'in ft_call):
                ###if both occur
                if (('stimulat' in ft_call) or ('activat' in ft_call)) and ('inhibit' in ft_call):
                    if 'inhibit aden' in ft_call: g_type.append('Gi')
                    if 'stimulate aden' in ft_call or 'activate aden' in ft_call: g_type.append('Gs')
                elif ('stimulat' in ft_call) or ('activat' in ft_call): g_type.append('Gs')
                elif ('inhibit' in ft_call): g_type.append('Gi')
        if ('cAMP' in ft_call):
            if ('stimulat' in ft_call) or ('activat' in ft_call): g_type.append('Gs')
            if ('inhibit' in ft_call): g_type.append('Gi')
        if ('G(s)' in ft_call): g_type.append('Gs')
        if ('G(i)' in ft_call): g_type.append('Gi')
        if ('pertussis' in ft_call and 'insensitive' not in ft_call): g_type.append('Gi')
        if ('G(i/0)' in ft_call) or ('G(i/o)' in ft_call): g_type.append('Gi')
        if ('G(o)' in ft_call): g_type.append('Go')
        if ('G(alpha)q' in ft_call): g_type.append('Gq')
        if ('G(11)' in ft_call): g_type.append('G11')
        if ('G(12)' in ft_call): g_type.append('G12')
        if ('G(13)' in ft_call): g_type.append('G13')
        if ('mobiliz' in ft_call and 'calcium' in ft_call and 'without formation' not in ft_call): g_type.append('Gq')
        if ('phosphatidyl' in ft_call and 'inositol' in ft_call) or ('G(q)' in ft_call) or ('phospholipase C' in ft_call):
                g_type.append('Gq')
        if ('inositol phos' in ft_call) or ('phosphoinositide' in ft_call) or ('PKC'in ft_call) or ('PLC' in ft_call):
            g_type.append('Gq')
        if ('intracellular' in ft_call and 'calcium' in ft_call) and 'nor induced' not in ft_call: g_type.append('Gq')
        if 'G-alpha-11' in ft_call: g_type.append('G11')
        if 'Orphan' in ft_call or 'orphan' in ft_call: g_type.append('orphan')
        if 'taste' in ft_call or 'Taste' in ft_call: g_type.append('taste')
        if 'vision' in ft_call or 'Vision' in ft_call: g_type.append('vision')
        if 'odorant' in ft_call or 'Odorant' in ft_call: g_type.append('oderant')
        if 'influx of extracellar calcium' in ft_call: g_type.append('Gq')
        if 'pheromone receptor' in ft_call or 'Pheromone receptor' in ft_call: g_type.append('pheromone')
        g_protein_list = unique.unique(g_type); g_protein_str = string.join(g_protein_list,'|')
        class_def.append('GPCR(%s)' % g_protein_str)
    elif 'receptor' in sm or 'Receptor' in go: class_def.append('receptor')
    if len(ft_string)>0: ### Add cellular component annotations
        if 'ecreted' in sm: k = 1; annotation.append('extracellular')
        if 'Extraceullar space' in sm: k = 1; annotation.append('extracellular')
        if 'ecreted' in go: k = 1; annotation.append('extracellular')
        if 'xtracellular' in go: k = 1; annotation.append('extracellular')
        if 'Membrane' in sm: k = 1; annotation.append('transmembrane')
        if 'TRANSMEM' in ft_string: k = 1; annotation.append('transmembrane')
        if 'integral to membrane' in go: k = 1; annotation.append('transmembrane')
        if 'Nucleus' in sm: k = 1; annotation.append('nucleus')
        if 'nucleus' in go: k = 1; annotation.append('nucleus')
        if 'Cytoplasm' in sm: k = 1; annotation.append('cytoplasm')
        if 'Mitochondrion' in sm: k = 1; annotation.append('mitochondrion')
        if 'SIGNAL' in ft_string: k = 1; annotation.append('signal')
        ###Generate probably secreted annotations
        if 'signal' in annotation and 'transmembrane' not in annotation:
            for entry in annotation:
                if entry != 'signal': cellular_components.append(entry)
            cellular_components.append('extracellular');annotation = cellular_components
        elif 'signal' in annotation and 'transmembrane' in annotation:
            for entry in annotation:
                if entry != 'signal': cellular_components.append(entry)
            annotation = cellular_components
            
    cellular_components = string.join(annotation,'|')
    class_def = string.join(class_def,'|')
    return class_def, cellular_components

 def test_unique(self):
     self.assertEqual(unique('programmare'), 4)

def associateQueryGenesWithInteractions(query_db,query_interactions,dir_file):
    suffix=''
    if dir_file!=None:
        if len(dir_file)!=0:
            suffix='-'+intNameShort+'_'+export.findFilename(dir_file)[:-4]
    if len(suffix)==0:
        try: suffix = '_'+FileName
        except Exception: None
    file_name = 'AltAnalyze-network'+suffix
    
    query_interactions_unique={}
    interacting_genes={}
    connections = 1
    primary=0
    secondary=0
    terciary=0
    for ensemblGene in query_db:
        if ensemblGene in interaction_db:
            for interacting_ensembl in interaction_db[ensemblGene]:
                if interacting_ensembl not in blackList:
                    ###Only allow direct interactions found in query
                    if interacting_ensembl in query_db:
                        try: query_interactions[ensemblGene].append(interacting_ensembl)
                        except KeyError: query_interactions[ensemblGene] = [interacting_ensembl]
                        try: query_interactions[interacting_ensembl].append(ensemblGene)
                        except KeyError: query_interactions[interacting_ensembl] = [ensemblGene]
                        primary+=1
                    if degrees == 2 or degrees == 'indirect':
                        try: interacting_genes[interacting_ensembl].append(ensemblGene)
                        except KeyError: interacting_genes[interacting_ensembl] = [ensemblGene]
                    elif degrees == 'allInteracting' or degrees == 'all possible':
                        try: query_interactions[ensemblGene].append(interacting_ensembl)
                        except KeyError: query_interactions[ensemblGene] = [interacting_ensembl]
                    if interacting_ensembl in secondaryQueryIDs: ### IDs in the expression file
                        secondary+=1 ### When indirect degrees selected, no additional power added by this (only for direct or shortest path)
                        try: query_interactions[ensemblGene].append(interacting_ensembl)
                        except KeyError: query_interactions[ensemblGene] = [interacting_ensembl]    
        if ensemblGene in second_degree_obligatory:
            for interacting_ensembl in second_degree_obligatory[ensemblGene]:
                try: interacting_genes[interacting_ensembl].append(ensemblGene)
                except KeyError: interacting_genes[interacting_ensembl] = [ensemblGene]

    ### Include indirect interactions to secondaryQueryIDs from the expression file
    if degrees == 2 or degrees == 'indirect':
        for ensemblGene in secondaryQueryIDs:
            if ensemblGene in interaction_db:
                for interacting_ensembl in interaction_db[ensemblGene]:
                    if interacting_ensembl not in blackList:
                        try:
                            interacting_genes[interacting_ensembl].append(ensemblGene)
                            terciary+=1#; print interacting_ensembl
                        except KeyError: None ### Only increase the interacting_genes count if the interacting partner is present from the primary query list
    #print primary,secondary,terciary
    
    ### Report the number of unique interacting genes
    for interacting_ensembl in interacting_genes:
        if len(interacting_genes[interacting_ensembl])==1:
            interacting_genes[interacting_ensembl] = 1
        else:
            unique_interactions = unique.unique(interacting_genes[interacting_ensembl])
            interacting_genes[interacting_ensembl] = len(unique_interactions)
    
    query_indirect_interactions={}; indirect_interacting_gene_list=[]; interacting_gene_list=[]; added=[] 
    if degrees=='shortestPath' or degrees=='shortest path': ### Typically identifying the single smallest path(s) between two nodes.
        query_indirect_interactions, indirect_interacting_gene_list, interacting_gene_list = evaluateShortestPath(query_db,interaction_db,10)
        
    else:
        if degrees==2 or degrees=='indirect' or len(secondDegreeObligatoryCategories)>0:
            for ensembl in interacting_genes:
                if interacting_genes[ensembl] > connections:
                    if ensembl in interaction_db: ### Only nodes removed due to promiscuity will not be found
                        for interacting_ensembl in interaction_db[ensembl]:
                            if interacting_ensembl in query_db or interacting_ensembl in secondaryQueryIDs:
                                try: query_indirect_interactions[interacting_ensembl].append(ensembl)
                                except KeyError: query_indirect_interactions[interacting_ensembl] = [ensembl]
                        ###Record the highest linked nodes
                        indirect_interacting_gene_list.append((interacting_genes[ensembl],ensembl)) 
        if len(obligatory_interactions)>0: ### Include always
            all_reported_genes = combineDBs(query_interactions,query_indirect_interactions) ### combinesDBs and returns a unique list of genes
            for ensemblGene in all_reported_genes: ###This only includes genes in the original input list
                if ensemblGene in obligatory_interactions:
                    for interacting_ensembl in obligatory_interactions[ensemblGene]:
                        #symbol = ensembl_symbol_db[ensemblGene]                    
                        try: query_interactions[ensemblGene].append(interacting_ensembl)
                        except KeyError: query_interactions[ensemblGene] = [interacting_ensembl]
    
    z = dict(query_interactions.items() + query_indirect_interactions.items())
    interaction_restricted_db={}
    for ensembl in z:
        interacting_nodes = z[ensembl]
        for node in interacting_nodes:
            if ensembl in interaction_restricted_db:
                db = interaction_restricted_db[ensembl]
                db[node] = 1
            else: interaction_restricted_db[ensembl] = {node:1}

            if node in interaction_restricted_db:
                db = interaction_restricted_db[node]
                db[ensembl] = 1
            else: interaction_restricted_db[node] = {ensembl:1}
#.........这里部分代码省略.........