Help to understand the code for dipeptide composition calculation (in python)

Question

Dipeptide composition of a protein sequence is the number of times a particular dipeptide (e.g. Arginine-Histidine) occurs in a sequence divided by the total number of dipeptides in the sequence (which is the length of the sequence - 1)

I have found the following code for this:

import re

def DPC(fastas):
    AA = 'ACDEFGHIKLMNPQRSTVWY'
    encodings = []
    diPeptides = [aa1 + aa2 for aa1 in AA for aa2 in AA]
    header = ['#'] + diPeptides
    encodings.append(header)

    AADict = {}
    for i in range(len(AA)):
        AADict[AA[i]] = i

    for i in fastas:
        name, sequence = i[0], re.sub('-', '', i[1])
        code = [name]
        tmpCode = [0] * 400
        for j in range(len(sequence) - 2 + 1):
            tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] = tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] +1
        if sum(tmpCode) != 0:
            tmpCode = [i/sum(tmpCode) for i in tmpCode]
        code = code + tmpCode
        encodings.append(code)
    return encodings

There is another file (main) from where the arguments fastas come from which is basically a list of strings where every string is a protein sequence. I am finding it really difficult to understand the nested loops that were used:

    for i in fastas:
        name, sequence = i[0], re.sub('-', '', i[1])
        code = [name]
        tmpCode = [0] * 400
        for j in range(len(sequence) - 2 + 1):
            tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] = tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] +1
        if sum(tmpCode) != 0:
            tmpCode = [i/sum(tmpCode) for i in tmpCode]
        code = code + tmpCode
        encodings.append(code)

I have spent 4-5 hours on this but I am still finding it really difficult to follow. I would really appreciate if anyone can explain the steps involved in these loops. Thank you!

Edit: The github link to the files are :

1. DPC.py
2. iFeature.py

The github link to the python toolkit: iFeature

Answer 1

A few additional hints that may help clarify the algorithm:

The 20 magic number is because there are twenty amino acids. The 400 magic number is because there are 20 * 20 potential dipeptides.

These peptides are then hashed by converting their letter to a number in the range [0, 20). You can then multiple one of them by 20 and add the other to get a number between [0, 400). In the case of this function the first amino acid "code" is multiplied by 20. This number can be converted back to the dipeptide sequence when outputting the data. The translation table they use is generated here diPeptides = [aa1 + aa2 for aa1 in AA for aa2 in AA].

The nested loops are to loop through each entry and then each sequence to get every dipeptide sequence, which can be done in linear time.

Just as a quick note on hashing DNA and peptide sequences. You can generally convert them to an array index by using increasing powers of the alphabet size multiplied by the letter index of the alphabet. For example, for DNA you can convert A=0, C=1, G=2, T=3, and you can hash fourmers by using 4^4 * a[i] + 4^3 * a[i + 1] + 4^2 * a[i + 2] + 4^1 * a[i + 3]. In code you may see this performed with bitshifts for DNA since the alphabet size is a power of 2.

The advantage of this hashing scheme is. thateach sequence is guaranteed to have a unique index. The problem with this method is that as if you wish to use longer k-mers the memory usage will grow exponentially. For example, for 16-mers you need 4 billion combinations, which is the size of a 32 bit unsigned integer.