<jupyter_start><jupyter_text>Bring in Alignment for mappingThis program will map TFBS using the Biopython's motif package.**Inputs**: 1. before alignment (fasta) 2. after alignment (fasta) 3. TFBS Position Frequency Matrix.<jupyter_code>from Bio import motifs
from Bio import SeqIO 
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC, generic_dna, generic_protein
import pandas as pd
import numpy as np
import os, sys
## Alignment Input

# read in alignment as a list of sequences
alignment = list(SeqIO.parse("../data/fasta/output_ludwig_eve-striped-2.fa", "fasta"))

# Check
print("Found %i records in alignment file" % len(alignment))

## Turn a sequences into a list of strings
alignment_string_list = []
for seq in alignment:
    alignment_string_list.append(str(seq.seq))

alignment_id = []
for seq in alignment:
    alignment_id.append(str(seq.id))

## Do I even need turn int a dataframe?
alignment_df = pd.DataFrame(
    {'id': alignment_id,
     'align_seq': alignment_string_list
    })


## Check
print list(alignment_df)
print type(alignment_df)

##############################
# [x] I feel weird going ahead with out seq ids. Maybe at this point turn into a dictionary?
# [x] No, what I really need is to turn this into a dataframe
# [ ] Now I need to make sure I can use the seqences to continue on. Since I turned it into a dataframe.
## Raw Sequences Input

raw_sequences = list(SeqIO.parse("../data/fasta/ludwig_eve-striped-2.fasta", "fasta"))
print("Found %i records in raw sequence file" % len(raw_sequences))

# make all IUPAC.IUPACUnambiguousDNA()
raw_sequences_2 = []

for seq in raw_sequences:
    raw_sequences_2.append(Seq(str(seq.seq), IUPAC.IUPACUnambiguousDNA()))

# Check

#print raw_sequences_2
#print type(raw_sequences_2)

# Check
#for seq in raw_sequences_2:
    #print(seq.alphabet)
    #print(type(seq))
## Motif Input

bcd = motifs.read(open("../data/PWM/transpose_fm/bcd_FlyReg.fm"),"pfm")
print(bcd.counts)
pwm = bcd.counts.normalize(pseudocounts=0.0)
pssm = pwm.log_odds()
print(pssm.alphabet)
print(type(raw_sequences_2))
## Searching the Sequences
pssm_list = [ ]
for seq in raw_sequences_2:
    pssm_list.append(pssm.calculate(seq))

## Check
#print(pssm_list)
#for seq in pssm_list:
    #print("Background: %f" % bcd.pssm.mean(bcd.background))
    
################################
# [ ] Its the same background for all the sequences? That weird.  Right?
################################

# Patser Threshold

distribution = pssm.distribution(background=bcd.background, precision=10**4)
threshold = distribution.threshold_patser()

print("Patser Threshold %5.3f" % threshold) #automatically calulate Paster threshold. 
print("nothing")
position_list = []
score_list = []

###################################
# [ ] Need to reiterate over raw_sequences_2
# [ ] When reiterating over raw_sequences_2, attach id
##################################

for position, score in pssm.search(raw_sequences_2[0], threshold=6):
    position_list.append(position)
    score_list.append(score)

# Change position to positive
position_list_pos = []
for x in position_list:
    if x < 0:
       position_list_pos.append(905 + x)
    else:
       position_list_pos.append(x)
#print(position_list_pos)

strand = []
for x in position_list:
    if x < 0:
       strand.append("negative")
    else:
       strand.append("positive")
## get alignment position using `alignment_string_list`

remap_dict = {}
nuc_list = ['A', 'a', 'G', 'g', 'C', 'c', 'T', 't', 'N', 'n']
counter = 0

#######################
# [ ] Reiterate through all species?
# [ ] maybe create a list of dictionaries?
#######################


for xInd, x in enumerate(alignment_string_list[1]):    
    if x in nuc_list:
        remap_dict[counter] = xInd
        counter += 1
# Check
# print(remap_dict)

# Now find the value from the key??? Find the alignment posititon from raw position

align_pos = [remap_dict[x] for x in position_list_pos]

# check
print(align_pos)
# Make dataframe that has everything
pos_df = pd.DataFrame(
    {'raw_positition': position_list,
     'raw_position_pos_only': position_list_pos,
     'alignment_pos':position_list_pos,
     'strand_direction': strand,
     'score': score_list
    })

print(pos_df)<jupyter_output>raw_position_pos_only  raw_positition      score strand_direction
0                    145             145  10.457056         positive
1                    165            -740   8.946094         negative
2                    457            -448   9.909568         negative
3                    548            -357   9.909568         negative
4                    641             641  10.016483         positive
5                    823             -82   6.380240         negative
{'fraction_non_alphanumeric': 0.14312936692780012, 'fraction_numerical': 0.02477239042981156, 'mean_word_length': 3.604288499025341, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '512457', 'license': 'MIT', 'n_tokens_mistral': 1609, 'n_tokens_neox': 1494, 'n_words': 424, 'path': 'py/.ipynb_checkpoints/TFBSmapping_draft-checkpoint.ipynb', 'repo_name': 'iamciera/montium_analyses_2_conservation', 'size': '4723'}