0
$\begingroup$

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

$\endgroup$
2
  • 2
    $\begingroup$ One thing that may help when trying to understand other peoples code is to add print statements and possibly break statements at various parts of the code. You can print out exactly what is happening in the loop; for example, in the nested loop print j, print the tmpCode index, print the value at the tmp code etc. $\endgroup$
    – GWW
    Sep 11 at 14:50
  • $\begingroup$ @user438383, deleted $\endgroup$
    – Noob
    Sep 12 at 12:06
3
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Thank you so much. I am trying to follow what your suggestion. I am getting all sorts of error messages. I have very little programming skills. I only learnt to work with python on jupyter notebook. This toolkit (or any toolkit, I presume). works across multiple files at a time I presume. I am completely lost. $\endgroup$
    – Noob
    Sep 11 at 15:27
  • $\begingroup$ I have understood the concept of hashing i.e. AADict={'A': 0, 'C': 1, 'D': 2, 'E': 3, 'F': 4, 'G': 5, 'H': 6, 'I': 7, 'K': 8, 'L': 9, 'M': 10, 'N': 11, 'P': 12, 'Q': 13, 'R': 14, 'S': 15, 'T': 16, 'V': 17, 'W': 18, 'Y': 19} $\endgroup$
    – Noob
    Sep 11 at 15:30
  • $\begingroup$ I understood multiplying one of them by 20 and add the other to get a number between [0, 400). But still these 2 line if sum(tmpCode) != 0: tmpCode = [i/sum(tmpCode) for i in tmpCode] are really bothering me $\endgroup$
    – Noob
    Sep 11 at 15:49
  • 1
    $\begingroup$ They are checking if the sum is larger than zero. Then they are normalizing the totals to sum to 1.0; ie. calculating the frequency of each dipeptide. $\endgroup$
    – GWW
    Sep 11 at 21:25
  • $\begingroup$ It's working now. Thank you so much. I really appreciate $\endgroup$
    – Noob
    Sep 12 at 12:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.