make fasta sequence which is multiple of three
2
0
Entering edit mode
3.5 years ago
alim.hcu ▴ 10

I have fasta sequence of different length. But I want to make number of nucleotide in fasta sequences multiple of three by adding AA or -- at the end. Like if the sequence length is 149 then i want to make it 150 by adding -- or any nucleotides.

Thank You

alignment • 2.1k views
ADD COMMENT
3
Entering edit mode

This looks like an XY Problem. Why do you want to do this multiple of 3 thing? The 3 part tells me if probably has to do with codon usage, but in no scenario can I imagine needing this functionality. Can you explain your end goal please?

ADD REPLY
0
Entering edit mode

Hello, if you can program in Java, you may use the SEDA API (https://github.com/sing-group/seda) in order to create this operation that you need. You can even use the new Java 9 shell to easily do this as I show in this post:

Regards.

ADD REPLY
0
Entering edit mode

Thanks for you reply. I downloaded this tools. But the video is not much clear. which option i should use to make multiple of three nucleotides.

ADD REPLY
0
Entering edit mode

Actually I have to calculate dn/ds using paml. And it does not take any sequences which is not multiplle of three.

ADD REPLY
3
Entering edit mode

I have to calculate dn/ds using paml

Then I strongly advise against adding random nucleotides.

ADD REPLY
0
Entering edit mode

You should have said so in the original question to avoid half the comments you received here :-).

ADD REPLY
0
Entering edit mode

Please use ADD COMMENT or ADD REPLY to answer to previous reactions, as such this thread remains logically structured and easy to follow. Your reaction was moved but as you can see it's not optimal. Adding an answer should only be used for providing a solution to the question asked.

ADD REPLY
0
Entering edit mode

If you are planning on using PAML then you'll need to pairwise align those sequences anyway. So the more cleaver thing to do is to "codon-aware" align the nucleotide sequence and strip off anything that is not aligned (== will automatically result in modulo3 'sequences')

ADD REPLY
0
Entering edit mode

Ok.. I got it..is there anything to do this in paml. I am aligning those sequences using maft. How should i trim those sequences.

ADD REPLY
0
Entering edit mode

Start by searching the forum for MAFTT posts and if you don't find an appropriate post, ask a new question.

ADD REPLY
0
Entering edit mode

ah, no you will have to do this prior to providing them to PAML.

In my lab we have custom script to do this (not helping much :-) ), but here is a tool that does it as well : trimAl and there will be others around

make sure that you do codon-align in mafft (this is crucial for dN/dS estimations!!)

ADD REPLY
0
Entering edit mode

Specifically "AA" or any nucleotide?

ADD REPLY
0
Entering edit mode

Thank You very much..it worked..

ADD REPLY
0
Entering edit mode

Wouter literally just asked you to not add answers. Please be more careful. This should be a comment or comment reply. I've moved it to a comment on the top level post, but it ideally belongs elsewhere. Please read the posts under http://biostars.org/t/how-to to structure your posts better.

ADD REPLY
1
Entering edit mode
3.5 years ago

Beside the correct question from Ram why you want to do this, here is a python solution:

from Bio import SeqIO

original_file = "./original.fasta"
extended_file = "./extended.fasta"

with open(original_file) as original, open(extended_file, 'w') as extended:
    records = SeqIO.parse(original_file, 'fasta')
    for record in records:
        if len(record.seq)%3 == 1:
            record.seq = record.seq + "--"
        elif len(record.seq)%3 == 2:
            record.seq = record.seq + "-"
        SeqIO.write(record, extended, 'fasta')

fin swimmer

ADD COMMENT
1
Entering edit mode

Minor comment, but there is no need to with open(...) your input file when using SeqIO

ADD REPLY
0
Entering edit mode

I want to run this script in loop..Please suggest me to modify it for running in this script in loop..

Thank You

ADD REPLY
3
Entering edit mode

why are you still insisting in fiddling about with your input sequences? while the solution offered by finswimmer will work perfectly , several other people have pointed out to you that this is not the way you need to go!

ADD REPLY
0
Entering edit mode

Actually i am working on transposons. And i have to calculate insertion time for LTR pairs present in LTR elements. So basically these are the nucleotide sequences.

ADD REPLY
1
Entering edit mode

Actually i am working on transposons.

Why did we have to wait 20 hours and so many comments before you tell us what you are working on? Don't you think all of this thread is a big waste of time?

ADD REPLY
0
Entering edit mode

alright, I was kinda expecting this already.

anyway the point remains: you should not add bases to get to modulo3 sequences, You first have to align them and subsequently remove unaligned parts (rather than adding bases!!).

Moreover have you looked into the literature on estimating transposon insertion time? There will be software around that specifically does this kind of analysis (though I can not think of one from the top of my head)

ADD REPLY
0
Entering edit mode

"You first have to align them and subsequently remove unaligned parts (rather than adding bases!!" I would be interested in doing this how to do this? Since i was trying to do this "If you're intending to run it through PAML to do dNdS " but again the multiple of 3 errors comes.

ADD REPLY
0
Entering edit mode

have a look at TrimAl , or such ... in general tools to clean multiple seq alignments

or manual if there are not that much sequences?

ADD REPLY
0
Entering edit mode

thanks will look into it trimAI

ADD REPLY
1
Entering edit mode
3.5 years ago
Joe 19k

It can be done pretty easily with BioPython, as FinSwimmer has shown, particularly with the MutableSeq objects.

If you're intending to run it through PAML to do dNdS etc, then a far more common approach, rather than hacking at your sequences, is just to use a dedicated codon-aware aligner like the imaginatively named CodonAlign

EDIT:

For clarity from the comments we've all made, this should be your workflow:

  1. I'd strongly suggest creating files with all pairwise comparison sequences you're interested in. If you have many sequences, this may rapidly become untenable though, so think very carefully about what you're trying to show/achieve.

    • Its an n choose k problem so, all pairwise comparisons of even just 100 sequences is 4950.
  2. CRUCIALLY, use a codon-aware aligner as we've mentioned. Do not try to hack your sequences to a multiple of 3 by adding random characters as this will affect your results.

  3. Santise your alignments if necessary. This you may need to do manually.

  4. Run PAML.

ADD COMMENT
0
Entering edit mode

Can You describe the pipe line. I have thousand of sequences in pair wise alignment. I want to calculate dn/ds of each pairwise alignment in paml. If any script you can suggest, it would be very helpful for me.

Thank You

ADD REPLY
0
Entering edit mode

The program is supposed to be downloadable from here: http://homepage.mac.com/barryghall/CodonAlign.html but the webpage seems to be down.

EDIT: You can download it here: https://github.com/ozagordi/CodonAlign OOPS

ADD REPLY
0
Entering edit mode

To my knowledge (limited as it is in this instance), you can't do this with pairwise alignments, they will need to be multiple sequence alignments, else there isn't sufficient information for the dNdS rates to be meaningful. I could be totally wrong here though.

As for scripts, no I don't know I'm afraid. But CodonAlign, and tools like it, are self contained executables, so you shouldn't need much if any scripting in order to use them.

If I am wrong, then it will probably be as simple as concatenating all your sequence pairs in to files of their own, and then running CodonAlign and PAML for each file in a loop or parallel.

ADD REPLY
0
Entering edit mode

I have to disagree here: it's perfectly OK to calculate dN and/or dS from pairwise alignments (this is how all these Ks plots are being generated).

Otherwise the approach you describe is exactly what I would do :

  1. create set of files with each a pair of sequences in it to be aligned
  2. align the sequences with something like CodonAlign (codon-aware aligner is key here!!)
  3. clean the alignments (to remove gaps etc)
  4. run PAML (codeml) on each alignment result file
ADD REPLY
0
Entering edit mode

I stand corrected on the pairwise in that case (I've never done it myself), but this is absolutely the workflow, so we all agree there!

I'm curious though, if you're just using pairwise sequences, how can you decide which sequence is 'right' relative to the other? You'd have no other information about what the most common codon in that position is to know whether a particular sequence is drifting/being selected for or against?

(this is my pure naivety though I'm sure)

ADD REPLY
1
Entering edit mode

Indeed, in a pairwise context you can't tell, but what you are referring to is already step 2 of this kind of procedure. step 1 is to calculate (estimate rather) the dS and dN for each pair. We're not looking for the most common codon or such yet, we're simply looking at the divergence on nucleotide level between 2 genes that once shared a common ancestor (origin).

In step 2 we can start looking at the ration of dS over dN which (coming to think of it) we can actually also do in a pairwise manner. The thing we are looking for is to see whether there are more non-synonymous or more synonymous substitutions happened dN/dS > 1 is positive (driving) selection , = 1 is neutral, < 1 is purifying selection . There should in "theory" not be any bias between codons used , we assume that mutations are happening at a constant rate (kuch ;) ) and thus we are only interested to see if the mutations we observe are more dN or are more dS

ADD REPLY

Login before adding your answer.

Traffic: 1884 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6