Python for Bioinformatics: 16S rRNA V regions, continued

Tuesday, March 8, 2011

16S rRNA V regions, continued

As promised last time (here), this is the code to plot V regions from the sequences in the Enterobacteriales, obtained from RDP. If you compare to this post from the other day, you'll see we match pretty well.

The figure looks better with a wider window, but the limits are established more accurately with a smaller window. The threshold T needs to be adjusted based on the window size.

One detail: I did the redirect below to save the extreme values and then parsed them with the code at the end of the post.

python script.py > extreme_values.txt

script.py

import utils as ut
from info import shannon
import matplotlib.pyplot as plt

data = ut.load_data('entero.txt')
data = data.strip().split('>')[1:]
EC = [e for e in data if 'X80725' in e][0]
EC = ut.clean_fasta(EC)[1]
data = [ut.clean_fasta(e)[1] for e in data]
L = ut.make_count_list(data)

pos = 0
R = range(1,1451)
iL = list()

for i,c in enumerate(EC):
    if c == '-':  continue
    pos += 1
    if not pos in R:  continue
    cD = L[i]
    e = shannon(cD,'ACGT')
    iL.append(e)
    #print str(pos).ljust(3), c + '  ', 
    #print ''.join([str(cD[k]).ljust(5) for k in 'ACGT']),
    #print round(e,2)

aL = list()
w = 20
T = 1.8
for i in range(len(iL)):
    j, k = i - w, i + w + 1
    if j < 0:  j = 0
    if k > len(iL):  k = len(iL)
    m = ut.mean(iL[j:k])
    if m < T:  print i+1
    aL.append(m)
    
plt.plot((1,1451),(T,T),lw=2,color='r',zorder=0)
plt.scatter(R,aL)
ax = plt.axes()
ax.set_xlim(-5,1455)
ax.set_ylim(0.8,2.05)
plt.savefig('example.pdf')

analyze.py

import utils as ut
data = ut.load_data('extreme_values.txt')
L = [int(n) for n in data.strip().split('\n')]

current = L.pop(0)
print current, '-',
while L:
    next = L.pop(0)
    if next != current + 1:
        print current
        print next, '-',
    current = next
print next