I am looking for a tool that I can use to remove sequences from my protein alignments that do not contain highly conserved residues. For example, I would be able to direct this tool to remove any sequence that does not contain an arginine at position 124 in my supplied alignment. Is anyone aware of a tool that can help me with this goal?
Just a bit of additional clarification: I am sequencing an amplicon of a functional gene and generating protein sequence alignments using RDP's fungene pipeline. I want to further screen the sequences by eliminating any that do not contain a selection of conserved residues at various positions.
Thanks!
I would strongly recommend to use Python and Biopython. In Biopython, there are some libraries to manipulate alignment, count sites, assign position. In a few lines, you could be able to do what you want. And even more.
You can do this with a simple selections, depending on your platform the code will vary.
Here are the steps I would take.
1st: renumber all amino acids based on the alignment, to define a parameter that establishes an easy to use residue by residue correspondence.
2nd: select all sequences that do not contain an arginine at 174, or any residue you do not want. Your software should have documentation for it's particular selection syntax which should use logical operators such as "and", "not", "or" etc
3rd: delete selected sequences
I would strongly recommend to use Python and Biopython. In Biopython, there are some libraries to manipulate alignment, count sites, assign position. In a few lines, you could be able to do what you want. And even more.
I strongly recommend to use a simple code in Perl, Python (even c++), too.
In this way, you do not face to any problem in order to manipulate big data.
Perl will do it nicely. Here follows a possible logic.
#1 get the aln into a hash, %ALN (seqname => sequence)
$x = n; # n = position in aln of residue of interes.
$reswanted = 'X'; # where X = (A .. Y) residue you want to be present.
foreach $seqname (keys %ALN){
my $sequence = $ALN{$seqname};
my $res = substr ($sequence, $x-1,1);
if ($res neq $reswanted ){
delete $ALN{$seqname};
}
}
## Then you can print the hash
Thank you all for the replies.
I used biopython to solve my problem. Here is the python program I wrote (it's my first attempt at python so any critiques are welcome):
from Bio import AlignIO
from Bio.Align import MultipleSeqAlignment
alignment = AlignIO.read("MyProt.fasta", "fasta")
goodseqs1 = MultipleSeqAlignment([])
goodseqs2 = MultipleSeqAlignment([])
badseqs = MultipleSeqAlignment([])
for sequence in alignment:
if sequence.seq[48] == "F":
goodseqs1.append(sequence)
elif sequence.seq[48] == "Y":
goodseqs1.append(sequence)
else:
badseqs.append(sequence)
for sequence in goodseqs1:
if sequence.seq[46] == "Q":
goodseqs2.append(sequence)
elif sequence.seq[46] == "R":
goodseqs2.append(sequence)
elif sequence.seq[46] == "G":
goodseqs2.append(sequence)
elif sequence.seq[46] == "I":
goodseqs2.append(sequence)
else:
badseqs.append(sequence)
AlignIO.write(goodseqs2, "SCREENED_SEQS.FASTA", "fasta")
AlignIO.write(badseqs, "Discard.fasta", "fasta")
print("Alignment length %i" % alignment.get_alignment_length())
print("initial # of seqs", len(alignment))
print("seqs passing first screen:", len(goodseqs1))
print("final # of seqs", len(goodseqs2))
print("output files are: SCREENED_SEQS.FASTA and Discard.fasta")
You can do in one screen I think (I didn't try):
from Bio import AlignIO
from Bio.Align import MultipleSeqAlignment
alignment = AlignIO.read("MyProt.fasta", "fasta")
goodseqs = MultipleSeqAlignment([])
badseqs = MultipleSeqAlignment([])
for sequence in alignment:
if sequence.seq[46] in ["Q","R","G","I"] and sequence.seq[48] in ["F","Y"] :
goodseqs.append(sequence)
else:
badseqs.append(sequence)
AlignIO.write(goodseqs, "SCREENED_SEQS.FASTA", "fasta")
AlignIO.write(badseqs, "Discard.fasta", "fasta")
print("Alignment length %i" % alignment.get_alignment_length())
print("initial # of seqs", len(alignment))
print("final # of seqs", len(goodseqs))
print("output files are: SCREENED_SEQS.FASTA and Discard.fasta")
Yes you are correct Romain! Your version worked perfectly. Thank you very much for your help.
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics
- Bioinformatic Tools