Bring in Alignment for mapping

This program will map TFBS using the Biopython's motif package.

Inputs:

  1. before alignment (fasta)
  2. after alignment (fasta)
  3. TFBS Position Frequency Matrix.

In [8]:
from Bio import motifs
from Bio import SeqIO 
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC, generic_dna, generic_protein
import re
import pandas as pd
import numpy as np
import os, sys

In [9]:
## Input 1 - Alignment Input 

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

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

## Turn sequences into a list of strings
## They are no longer bio.seq.seq objects though
alignment_string_list = []
for seq in alignment:
    alignment_string_list.append(str(seq.seq))


Found 9 records in alignment file

In [21]:
## Input 2 - 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()))


Found 9 records in raw sequence file
<class 'Bio.SeqRecord.SeqRecord'>
<class 'Bio.SeqRecord.SeqRecord'>
<class 'Bio.SeqRecord.SeqRecord'>
<class 'Bio.SeqRecord.SeqRecord'>
<class 'Bio.SeqRecord.SeqRecord'>
<class 'Bio.SeqRecord.SeqRecord'>
<class 'Bio.SeqRecord.SeqRecord'>
<class 'Bio.SeqRecord.SeqRecord'>
<class 'Bio.SeqRecord.SeqRecord'>

In [11]:
## Input 3 - Motif Input
## [ ] This is where I need to loop through all PWMs in a file

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()


        0      1      2      3      4      5      6      7
A:   0.19   0.17   0.88   0.92   0.04   0.04   0.06   0.12
C:   0.37   0.08   0.04   0.02   0.02   0.87   0.52   0.25
G:   0.08   0.04   0.04   0.04   0.33   0.02   0.08   0.37
T:   0.37   0.71   0.04   0.02   0.62   0.08   0.35   0.27


In [19]:
## Searching the Motifs in Sequences
## Returns a list of arrays with a score for each position

pssm_list = [ ]
for seq in raw_sequences_2:
    pssm_list.append(pssm.calculate(seq))

## Approximate calculation of appropriate thresholds for motif finding 
## Patser Threshold
## It selects such a threshold that the log(fpr)=-ic(M) 
## note: the actual patser software uses natural logarithms instead of log_2, so the numbers 
## are not directly comparable. 

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

print("Patser Threshold is %5.3f" % patser_threshold) # Calculates Paster threshold.


Patser Threshold 3.262

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

## This just gets all the scores for the chance of bcd in each position of the raw sequences
score_list = []
position_list = []

for position, score in pssm.search(raw_sequences_2[0], threshold= patser_threshold):
    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)
    
## Check   
## print(position_list_pos)

strand = []

for x in position_list:
    if x < 0:
       strand.append("negative")
    else:
       strand.append("positive")

In [14]:
## 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)
print(alignment_id[1])
print(type(alignment_id[1]))


[224, 221, 403, 596, 712, 881, 869, 903, 1114]
ludwig_eve-striped-2||MEMB002A|+
<type 'str'>

In [15]:
## Make dataframe that has everything
pos_df = pd.DataFrame(
    {'raw_position': position_list,
     'raw_position_pos_only': position_list_pos,
     'alignment_position':position_list_pos,
     'strand_direction': strand,
     'score': score_list,
     'species': alignment_id[1]
    })

#############
## [ ] this needs to be a value for species!!
##     like pos_df = alignment_id[1]
############

## pos_df['species'] = alignment_id[1]
print(pos_df)


   alignment_position  raw_position  raw_position_pos_only      score  \
0                 157           157                    157  10.457056   
1                 154          -751                    154   8.946094   
2                 307          -598                    307   6.417715   
3                 455          -450                    455   9.909568   
4                 567          -338                    567   9.909568   
5                 675           675                    675  10.016483   
6                 663          -242                    663   8.959556   
7                 690           690                    690   6.349059   
8                 846           -59                    846   6.380240   

                            species strand_direction  
0  ludwig_eve-striped-2||MEMB002A|+         positive  
1  ludwig_eve-striped-2||MEMB002A|+         negative  
2  ludwig_eve-striped-2||MEMB002A|+         negative  
3  ludwig_eve-striped-2||MEMB002A|+         negative  
4  ludwig_eve-striped-2||MEMB002A|+         negative  
5  ludwig_eve-striped-2||MEMB002A|+         positive  
6  ludwig_eve-striped-2||MEMB002A|+         negative  
7  ludwig_eve-striped-2||MEMB002A|+         positive  
8  ludwig_eve-striped-2||MEMB002A|+         negative  

In [ ]: