Shell Programming and Scripting

Hi,

I'm working with DNA sequences, which as you might know are composed of "codons", which are the "words" if you like, and all codons are 3 letters long. So something I frequently need in my scripts, is pattern searching that respects this typical word size.

E.g. : searching "TAC" in the string "ACGTACGTACGT" with sed would return 2 matches, while there is only one codon "TAC". Respecting the "reading frame" means that words can only start at positions 0, 3, 6, 9, etc.

As the sequences can be enormously long, is there an efficient way that command-line tools like sed or awk can handle this kind of reading ? GAWK for example has a FIELDWIDTHS variable, but it needs a list of fieldwidths, just specifying "3" doesn't work, does it ?

Otherwise, I have to turn to more specific bioinformatic software, which is less easy to incorporate in my scripts.

Thank you in advance for any thoughts on this.
Jos

Do you just need the count?

Code:

$ echo "ACGTACGTACGT" | perl -nle 'while(/(.{3})/g){ $total++ if($1 eq "TAC");} END{print $total}'
1
$ echo "ACGTACGTACGTTAC" | perl -nle 'while(/(.{3})/g){ $total++ if($1 eq "TAC");} END{print $total}'
2

Code:

$ echo "ACGTACGTACGT" | fold -w3 | grep -c TAC
1

From a file:

Code:

fold -w3 file | grep -c TAC

This might work too:

Code:

perl -0777 -ne '@m=/.{3}TAC/g;push @m,"1" if (/^TAC/); print $#m+1' file

Thanks folks !

Actually, counting is just one of the things I'd like to do frequently ; replacing is another.

I can see quite easily how I could do this using the "fold" command (thanks for bringing it to my attention!)

Code:

echo "ACGTACGTACGT" | fold -w3 | sed -e 's/TAC/NNN/g' | tr -d '\n'

Maybe the perl solutions are quicker ? I'll have to test that, but I'll have to know perl a bit better.

Thanks again !
jos

Depending on what your overall goal is, you might want to take a look at BioPerl.

Yes, I know BioPerl exists, I even installed it once.

For the moment, I hardly know perl, and simple tasks like pattern matching, changing file formats, etc are things one can do with awk, sed. But yes, eventuelly I'll get to try bioperl.

jos