1
$\begingroup$

I have 2 lists, list1 and list2, of protein sequences of the same given gene in different strains. In list1 are reported amino-acid sequences from resistant strains, and in list2 are reported amino-acid sequences from susceptible strains. Each item of each list is an amino-acid sequence.

I want to filter out (in Python 3) from list1 all of the variants that are also present in list2. The final aim is to have a list of the sequences that have unique mutations that are not in the 2nd list. The mutations can be of any type (e.g., indels, mismatches).

All the sequences are of at least 90% identity to some reference sequence (filtered out the low matches using seq alignment).

Any enlightenment / package / tool for this job will be appreciated!

Toy example:

a, b, c : amino-acids

For this input:

  1. Reference seq: aaabbb
  2. Resistant strains: list1 = [aacbbc, ababbb]
  3. Susceptible strains: list2 = [aacbbb, aaabbc]

The expected output would be: [ababbb], since the mutations in list2 cancels out the item aacbbc.

Of course this is a simple example and it is possible to have more complex scenarios (e.g., longer indels and mismatches).

|improve this question
$\endgroup$
  • $\begingroup$ Is there only one mutation per sequence entry in the lists? Could you align the two lists to each other and look for perfect matches? $\endgroup$ – Bioathlete Nov 10 '17 at 15:23
  • $\begingroup$ For the first step, I indeed filtered out the exact matches. But I also need to filter out other cases. For example if I have a sequence in list1 with 2 mutations and I have 2 sequences in list2 that each one of them has one of those mutations, then I want to filter this sequence from list 1 as well. $\endgroup$ – Kobi T Nov 10 '17 at 15:26
  • $\begingroup$ In which format is each list? It seems in a Fasta format but I am unsure because usually the mutations are not stored in these format. $\endgroup$ – llrs Nov 10 '17 at 15:41
  • $\begingroup$ The lists contain strings (extracted from transcriptome fasta files) of the amino acids that got score of at least 90% of the max score against some reference sequence. $\endgroup$ – Kobi T Nov 10 '17 at 15:43
  • 1
    $\begingroup$ You will have to look at the alignment tool you are uses and see if it supports it directly or if it will write out a BAM file then you can use samtools to generate the VCF. $\endgroup$ – Bioathlete Nov 10 '17 at 19:00
1
$\begingroup$

Samtools/VCF is overkill for comparing amino acid sequences from a single protein. I would suggest making pairwise alignments between each sequence in list2 and the reference and making a list of mutations (2c and 6c in your example). Then repeat the procedure for list1, discarding sequences if they contain mutations in your stored list. This can all be done using biopython. Here's how to make the alignments to get you started:

from Bio import pairwise2
from Bio import SeqIO
seq1 = SeqIO.read("alpha.faa", "fasta")
seq2 = SeqIO.read("beta.faa", "fasta")
alignments = pairwise2.align.globalxx(seq1.seq, seq2.seq)
|improve this answer
$\endgroup$
  • $\begingroup$ How can I get a list of differences between 2 proteins in a way that is convenient to work with? Is there a tool for that? It's not as simple as comparing strings (even after alignment). There are many possibilities (gap in one/gap in other/mismatch of different types/...) $\endgroup$ – Kobi T Nov 14 '17 at 15:16
  • $\begingroup$ Actually I think insertion deletion and mismatch are the only possibilities. Just store the position of the mutation in the reference sequence (or positions in the case of a deletion) and the sequence of the reference at that position/positions. You could also make a multiple sequence alignment of all sequences and discard sequences that from list1 that share non-reference characters with sequences in list2 without needing to keep track of reference positions, but that would require an external alignment program $\endgroup$ – heathobrien Nov 14 '17 at 15:26

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.