<p>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. <a href="http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc86" rel="nofollow noreferrer">Here's</a> how to make the alignments to get you started:</p>

<pre><code>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)
</code></pre>

CC BY-SA 3.0
2017-11-14T15:16:23.367Z
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/...)
Kobi T
CC BY-SA 3.0
2017-11-14T15:26:42.393Z
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
heathobrien
CC BY-SA 3.0
2017-11-14T15:06:42.967Z
2017-11-14T15:06:42.967Z
heathobrien
<p>I have 2 lists, <code>list1</code> and <code>list2</code>, of protein sequences of the same given gene in different strains. In <code>list1</code> are reported amino-acid sequences from resistant strains, and in <code>list2</code> are reported amino-acid sequences from susceptible strains. Each item of each list is an amino-acid sequence. </p>

<p>I want to filter out (in Python 3) from <code>list1</code> all of the variants that are also present in <code>list2</code>. 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).</p>

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

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

<p>Toy example:</p>

<p><code>a</code>, <code>b</code>, <code>c</code> : amino-acids</p>

<p>For this input: </p>

<ol>
<li>Reference seq: <code>aaabbb</code></li>
<li>Resistant strains: <code>list1</code> = [aa<strong>c</strong>bb<strong>c</strong>, a<strong>b</strong>abbb]</li>
<li>Susceptible strains: <code>list2</code> = [aa<strong>c</strong>bbb, aaabb<strong>c</strong>]</li>
</ol>

<p>The expected output would be: [a<strong>b</strong>abbb], since the mutations in <code>list2</code> cancels out the item aa<strong>c</strong>bb<strong>c</strong>.</p>

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

CC BY-SA 3.0
2017-11-10T15:23:15.980Z
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?
Bioathlete
CC BY-SA 3.0
2017-11-10T15:26:22.907Z
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.
Kobi T
CC BY-SA 3.0
2017-11-10T15:41:26.953Z
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.
llrs
CC BY-SA 3.0
2017-11-10T15:43:40.233Z
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.
Kobi T
CC BY-SA 3.0
2017-11-10T15:48:16.130Z
I think the best thing to do would be to isolate the mutations themselves into a VCF and then there are tons of tools that do that filtering step.  Otherwise I think you will have to write something custom
Bioathlete
CC BY-SA 3.0
2017-11-10T15:51:47.223Z
Is there a python package that create VCF from 2 sequences? It's the first time I hear about this format and it sounds like the right direction.
Kobi T
CC BY-SA 3.0
2017-11-10T16:36:17.777Z
Wait, so you have two lists of strings (sequence) and just want to find those strings present in `file1` that are not in `file2`? What format are your sequences in? Fasta? One per line? Single letter amino acid code? Please show us an example of each of your two files and the output you would expect from that example.
terdon
CC BY-SA 3.0
2017-11-10T16:48:10.643Z
All of those questions have been answered in comments.
Kobi T
CC BY-SA 3.0
2017-11-10T16:58:49.770Z
added an example.
Kobi T
CC BY-SA 3.0
2017-11-10T18:11:41.010Z
A VCF is usually output from a variant caller.  I would suggest aligning all the sequences to a single reference and then you can use samtools to report the differences in a VCF format.
Bioathlete
CC BY-SA 3.0
2017-11-10T18:31:36.830Z
@Bioathlete is there a python package that creates VCF? I use Biopython for alignments.
Kobi T
CC BY-SA 3.0
2017-11-10T19:00:04.383Z
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.
Bioathlete
CC BY-SA 3.0
2017-11-10T19:08:09.277Z
The point is that I want it to be a completely automated tool. As I understand it, samtools is an external tool and does not have a python API.
Kobi T
CC BY-SA 3.0
2017-11-11T02:30:45.883Z
samtools is external but has a python interface via the pysam library.  I am suggesting that you either write something custom or you will have to branch out and use something in addition to BioPython.  The whole process can be automated via a single python script and even as a worst case you can use the python subprocess to access third party tools.
Bioathlete
CC BY-SA 3.0
2017-11-11T10:13:49.297Z
Ok, thanks. Last question: How can I generate a BAM file from 2 sequences in python? Couldn't find it in BioPython or any other package. Is there a package for that?
Kobi T
CC BY-SA 3.0
2017-11-11T17:36:29.430Z
It would be an output format from the aligner you are using if it is supported.
Bioathlete
CC BY-SA 3.0
2017-11-11T21:22:52.010Z
Could you please give a name of package that has this functionality?
Kobi T
<python><proteins><phylogenetics><filtering>
Filtering out all seqs with mutations of list2 from list1
# Filtering out all seqs with mutations of list2 from list1

(Asked by: username_0 on 2017-11-10)

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` = [aa**c**bb**c**, a**b**abbb]
3. Susceptible strains: `list2` = [aa**c**bbb, aaabb**c**]

The expected output would be: [a**b**abbb], since the mutations in `list2` cancels out the item aa**c**bb**c**.

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

**Tags:** python, proteins, phylogenetics, filtering
> **Comments:**
> **username_1**: 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?  
> **username_0**: 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.  
> **username_2**: 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.  
> **username_0**: 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.  
> **username_1**: I think the best thing to do would be to isolate the mutations themselves into a VCF and then there are tons of tools that do that filtering step.  Otherwise I think you will have to write something custom  
> **username_0**: Is there a python package that create VCF from 2 sequences? It's the first time I hear about this format and it sounds like the right direction.  
> **username_3**: Wait, so you have two lists of strings (sequence) and just want to find those strings present in `file1` that are not in `file2`? What format are your sequences in? Fasta? One per line? Single letter amino acid code? Please show us an example of each of your two files and the output you would expect from that example.  
> **username_0**: All of those questions have been answered in comments.  
> **username_0**: added an example.  
> **username_1**: A VCF is usually output from a variant caller.  I would suggest aligning all the sequences to a single reference and then you can use samtools to report the differences in a VCF format.  
> **username_0**: @Bioathlete is there a python package that creates VCF? I use Biopython for alignments.  
> **username_1**: 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.  
> **username_0**: The point is that I want it to be a completely automated tool. As I understand it, samtools is an external tool and does not have a python API.  
> **username_1**: samtools is external but has a python interface via the pysam library.  I am suggesting that you either write something custom or you will have to branch out and use something in addition to BioPython.  The whole process can be automated via a single python script and even as a worst case you can use the python subprocess to access third party tools.  
> **username_0**: Ok, thanks. Last question: How can I generate a BAM file from 2 sequences in python? Couldn't find it in BioPython or any other package. Is there a package for that?  
> **username_1**: It would be an output format from the aligner you are using if it is supported.  
> **username_0**: Could you please give a name of package that has this functionality?  

---

## 1 answers

### Answer by username_4 on 2017-11-14. Score: 1

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)

> **Comments:**
> **username_0**: 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/...)  
> **username_4**: 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  

---

