0

I need to represent many gene sequences using chaos game representation I got this python code from Boštjan Cigan's blog (https://bostjan-cigan.com/chaos-game-representation-of-gene-structure-in-python/)

Author: Bostjan Cigan

https://bostjan-cigan.com

    import collections
    from collections import OrderedDict
    from matplotlib import pyplot as plt
    from matplotlib import cm

    import pylab
    import math

    f = open("ensemblSeq.fa")
    s1 = f.read()
    data = "".join(s1.split("\n")[1:])

    def count_kmers(sequence, k):
        d = collections.defaultdict(int)
        for i in xrange(len(data)-(k-1)):
            d[sequence[i:i+k]] +=1
        for key in d.keys():
             if "N" in key:
                 del d[key]
        return d

    def probabilities(kmer_count, k):
        probabilities = collections.defaultdict(float)
        N = len(data)
        for key, value in kmer_count.items():
            probabilities[key] = float(value) / (N - k + 1)
        return probabilities

    def chaos_game_representation(probabilities, k):
        array_size = int(math.sqrt(4**k))
        chaos = []
        for i in range(array_size):
            chaos.append([0]*array_size)

        maxx = array_size
        maxy = array_size
        posx = 1
        posy = 1
         for key, value in probabilities.items():
             for char in key:
                 if char == "T":
                    posx += maxx / 2
                 elif char == "C":
                    posy += maxy / 2
                 elif char == "G":
                     posx += maxx / 2
                     posy += maxy / 2
                 maxx = maxx / 2
                 maxy /= 2
             chaos[posy-1][posx-1] = value
             maxx = array_size
             maxy = array_size
             posx = 1
             posy = 1

         return chaos

      f3 = count_kmers(data, 3)
      f4 = count_kmers(data, 4)

      f3_prob = probabilities(f3, 3)
      f4_prob = probabilities(f4, 4)

      chaos_k3 = chaos_game_representation(f3_prob, 3)
      pylab.title('Chaos game representation for 3-mers')
      pylab.imshow(chaos_k3, interpolation='nearest', cmap=cm.gray_r)
      pylab.show()

      chaos_k4 = chaos_game_representation(f4_prob, 4)
      pylab.title('Chaos game representation for 4-mers')
      pylab.imshow(chaos_k4, interpolation='nearest', cmap=cm.gray_r)
      pylab.show()

This code works fine but I have many sequence files I need to iterate through each fasta file in the folder and get individual plots stored in a folder with the name of the image file corresponding to the name of the fasta file how can I modify the code according to my need

I am new to python as well as StackOverflow if any mistake is there kindly ignore

Thanks in advance

Lakshmi KrishnaKumaar
  • 39
  • 1
  • 1
  • 9

1 Answers1

0

So, if you want to apply your code on every file in your directory a very simple way to do this is calling all files inside of a for-loop. I suggest the following:

import collections
import os
from collections import OrderedDict
from matplotlib import pyplot as plt
from matplotlib import cm

import pylab
import math

def count_kmers(sequence, k):
    d = collections.defaultdict(int)
    for i in xrange(len(data)-(k-1)):
        d[sequence[i:i+k]] +=1
    for key in d.keys():
         if "N" in key:
             del d[key]
    return d

def probabilities(kmer_count, k):
    probabilities = collections.defaultdict(float)
    N = len(data)
    for key, value in kmer_count.items():
        probabilities[key] = float(value) / (N - k + 1)
    return probabilities

def chaos_game_representation(probabilities, k):
    array_size = int(math.sqrt(4**k))
    chaos = []
    for i in range(array_size):
        chaos.append([0]*array_size)

    maxx = array_size
    maxy = array_size
    posx = 1
    posy = 1
    for key, value in probabilities.items():
        for char in key:
            if char == "T":
                posx += maxx / 2
            elif char == "C":
                posy += maxy / 2
            elif char == "G":
                 posx += maxx / 2
                 posy += maxy / 2
            maxx = maxx / 2
            maxy /= 2
        chaos[posy-1][posx-1] = value
        maxx = array_size
        maxy = array_size
        posx = 1
        posy = 1
    return chaos

if __name__ == "__main__":
    PATH = os.getcwd()
    filelist = sorted([os.path.join(PATH, f) for f in os.listdir(PATH) if f.endswith('.fa')])
    for file in filelist:
        f = open(file)
        s1 = f.read()
        data = "".join(s1.split("\n")[1:])
        f3 = count_kmers(data, 3)
        f4 = count_kmers(data, 4)

        f3_prob = probabilities(f3, 3)
        f4_prob = probabilities(f4, 4)

        chaos_k3 = chaos_game_representation(f3_prob, 3)
        pylab.title('Chaos game representation for 3-mers')
        pylab.imshow(chaos_k3, interpolation='nearest', cmap=cm.gray_r)
        pylab.savefig(os.path.splitext(file)[0]+'chaos3.png')
        pylab.show()

        chaos_k4 = chaos_game_representation(f4_prob, 4)
        pylab.title('Chaos game representation for 4-mers')
        pylab.imshow(chaos_k4, interpolation='nearest', cmap=cm.gray_r)
        pylab.savefig(os.path.splitext(file)[0]+'chaos4.png')
        pylab.show()

I just wrapped a loop around and added a pylab.savefig() call. Furthermore I used os to get the filenames from your directory. It should work now.

Franz
  • 623
  • 8
  • 14
  • I am getting an unindent error near maxx = maxx /2 – Lakshmi KrishnaKumaar Jun 13 '17 at 10:08
  • The indentation was messed up a little. I fixed it. By the way, I did not change your functions. So you could also use the original indentation. – Franz Jun 13 '17 at 10:20
  • can you explain what this actually does if __name__ == "__main__" – Lakshmi KrishnaKumaar Jun 13 '17 at 11:01
  • This is just a check, whether you are executing your script directly. If you just import your `your_file.py` with `import your_file` to another script its just importing the classes and functions but not the code inside of `main`. It makes it more simple to reuse your code. I would recommend you to read https://stackoverflow.com/questions/419163/what-does-if-name-main-do – Franz Jun 13 '17 at 11:06
  • ' chaos_k3 = chaos_game_representation(f3_prob, 3) matrix1 = np.array(chaos_k3).reshape(8,8) print matrix1' I used this code in the above program to get the output as array but I get the matrix as [[ 0.03061224 0.01391466 0.01113173 0.02319109 0.01669759 0.01020408 0.01669759 0.02690167] (I am getting 8 similar lines pasting one for sample due to space issue – Lakshmi KrishnaKumaar Jun 16 '17 at 05:03