Trying to write a script to find ORF from 6 different frames

Question

So basically, you have a sequence, say

AAAGCATATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAG 

You need to find all the 'ORFs' in the sequence, an ORF is any region that starts with ATG and ends in TGA or TAA or TAG.

So, for example, in the above sequence, the ORFs are

AAAGCAT**ATGCTAG**CCCGTATAGCGATACTAGCTATACGATATAT**ATGATCAATGCCCGTATAG** 

You'll notice, in the second ORF, there is an ATG inside the ORF, this happens when there isn't a TGA or TAA or TAG before it ends.

So yeah, basically that's the question. I know 2-3 ways of doing it on C++, but in terms of syntax in bash I have no clue.

I can't use any libraries or perl or anything like that, no special functions, just stuff like grep, awk, sed and loops.

Answer 1

You can use grep if you are using the GNU version of grep, which has the -P option for perl-compatible regular expressions (PCRE).

Alternatively, you could use pcregrep (aka pgrep) by the author of the PCRE library. There's not much reason to use this these days unless you are using an old version of GNU grep, or a non-GNU grep, that doesn't suppoer the -P option and can't be upgraded or replaced.

e.g. (assuming that the sequences are in a file called input.txt):

$ grep -oP 'ATG.*?TA[AG]'  input.txt 
ATGCTAG
ATGATCAATGCCCGTATAG

The -o option tells GNU grep to output only the matching text, not the entire line, and the -P tells it to use perl-compatible regular expressions.

Or, if you're using pcregrep:

$ pcregrep -o 'ATG.*?TA[AG]'  input.txt 
ATGCTAG
ATGATCAATGCCCGTATAG

The non-greedy ? modifier in the regexp (.*?) ensures that it captures ALL matching patterns, not just the longest one. In the context of regular expressions, "greedy" means "try to match as MUCH as possible" (the default), and "non-greedy" means "try to match as LITTLE as possible".

There's a good explanation of this at https://www.regular-expressions.info/repeat.html. BTW, the rest of that site is a pretty good place to learn about regular expressions, with lots of tutorials and examples.

Note that most regular expression libraries do not implement non-greedy matches, it is a perl extension that has also been adopted by GNU grep. and programs linked with PCRE, of course.

BTW, This is what the output would look like without the non-greedy modifier:

$ grep -oP 'ATG.*TA[AG]'  input.txt 
ATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAG

Answer 2

The searching and matching to a regular expression can be done with perl (my favourite sed does not support the needed non-greedy regex(p) match):

# echo AAAGCATATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAG | perl -pe 's|(.*?)(ATG.*?TA[AG])(.*?)|\2\n|g'
ATGCTAG
ATGATCAATGCCCGTATAG
AA

where the perl command prints the sequences you requested - and one last line with the rest of the line. In order not to make the regular expression more complex this can easily be removed by hand or with head:

echo AAAGCATATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAGAA | perl -pe 's|(.*?)(ATG.*?TA[AG])(.*?)|\2\n|g' | head -n -1
ATGCTAG
ATGATCAATGCCCGTATAG

These echo commands show the correct result for your sequence. If you wish to read the contents from a file, and send the result to a second file, do this:

cat original_file | perl -pe 's|(.*?)(ATG.*?TA[AG])(.*?)|\2\n|g' | head -n -1 > new_file

where original_file is your source file and new_file the target with filtered patterns.

The above uses head from GNU coreutils with the "-n -1" syntax. If that does not work for you, please try

cat original_file | perl -pe 's|(.*?)(ATG.*?TA[AG])(.*?)|\2\n|g' | awk 'NR>1 {print prev} {prev=$0}' > new_file

Answer 3

This is possible in bash, but it's really not a good tool for it:

#!/bin/bash 

# Read the sequence into the variable $seq
seq=$1

## Check all three frames
for ((frame=0; frame<=3; frame++)); do
  ## Read the sequence in groups of 3
  for ((i=$frame;i<=${#seq};i+=3)); do
    ## The codon: three nucleotides starting from the current position.
    codon=${seq:i:3}
    ## set isORF to 1 if this is an ATG
    if [[ ${seq:i:3} = "ATG" ]]; then 
      isORF=1
    fi
    ## If we're in an ORF
    if [[ $isORF = 1 ]]; then
      ## Add this codon to the ORF's sequence
      ORF="${ORF}${codon}"
      ## Is this a STOP?
      if [[ ${seq:i:3} = "TGA" ||
            ${seq:i:3} = "TAA" ||
            ${seq:i:3} = "TAG" ]];
      then
        ## Print the ORF
        echo "ORF: $ORF"
        ## Set isORF to 0 and empty the ORF variable to repeat the process
        isORF=0
        ORF=""
    fi
  fi
  done
done

Save that as foo.sh, make it executable (chmod a+x foo.sh) and run it like this:

/path/to/foo.sh AAAGCATATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAG