This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to handle the phylip sequence file through amino acid characters?

I am using python 3 and biopython 1.72 and I have a protein sequence file in phylip format. This is the example file:

     6    200
Human      ---------- ---------- ---------- ---------- ---------- 
Chimpanzee ---------- ---------- ---------- ---------- ---------- 
Dog        ---------- ---------- ---------- ---------- ---------- 
Mouse      ---------- ---------- ---------- ---------- ---------- 
Xenopus    ---------- ---------- ---------- ---------- ---------- 
Amphioxus  MQWTGFRVSM TTLMMIMGVV AVLIALLPAK AQQPHDKSLR TTSTLTDTGA 

           ---------- ---------- ---------- ---------- ---------- 
           ---------- ---------- ---------- ---------- ---------- 
           ---------- ---------- ---------- ---------- ---------- 
           ---------- ---------- ---------- ---------- ---------- 
           ---------- ---------- ---------- ---------- ---------- 
           SADEADMGSA HVELLDGDDD VGNGSDQMMV TLHLQSIFQC IRRPCEKVDR 

           ---------- ---------- ---------- --------MR LRVRLLKRTW 
           ---------- ---------- ---------- --------MR LRVRLLKRTW 
           ---------- ---------- ---------- --------MK LRVRLQKRTW 
           ---------- ---------- ---------- --------MK LRVRLQKRTQ 
           ---------- ---------- ---------- --------MK LRVRVRKQTN 
           AIDPVTQRWR TANTRNDYQK INVCVVPAYD VSLSTGVRMK LRVKISGQKT 

           PLEVPETEPT LGHLRSHLRQ SLLCTWGYSS NTRFTITLNY KDPLTGDEET 
           PLEVPETEPT LGHLRSRLRQ SLLCTWGYSS NTRFTITLNY KDPLTGDEET 
           PLDLPDAEPT LGQLRAHLSQ ALLPSWGFGS DTRFAITLNN KDALTGDEET 
           PLEVPESEPT LGQLRAHLSQ VLLPTLGFSS DTRFAITLNN KDALTGDEET 
           RLELEAESPT LGDLRSKLSS VTLPALGYST EANFTITLNG KDALTGDQNT 
           RVDVGQDCHT LGTLRTLLAP VLGEQYGLGD DMPFEISLNG RDALLGDDKP

I want to navigate through the file in a way that ;

  1. If the columns of the file have same amino acid, it saves and prints the amino acid that exists in the column and puts a '#' in the file at the end of that column.
  2. I need to know how can I search the columns by using the amino acids, if particular set of amino acids exist in that column or not. For example, if 'STA'or 'HY' or 'FVLIM' or 'NDEQHK' exists in any of the columns, then put a '@' in the file at the end of that column.

Following is the code I have been trying to manipulate the file with:

alignment = AlignIO.read(open("example.phy"), "phylip")
v1 = v2 = v3 = 0
for col in zip(*alignment):
    num_unique = len(set(col))
    if num_unique == 1 and col[0] != '-':
        print (num_unique)
        v1 += 1

    elif num_unique > 1 and '-' not in col:
        print (num_unique)
        v2 += 1
    elif '-' in col:  # assumes 1 or more dashes
        v3 += 1

print('Number of columns with the same amino acid: {}\n'
      'Number of columns with at least 2 amino acids (no gaps): {}\n'
      'Number of columns with one/more gaps: {}'
       .format(v1, v2, v3))

These variables return the occurrences of amino acids in the columns of the file but I don't understand how to search using the amino acids. Moreover how can I put a character of my choice in a desired column of the file ???

python biopython

Rather than reinventing the wheel, why not just use BioPython's built in consensus-deriving tools?

See for example: http://biopython.org/DIST/docs/api/Bio.Align.AlignInfo.SummaryInfo-class.html

You can easily get a so-called 'dumb-consensus', and just set your threshold to 1.

The syntax for accessing a column within the alignment is:

aln[:, 1]   # if 'aln' was the handle for AlignIO.read()

So you can loop through all columns as:

[aln[:,col] for col in xrange(len(aln[0]))]

Alright, I am gonna try using this to access my point of interest. Thankyou for the response.

I'm extremely sorry for this mistake. Its just that I was a little stressed about my work and I have done this by mistake. I'll remove this at once. This is so disappointing I know, but I'm sorry, I could'nt get how I have done this.

Hey, it's OK. Don't worry about it.

Hello mdsiddra!

We believe that this post does not fit the main topic of this site.

This is being closed because this has been asked earlier somewhere else.

For this reason we have closed your question. This allows us to keep the site focused on the topics that the community can help with.

If you disagree please tell us why in a reply below, we'll be happy to talk about it.

Cheers!

0 answers

No answers yet.

Log in to answer this question.