<jupyter_start><jupyter_text>Use YGOB and Trevors intergenics to make orthologous promoters .<jupyter_code>#For a given gene: 

#Given a species list from the YGOB Pillars

    #get set of orthologous genes

#make fasta of orthologous promoters. 
#export file. 


#count no of STREs, PDS's.  Output surrounding 2 bases and locations
#Align and get all conserved regions

#visualize list.
#Make a dataframe based on the YGOB pillars file

fname_pillars = os.path.normpath(data_processing_dir + 'ortholog_files_YGOB\YGOB_pillars_bmh.txt')
YGOB_pillars = pd.read_table(fname_pillars, header = None)

post_wgd = ['Vpol', 'Tpha', 'Tbla', 'Ndai', 'Ncas', 'Knag', 'Kafr', 'Cgla', 'Suva', 'Skud', 'Smik','Scer']
post_wgd_p1 = [spec + '_p1' for spec in post_wgd]
post_wgd_p2 = [spec + '_p2' for spec in post_wgd]
post_wgd_p2.reverse()
pre_wgd = ['Zrou', 'Tdel','Klac','Egos','Ecym','Lklu','Lthe','Lwal']
columns = post_wgd_p1 + ['Anc'] + pre_wgd + post_wgd_p2

YGOB_pillars.columns = columns


#    0    V. polyspora Position 1
#    1    T. phaffii Position 1
#    2    T. blattae Position 1
#    3    N. dairenensis Position 1
#    4    N. castellii Position 1
#    5    K. naganishii Position 1
#    6    K. africana Position 1
#    7    C. glabrata Position 1
#    8    S. bayanus var. uvarum Position 1
#    9    S. kudriavzevii Position 1
#    10   S. mikatae Position 1
#    11   S. cerevisiae Position 1
#    12   Ancestral Gene Order
#    13   <NAME>
#    14   <NAME>
#    15   <NAME>
#    16   <NAME> 
#    17   E. cymbalariae
#    18   <NAME>
#    19   L. thermotolerans
#    20   L. waltii
#    21   S. cerevisiae Position 2
#    22   S. mikatae Position 2
#    23   S. kudriavzevii Position 2
#    24   S. bayanus var. uvarum Position 2
#    25   C. glabrata Position 2
#    26   K. africana Position 2
#    27   K. naganishii Position 2
#    28   N. castellii Position 2
#    29   N. dairenensis Position 2
#    30   T. blattae Position 2
#    31   T. phaffii Position 2
#    32   V. polyspora Position 2
gois_gr1_m24 = pd.read_csv(data_processing_dir + '20170817_klscpka/20180327_msn24_gis1_rph1_overlap.csv')
genes = io_library.SC_orf_lookup_by_name(gois_gr1_m24['SC_common_name'])
#for a given list of genes:
#Searches for exact matches of defined motifs in various sequences
motif_dict = {'STRE':'AGGGG',
             'Msn2_pref': 'CGGGG', 
              #'PDS1': 'TTAGGGAT',
              #'PDS2': 'TAAGGGAT',
              'PDS core': 'AGGGA'
             }
output_format = 'full'
sequence_context = 2

output_dict = {}
gene_names = io_library.SC_common_name_lookup(genes) #list(gois_gr1_m24['SC_common_name'])
for ind, gene in enumerate(genes): 
    gene_name = gene_names[ind]
    print(gene_name)
    #Lookup orthologs based on the YGOB Pillars dataframe
    gene_pos = ''
    for pos in ['p1','p2']: 
        Scer_pos = YGOB_pillars['Scer_' + pos]
        Scer_pos_find = Scer_pos[Scer_pos==gene]
        if len(Scer_pos_find)==1:
            if gene_pos != '':
                print('Warning: found gene in two pillars')
            gene_pos = pos
            gene_index = Scer_pos_find.index[0]
        elif len(Scer_pos_find)>1: 
            print('Warning: found gene twice in same pillar')

    #extract intergenic regions from trevor's file

    specs = ['Kafr','Cgla','Suva', 'Skud', 'Smik','Scer']
    ##to include pre wgd species would need to update
    specs_pillar = [spec + '_' + gene_pos for spec in specs]

    orthologs = YGOB_pillars.loc[gene_index, specs_pillar]


    seqs = {}

    for spec in specs: 
        fname_intergenics = os.path.normpath(data_processing_dir + 'promoter_phylogenies/intergenics/' + spec + '_intergenics.faa' )
        record_iterator = SeqIO.parse(fname_intergenics, "fasta")

        for seq in record_iterator: 
            if seq.id == orthologs[spec + '_' + gene_pos]:
                seqs[spec] = seq

    #initial dataframe has species for index, but reindex on seq id before outputting
    seq_ids = [seq.id for seq in seqs.values()]
    output_data_frame = pd.DataFrame(index = seq_ids)
    output_data_frame['spec'] = seqs.keys()

    for motif_name, motif in motif_dict.items():
        L_motif = len(motif)
        output_motif = []

        for spec, seq in seqs.items(): 

            prom_seq = seq.seq
            if len(prom_seq)<no_prom_thresh: 
                output_motif.append('no_promoter')
            else: 
                output_motif_seq = io_library.exact_promoter_scan(motif, prom_seq, output_format, sequence_context)
                output_motif.append(output_motif_seq)
        output_data_frame[motif_name] = output_motif

    output_dict[gene_name] = output_data_frame

#Load my data and Msn2/4 delete data

full_output = pd.DataFrame()
category_dict = dict(zip(gois_gr1_m24['SC_common_name'],gois_gr1_m24['category']))
kl_sc_PKA_data_solis = pd.read_csv(data_processing_dir + os.path.normpath('20170817_klscpka/kl_sc_PKA_data_solis.csv'))
kl_genename_dict = dict(zip(kl_sc_PKA_data_solis['SC_common_name'],kl_sc_PKA_data_solis['kl_genename']))

for gene_name in gois_gr1_m24['SC_common_name']: 
       
    new_output = output_dict[gene_name]
    new_output['SC_common_name'] = gene_name
    new_output['category'] = category_dict[gene_name]
    new_output['kl_genename'] = kl_genename_dict[gene_name]

    full_output = pd.concat([full_output, new_output])

full_output.to_csv(data_processing_dir + 'promoter_phylogenies/gis1_rph1_STRE_PDS_hits.csv')
full_output = pd.read_csv(data_processing_dir + 'promoter_phylogenies/gis1_rph1_STRE_PDS_hits.csv')
new_output
full_output = pd.DataFrame()
pd.concat([full_output, new_output])
motif_dict = {'STRE':'AGGGG',
             'Msn2_pref': 'CGGGG', 
              'PDS1': 'TTAGGGAT',
              'PDS2': 'TAAGGGAT',
              'PDS core': 'AGGGA'
             }
output_format = 'full'
sequence_context = 2
no_prom_thresh = 20 #if length of intergenic region is less than this, call no promoter

gene_name = 'SOL4'
gene = 'YGR248W'

#Lookup orthologs based on the YGOB Pillars dataframe
gene_pos = ''
for pos in ['p1','p2']: 
    Scer_pos = YGOB_pillars['Scer_' + pos]
    Scer_pos_find = Scer_pos[Scer_pos==gene]
    if len(Scer_pos_find)==1:
        if gene_pos != '':
            print('Warning: found gene in two pillars')
        gene_pos = pos
        gene_index = Scer_pos_find.index[0]
    elif len(Scer_pos_find)>1: 
        print('Warning: found gene twice in same pillar')

#extract intergenic regions from trevor's file

specs = ['Kafr','Cgla','Suva', 'Skud', 'Smik','Scer']
##to include pre wgd species would need to update
specs_pillar = [spec + '_' + gene_pos for spec in specs]

orthologs = YGOB_pillars.loc[gene_index, specs_pillar]


seqs = {}

for spec in specs: 
    fname_intergenics = os.path.normpath(data_processing_dir + 'promoter_phylogenies/intergenics/' + spec + '_intergenics.faa' )
    record_iterator = SeqIO.parse(fname_intergenics, "fasta")

    for seq in record_iterator: 
        if seq.id == orthologs[spec + '_' + gene_pos]:
            seqs[spec] = seq




#promoter_counts = io_library.exact_promoter_scan(gene_list, motif_dict, promoter_database, output_format, sequence_context)

#initial dataframe has species for index, but reindex on seq id before outputting
seq_ids = [seq.id for seq in seqs.values()]
output_data_frame = pd.DataFrame(index = seq_ids)
output_data_frame['spec'] = seqs.keys()

for motif_name, motif in motif_dict.items():
    L_motif = len(motif)
    output_motif = []

    for spec, seq in seqs.items(): 

        prom_seq = seq.seq
        if len(prom_seq)<no_prom_thresh: 
            output_motif.append('no_promoter')
        else: 
            output_motif_seq = io_library.exact_promoter_scan(motif, prom_seq, output_format, sequence_context)
            output_motif.append(output_motif_seq)
    output_data_frame[motif_name] = output_motif

output_data_frame
gene_name = 'SOL3'
gene = 'YHR163W'

## Flag if no gene comes up
#Lookup orthologs based on the YGOB Pillars dataframe
gene_pos = ''
for pos in ['p1','p2']: 
    Scer_pos = YGOB_pillars['Scer_' + pos]
    Scer_pos_find = Scer_pos[Scer_pos==gene]
    if len(Scer_pos_find)==1:
        if gene_pos != '':
            print('Warning: found gene in two pillars')
        gene_pos = pos
        gene_index = Scer_pos_find.index[0]
    elif len(Scer_pos_find)>1: 
        print('Warning: found gene twice in same pillar')

#extract intergenic regions from trevor's file

specs = ['Kafr','Cgla','Suva', 'Skud', 'Smik','Scer']
##to include pre wgd species would need to update
specs_pillar = [spec + '_' + gene_pos for spec in specs]

orthologs = YGOB_pillars.loc[gene_index, specs_pillar]


seqs = {}

for spec in specs: 
    fname_intergenics = os.path.normpath(data_processing_dir + 'promoter_phylogenies/intergenics/' + spec + '_intergenics.faa' )
    record_iterator = SeqIO.parse(fname_intergenics, "fasta")

    for seq in record_iterator: 
        if seq.id == orthologs[spec + '_' + gene_pos]:
            seqs[spec] = seq

motif = 'AGGGG'
sequence_context = 2
output_format = 'full'
prom_seq = seq.seq

#exact_promoter_scan(motif, prom_seq, output_format='count', sequence_context=0):
#this function is for a single sequence - same name used to be for one that searched a list of genes 
# now called exact_promoter_scan_genelist
# 
# scans a sequence (biopython object) for a motif in both forward and backward directions 
# If output format is 'count', then just lists the number of hits. 
# If output format is 'full', gives the location, direction (relative to the gene), and
# sequence context (number of bases on either side of the hit)output_format


L_motif = len(motif)

#fwd search
prom_seq_fwd = str(prom_seq)
L_prom = len(prom_seq_fwd)
motif_sites_fwd = [m.start() for m in re.finditer(motif, prom_seq_fwd)]
#to find overlapping motifs
#[m.start() fr m in re.finditer('(?=' + motif + ')', prom_seq_rev)]

#reverse search
#rev search
prom_seq_rev = str(prom_seq.reverse_complement())
#location is in reference to the fwd seq
motif_sites_rev = [len(prom_seq_rev) - m.start() - L_motif for m in re.finditer(motif, prom_seq_rev)]

if output_format == 'count':  
    n_motifs_gene = len(motif_sites_fwd) + len(motif_sites_rev)
    output = n_motifs_gene
elif output_format == 'full': 
    loc_hitseq_gene = []
    for motif_site in motif_sites_fwd: 
        loc = L_prom-1-motif_site
        hitseq = prom_seq_fwd[motif_site-sequence_context:motif_site+L_motif+sequence_context]  
        hitdir = 'fwd'
        loc_hitseq_gene.append((loc,hitdir,hitseq))
    for motif_site in motif_sites_rev: 
        loc = L_prom - motif_site+L_motif
        hitseq = prom_seq_rev[(L_prom-motif_site-L_motif-sequence_context):(L_prom-motif_site+sequence_context)]
        hitdir = 'rev'
        loc_hitseq_gene.append((loc,hitdir,hitseq))

    if len(loc_hitseq_gene)==0: 
        output = None
    else:
        output = loc_hitseq_gene
else: 
    print("choose output_format: 'count' or 'full' ")
        

#suspect Smik gnd2 is actually Smik gnd1 - check intergenics
#RGI1 overlaps with dubious orf.  
#Smik 4.730 is missing from intergenics.  I think that intergenics is not as good as 700bp upstream  
#YER067W overlaps with dubious orf so isn't in trevor's file.  If you take the promoter as the 1kb upstream then it there are several STREs and one PDS. 


gene_names<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10314808301090397, 'fraction_numerical': 0.016795638410130143, 'mean_word_length': 3.447790379350802, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '660663', 'license': 'MIT', 'n_tokens_mistral': 4382, 'n_tokens_neox': 4079, 'n_words': 1017, 'path': 'scripts/BMH_code_for_trevors_intergenics_files.ipynb', 'repo_name': 'heineike02/y1000plus_tools', 'size': '11396'}