I'm currently trying to write a program that will calculate Mutual Information given text files of nucleotide distributions.

Below is the MI program, as seen I am currently stuck after finding marginal distribution, then dividing the values by the number of bins (8) in my Sort-Seq experiment:

import pandas as pd
import numpy as np
import sys

filename = sys.argv[1]

pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
pd.set_option('display.width', None)
pd.set_option('display.max_colwidth', None)

df = pd.read_csv(filename, header=0, sep= ',', skipinitialspace=True)
headers = ['A', 'T', 'G', 'C']
df.columns = headers

df.head(499)

df['max'] = df[['A', 'T', 'G', 'C']].max(axis=1)
df['sum'] = df[['A', 'T', 'G', 'C']].sum(axis=1)

df.loc[:,"A":"C"] = df.loc[:,"A":"C"].div(df["sum"], axis=0)
df['mutation_rate'] = (1-df['max']/df['sum'])
df['max2'] = df[['A', 'T', 'G', 'C']].max(axis=1)
df['sum2'] = df[['A', 'T',  'G', 'C']].sum(axis=1)
df['marginal_distribution']=(1-df['max2']/df['sum2'])

df.head()

#df.head()
#df['A/8'] = df['A'].div(8)
#df['T/8'] = df['T'].div(8)
#df['G/8'] = df['G'].div(8)
#df['C/8'] = df['G '].div(8)
#df.head() --> **8 bins of sequences** 

sys.stdout = open("stats_" + filename, "w")
print(df)
sys.stdout.close()

With the output

    A   T   G   C
0   0.000671    0.000471    0.00028     0.998578
1   0.000591    0.000319    0.000048    0.999042
2   0.999169    0.000351    0.000192    0.000288
3   0.000024    0.000351    0.000032    0.999593
4   0.999297    0.000184    0.000224    0.000296
5   0.00004     0.000184    0.000032    0.999744
6   0.999497    0.000064    0.000144    0.000295
7   0.999329    0.000256    0.000112    0.000303
8   0.000072    0.0002      0.000064    0.999665


 max     sum      mutation_rate
125032  125210  0.001422
125082  125202  0.000958
125107  125211  0.000831
125161  125212  0.000407
125122  125210  0.000703
125180  125212  0.000256
125149  125212  0.000503
125124  125208  0.000671
125170  125212  0.000335

max2     sum2
0.998578    1
0.999042    1
0.999169    1
0.999593    1
0.999297    1
0.999744    1
0.999497    1
0.999329    1
0.999665    1

marginal_distribution
0.001422
0.000958
0.000831
0.000407
0.000703
0.000256
0.000503
0.000671
0.000335

A/numberOfBins  T/numberOfBins  G/numberOfBins  C/numberOfBins
0.000084    0.000059    0.000035    0.124822
0.000074    0.00004     0.000006    0.12488
0.124896    0.000044    0.000024    0.000036
0.000003    0.000044    0.000004    0.124949
0.124912    0.000023    0.000028    0.000037
0.000005    0.000023    0.000004    0.124968
0.124937    0.000008    0.000018    0.000037
0.124916    0.000032    0.000014    0.000038
0.000009    0.000025    0.000008    0.124958

As of right now, each bin is printed to its own output file as seen above. I will need to get specific values from each nucleotide at each specified position across the 8 bins I am using to properly carry out the calculation.

For instance, in the example output file above, 0.000084 is represented at position 1 for A.

When completing this in Excel I use the following formula. Where H294 = 0.000084. The other values in row 294 are the other bins at this exact location.

=H294*LOG(H294/(0.125*SUM(H294,N294,T294,Z294,AF294,AL294,AR294,AX294)),2)

How should I set up the input file? Can this be done in df / called with one line? What line of code do I need here?



Source link