How to trim nucleotide sequences to length divisible by 3

I want to strip predicted nucleotide protein coding sequences of the last 1 or 2 bases which represent incomplete codon.

The input is multifasta consisting of nucleotide protein-coding sequences in reading frame 1. Some of them are trimmed at various positions at the 3'end resulting in 1-2 3'end nucleotides that does not form a codon (these are ignored when I translate the nucleotide sequences to protein sequences; the extra nucleotides are then source of problems in my further workflow).

I tried to extract the longest ORF from the nucleotide sequences to overcome this, however, it does not work in all cases (in some of the input sequences there are in-frame stop codons; other drawback was that using e.g. hmmer2go I could not restrict the search for ORFs only in frame +1 and thus for some of the fastas, alternative ORFs in reverse sense, +2 or +3 reading frames were retrieved as the longest ORFs)

I'd like to ideally use tools such as awk (and its parameter % - modulo) to trim down the sequences from 3'end to the longest sequence divisible by 3 - however, I did not manage to make anything to work. Any other solutions are, however, most welcomed as well.

awk incomplete codon

• 4,051 views • link

written 8.6 years ago by al-ash

1 answer

I'm not sure what is the correct protocol for answering my own question.

In the end, my input were multiple fasta files each including one sequence. I used bioawk and the awk modulo function (%) to create 3 cathegories of sequences: with sequence length divisible by 3, divisible by 3 with remainder 1 and with remainder 2. Sequences divisible with remainder were trimmed from the end by one respectively two characters.

#!/bin/bash
workdir=mydir

for file in ${workdir}/nt_renamed/* 
do filebasename=$(basename "$file" )
# if modulo 3 of the sequence length is 0, print the fasta header and the sequence and output it to file 
bioawk -c fastx 'length($seq) % 3 == 0  {print ">"$name; print $seq}' $file  >  ${workdir}/nt_renamed_bioawked/${filebasename}
# if modulo 3 of the sequence length is 1 or 2, print the fasta header and the sequence, trim lines which does not start with ">" from its end by one resp. 2 characters
bioawk -c fastx 'length($seq) % 3 == 1  {print ">"$name; print $seq}' $file | sed '/^>/! s/.$//'  >> ${workdir}/nt_renamed_bioawked/${filebasename} 
bioawk -c fastx 'length($seq) % 3 == 2  {print ">"$name; print $seq}' $file | sed '/^>/! s/..$//'  >> ${workdir}/nt_renamed_bioawked/${filebasename} 
done

• 0 views • link

written 8.6 years ago by al-ash

This is perfect. I will toggle this answer as accepted to provide closure to this thread.

• 0 views • link

written 8.6 years ago by GenoMax