JoshuaTree2.py
def run_synteny_pipeline(query_prot_id,fasta_path,synteny_path,gff_path, bed_path, gene_info, fasta_out_dir, bps_threshold, loci_threshold, circos, circos_inputs, circos_outputs):
    """Stitch together many pairwise syntenic blocks into multiple species syntenic blocks and output as fasta files for a multiple sequence alignment. If synteny files are not supplied, conduct pairwise synteny between all included strains."""
    query_protID = query_prot_id
    fasta_files = {fasta.split('_')[-2] : fasta for fasta in glob.glob(fasta_path+'/*.fa')+glob.glob(fasta_path+'/*.fasta')}
    gff_files = {gff.split('.')[-2] : gff for gff in glob.glob(gff_path+'/*.gff')+glob.glob(gff_path+'/*.gff3') }
    intersect_keys = set(fasta_files.keys()) & set(gff_files.keys())
    fasta_files = {protID:fasta for protID,fasta in fasta_files.items() if protID in intersect_keys}
    gff_files = {protID:gff for protID,gff in gff_files.items() if protID in intersect_keys}
    genomes = {}
    pairwise_syntenies = []
    for protID in intersect_keys:
        genomes[protID] = Genome(fasta_files[protID],bed_path+'/'+protID+'.bed',protID,gff_files[protID],gene_info)
        if circos:
            genomes[protID].export_karyotype(circos_inputs+'/'+protID+'.karyotype.txt')
    synteny_files = glob.glob(synteny_path+'/*.unout')+glob.glob(synteny_path+'/*.lifted.anchors')
    synteny_protIDs, synteny_protIDs2 = [], []
    remaining_protIDs, remaining_synteny = [], []
    if synteny_files:
        for synteny_file in synteny_files:
            if synteny_file.endswith('.unout'):
                coords = reduce(lambda x,y: x+y, sorted([[m.start(0),m.end(0)] for m in re.finditer('PAC2_0|PAC4GC',synteny_file)]))[1::2]
                q_protID, s_prot_ID = map(lambda x: synteny_file[x+1:x+4],coords)
            else:
                q_protID, s_prot_ID = tuple(synteny_file[synteny_file.rfind('/')+1:].split('.')[:2])
            if q_protID == query_protID or query_protID == 'all': # fixme for now... in future, implement -global option so all sequences can be included
                synteny_protIDs.append((q_protID, s_prot_ID))
                synteny_protIDs2.extend([(q_protID, s_prot_ID),(s_prot_ID, q_protID)])
                pairwise_synteny = PairwiseSynteny(genomes[q_protID],genomes[s_prot_ID],synteny_file,loci_threshold=loci_threshold)
                pairwise_synteny.generate_synteny_structure(synteny_path)
                pairwise_syntenies.append(pairwise_synteny)

        if len(pairwise_syntenies) < len(intersect_keys)-1 and query_protID != 'all':
            synteny_protIDs = set(np.array(synteny_protIDs)[:,1]).union({query_protID})
            remaining_protIDs = set(intersect_keys) - synteny_protIDs
        elif query_protID == 'all':
            remaining_synteny = set(list(combinations(intersect_keys,r=2))) - set(synteny_protIDs)
        else:
            remaining_protIDs = set()
    else:
        if query_protID == 'all':
            remaining_synteny = list(combinations(intersect_keys,r=2))
        else:
            remaining_protIDs = set(intersect_keys) - {query_protID}


    if list(remaining_protIDs) or remaining_synteny:
        def generate_CDS(protID):
            print(protID)
            genomes[protID].extract_CDS()
            return protID
        p = mp.ProcessingPool()
        r = p.amap(generate_CDS,list(set(reduce(lambda x,y: list(x)+list(y),remaining_synteny))) if query_protID == 'all' else remaining_protIDs.union({query_protID}))
        r.wait()
        protIDs = r.get()
        #for protID in remaining_protIDs.union({query_protID}):
        #    genomes[protID].extract_CDS()

        def p_synteny(protIDs):
            q_protID, s_prot_ID = protIDs
            pairwise_synteny = PairwiseSynteny(genomes[q_protID],genomes[s_prot_ID],loci_threshold=loci_threshold)
            pairwise_synteny.generate_synteny_structure(synteny_path)
            return pairwise_synteny

        r = p.amap(p_synteny,remaining_synteny if query_protID == 'all' else [(query_protID,s_prot_ID) for s_prot_ID in remaining_protIDs])#,callback=mycallback) # _async
        r.wait()
        pairwise_syntenies.extend(r.get())
        p.close()

        """
        for s_prot_ID in remaining_protIDs:
            pairwise_synteny = PairwiseSynteny(genomes[query_protID],genomes[s_prot_ID],loci_threshold=loci_threshold)
            pairwise_synteny.generate_synteny_structure(synteny_path)
            pairwise_syntenies.append(pairwise_synteny)"""

    if circos:
        for pairwise_synteny in pairwise_syntenies:
            pairwise_synteny.synteny_structure_2_link(circos_inputs+'/%s.%s.link.txt'%(pairwise_synteny.q_genome.protID,pairwise_synteny.s_genome.protID))
            circos_obj = Circos(pairwise_synteny)
            circos_obj.generate_config(ticks = circos_inputs+'/txticks.conf', ideogram = circos_inputs+'/txideogram.conf', links_and_rules = circos_inputs+'/linksAndrules.conf', config=circos_inputs+'/circos.conf')
            circos_obj.run_circos(circos_outputs+'/')

    super_synteny = SuperSynteny(pairwise_syntenies, bps_threshold, genomes[query_protID])
    super_synteny.generate_global_synteny_graph(fasta_out_dir)
Stitch together many pairwise syntenic blocks into multiple species syntenic blocks and output as fasta files for a multiple sequence alignment. If synteny files are not supplied, conduct pairwise synteny between all included strains.
Stitch together many pairwise syntenic blocks into multiple species syntenic blocks and output as fasta files for a multiple sequence alignment. If synteny files are not supplied, conduct pairwise synteny between all included strains.
Stitch together many pairwise syntenic blocks into multiple species syntenic blocks and output as fasta files for a multiple sequence alignment . If synteny files are not supplied conduct pairwise synteny between all included strains .
def run_synteny_pipeline(query_prot_id,fasta_path,synteny_path,gff_path, bed_path, gene_info, fasta_out_dir, bps_threshold, loci_threshold, circos, circos_inputs, circos_outputs):
    query_protID = query_prot_id
    fasta_files = {fasta.split('_')[-2] : fasta for fasta in glob.glob(fasta_path+'/*.fa')+glob.glob(fasta_path+'/*.fasta')}
    gff_files = {gff.split('.')[-2] : gff for gff in glob.glob(gff_path+'/*.gff')+glob.glob(gff_path+'/*.gff3') }
    intersect_keys = set(fasta_files.keys()) & set(gff_files.keys())
    fasta_files = {protID:fasta for protID,fasta in fasta_files.items() if protID in intersect_keys}
    gff_files = {protID:gff for protID,gff in gff_files.items() if protID in intersect_keys}
    genomes = {}
    pairwise_syntenies = []
    for protID in intersect_keys:
        genomes[protID] = Genome(fasta_files[protID],bed_path+'/'+protID+'.bed',protID,gff_files[protID],gene_info)
        if circos:
            genomes[protID].export_karyotype(circos_inputs+'/'+protID+'.karyotype.txt')
    synteny_files = glob.glob(synteny_path+'/*.unout')+glob.glob(synteny_path+'/*.lifted.anchors')
    synteny_protIDs, synteny_protIDs2 = [], []
    remaining_protIDs, remaining_synteny = [], []
    if synteny_files:
        for synteny_file in synteny_files:
            if synteny_file.endswith('.unout'):
                coords = reduce(lambda x,y: x+y, sorted([[m.start(0),m.end(0)] for m in re.finditer('PAC2_0|PAC4GC',synteny_file)]))[1::2]
                q_protID, s_prot_ID = map(lambda x: synteny_file[x+1:x+4],coords)
            else:
                q_protID, s_prot_ID = tuple(synteny_file[synteny_file.rfind('/')+1:].split('.')[:2])
            if q_protID == query_protID or query_protID == 'all': 
                synteny_protIDs.append((q_protID, s_prot_ID))
                synteny_protIDs2.extend([(q_protID, s_prot_ID),(s_prot_ID, q_protID)])
                pairwise_synteny = PairwiseSynteny(genomes[q_protID],genomes[s_prot_ID],synteny_file,loci_threshold=loci_threshold)
                pairwise_synteny.generate_synteny_structure(synteny_path)
                pairwise_syntenies.append(pairwise_synteny)
        if len(pairwise_syntenies) < len(intersect_keys)-1 and query_protID != 'all':
            synteny_protIDs = set(np.array(synteny_protIDs)[:,1]).union({query_protID})
            remaining_protIDs = set(intersect_keys) - synteny_protIDs
        elif query_protID == 'all':
            remaining_synteny = set(list(combinations(intersect_keys,r=2))) - set(synteny_protIDs)
        else:
            remaining_protIDs = set()
    else:
        if query_protID == 'all':
            remaining_synteny = list(combinations(intersect_keys,r=2))
        else:
            remaining_protIDs = set(intersect_keys) - {query_protID}
    if list(remaining_protIDs) or remaining_synteny:
        def generate_CDS(protID):
            print(protID)
            genomes[protID].extract_CDS()
            return protID
        p = mp.ProcessingPool()
        r = p.amap(generate_CDS,list(set(reduce(lambda x,y: list(x)+list(y),remaining_synteny))) if query_protID == 'all' else remaining_protIDs.union({query_protID}))
        r.wait()
        protIDs = r.get()
        def p_synteny(protIDs):
            q_protID, s_prot_ID = protIDs
            pairwise_synteny = PairwiseSynteny(genomes[q_protID],genomes[s_prot_ID],loci_threshold=loci_threshold)
            pairwise_synteny.generate_synteny_structure(synteny_path)
            return pairwise_synteny
        r = p.amap(p_synteny,remaining_synteny if query_protID == 'all' else [(query_protID,s_prot_ID) for s_prot_ID in remaining_protIDs])
        r.wait()
        pairwise_syntenies.extend(r.get())
        p.close()
    if circos:
        for pairwise_synteny in pairwise_syntenies:
            pairwise_synteny.synteny_structure_2_link(circos_inputs+'/%s.%s.link.txt'%(pairwise_synteny.q_genome.protID,pairwise_synteny.s_genome.protID))
            circos_obj = Circos(pairwise_synteny)
            circos_obj.generate_config(ticks = circos_inputs+'/txticks.conf', ideogram = circos_inputs+'/txideogram.conf', links_and_rules = circos_inputs+'/linksAndrules.conf', config=circos_inputs+'/circos.conf')
            circos_obj.run_circos(circos_outputs+'/')
    super_synteny = SuperSynteny(pairwise_syntenies, bps_threshold, genomes[query_protID])
    super_synteny.generate_global_synteny_graph(fasta_out_dir)
def run_synteny_pipeline ( query_prot_id , fasta_path , synteny_path , gff_path , bed_path , gene_info , fasta_out_dir , bps_threshold , loci_threshold , circos , circos_inputs , circos_outputs ) : query_protID = query_prot_id fasta_files = { fasta . split ( '_' ) [ - 2 ] : fasta for fasta in glob . glob ( fasta_path + '/*.fa' ) + glob . glob ( fasta_path + '/*.fasta' ) } gff_files = { gff . split ( '.' ) [ - 2 ] : gff for gff in glob . glob ( gff_path + '/*.gff' ) + glob . glob ( gff_path + '/*.gff3' ) } intersect_keys = set ( fasta_files . keys ( ) ) & set ( gff_files . keys ( ) ) fasta_files = { protID : fasta for protID , fasta in fasta_files . items ( ) if protID in intersect_keys } gff_files = { protID : gff for protID , gff in gff_files . items ( ) if protID in intersect_keys } genomes = { } pairwise_syntenies = [ ] for protID in intersect_keys : genomes [ protID ] = Genome ( fasta_files [ protID ] , bed_path + '/' + protID + '.bed' , protID , gff_files [ protID ] , gene_info ) if circos : genomes [ protID ] . export_karyotype ( circos_inputs + '/' + protID + '.karyotype.txt' ) synteny_files = glob . glob ( synteny_path + '/*.unout' ) + glob . glob ( synteny_path + '/*.lifted.anchors' ) synteny_protIDs , synteny_protIDs2 = [ ] , [ ] remaining_protIDs , remaining_synteny = [ ] , [ ] if synteny_files : for synteny_file in synteny_files : if synteny_file . endswith ( '.unout' ) : coords = reduce ( lambda x , y : x + y , sorted ( [ [ m . start ( 0 ) , m . end ( 0 ) ] for m in re . finditer ( 'PAC2_0|PAC4GC' , synteny_file ) ] ) ) [ 1 : : 2 ] q_protID , s_prot_ID = map ( lambda x : synteny_file [ x + 1 : x + 4 ] , coords ) else : q_protID , s_prot_ID = tuple ( synteny_file [ synteny_file . rfind ( '/' ) + 1 : ] . split ( '.' ) [ : 2 ] ) if q_protID == query_protID or query_protID == 'all' : synteny_protIDs . append ( ( q_protID , s_prot_ID ) ) synteny_protIDs2 . extend ( [ ( q_protID , s_prot_ID ) , ( s_prot_ID , q_protID ) ] ) pairwise_synteny = PairwiseSynteny ( genomes [ q_protID ] , genomes [ s_prot_ID ] , synteny_file , loci_threshold = loci_threshold ) pairwise_synteny . generate_synteny_structure ( synteny_path ) pairwise_syntenies . append ( pairwise_synteny ) if len ( pairwise_syntenies ) < len ( intersect_keys ) - 1 and query_protID != 'all' : synteny_protIDs = set ( np . array ( synteny_protIDs ) [ : , 1 ] ) . union ( { query_protID } ) remaining_protIDs = set ( intersect_keys ) - synteny_protIDs elif query_protID == 'all' : remaining_synteny = set ( list ( combinations ( intersect_keys , r = 2 ) ) ) - set ( synteny_protIDs ) else : remaining_protIDs = set ( ) else : if query_protID == 'all' : remaining_synteny = list ( combinations ( intersect_keys , r = 2 ) ) else : remaining_protIDs = set ( intersect_keys ) - { query_protID } if list ( remaining_protIDs ) or remaining_synteny : def generate_CDS ( protID ) : print ( protID ) genomes [ protID ] . extract_CDS ( ) return protID p = mp . ProcessingPool ( ) r = p . amap ( generate_CDS , list ( set ( reduce ( lambda x , y : list ( x ) + list ( y ) , remaining_synteny ) ) ) if query_protID == 'all' else remaining_protIDs . union ( { query_protID } ) ) r . wait ( ) protIDs = r . get ( ) def p_synteny ( protIDs ) : q_protID , s_prot_ID = protIDs pairwise_synteny = PairwiseSynteny ( genomes [ q_protID ] , genomes [ s_prot_ID ] , loci_threshold = loci_threshold ) pairwise_synteny . generate_synteny_structure ( synteny_path ) return pairwise_synteny r = p . amap ( p_synteny , remaining_synteny if query_protID == 'all' else [ ( query_protID , s_prot_ID ) for s_prot_ID in remaining_protIDs ] ) r . wait ( ) pairwise_syntenies . extend ( r . get ( ) ) p . close ( ) """
        for s_prot_ID in remaining_protIDs:
            pairwise_synteny = PairwiseSynteny(genomes[query_protID],genomes[s_prot_ID],loci_threshold=loci_threshold)
            pairwise_synteny.generate_synteny_structure(synteny_path)
            pairwise_syntenies.append(pairwise_synteny)""" if circos : for pairwise_synteny in pairwise_syntenies : pairwise_synteny . synteny_structure_2_link ( circos_inputs + '/%s.%s.link.txt' % ( pairwise_synteny . q_genome . protID , pairwise_synteny . s_genome . protID ) ) circos_obj = Circos ( pairwise_synteny ) circos_obj . generate_config ( ticks = circos_inputs + '/txticks.conf' , ideogram = circos_inputs + '/txideogram.conf' , links_and_rules = circos_inputs + '/linksAndrules.conf' , config = circos_inputs + '/circos.conf' ) circos_obj . run_circos ( circos_outputs + '/' ) super_synteny = SuperSynteny ( pairwise_syntenies , bps_threshold , genomes [ query_protID ] ) super_synteny . generate_global_synteny_graph ( fasta_out_dir )
Stitch together many pairwise syntenic blocks into multiple species syntenic blocks and output as fasta files for a multiple sequence alignment.
Stitch together many pairwise syntenic blocks into multiple species syntenic blocks and output as fasta files for a multiple sequence alignment .
"""Stitch together many pairwise syntenic blocks into multiple species syntenic blocks and output as fasta files for a multiple sequence alignment. If synteny files are not supplied, conduct pairwise synteny between all included strains.""" # fixme for now... in future, implement -global option so all sequences can be included #for protID in remaining_protIDs.union({query_protID}): #    genomes[protID].extract_CDS() #,callback=mycallback) # _async """
        for s_prot_ID in remaining_protIDs:
            pairwise_synteny = PairwiseSynteny(genomes[query_protID],genomes[s_prot_ID],loci_threshold=loci_threshold)
            pairwise_synteny.generate_synteny_structure(synteny_path)
            pairwise_syntenies.append(pairwise_synteny)"""
query_prot_id
fasta_path
synteny_path
gff_path
bed_path
gene_info
fasta_out_dir
bps_threshold
loci_threshold
circos
circos_inputs
circos_outputs