0

I'm new to bash, and I found this script to identify ORF, the problem I have is that it gives me three results (ATGATGGATCCTTGCGTGCAGTAA, ATGGGATAA and ATGATGGATCCTTGCGTGCAGTAA) and not two as it should be (ATGATGGATCCTTGCGTGCAGTAA and ATGGATCCTTGCGTGCAGTAA).

example sequence : TTAATGATGGATCCTTGCGTGCAGTAACCCATGGGA

#!/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

`

/home/ibioic/atg.sh TTAATGATGGATCCTTGCGTGCAGTAACCCATGGGA
ORF: ATGATGGATCCTTGCGTGCAGTAA
ORF: ATGGGATAA
ORF: ATGATGGATCCTTGCGTGCAGTAA

`

New contributor
Natalia is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
5
  • What's an ORF? It looks like you might be doing bioinformatics. If so please explain your requirements in computational terms so we can all understand what you're asking
    – roaima
    May 14 at 17:47
  • Hi, You must find all 'ORF's in the sequence, an ORF is any region that starts with ATG , is followed by any number of characters multiple of 3 and ends with TGA or TAA or TAG.
    – Natalia
    May 14 at 18:02
  • 2
    "You must find..." please edit your question to provide clarifications. Put all the information in one place and then it's easier for people to read. Make it easy for people to help you
    – roaima
    May 14 at 18:22
  • Note that there is a StackExchange site dedicated specifically to bioinformatics. It's at bioinformatics.stackexchange.com
    – Kusalananda
    May 15 at 10:24

1 Answer 1

0

Not clear what you're trying to do, but if it's about breaking the input into groups of 3 letters (codons?) and finding the sequences of those that start with ATG and end with the next occurrence after that of TAA, TAG or TGA as you sort-of clarified in comments, with pcregrep (perl compatible regular expression g/re/p):

$ pcregrep -o1 -o2 '\G(?:...)*?(ATG)(?=((?:...)*?T(?:AA|AG|GA)))' <<< TTAATGATGGATCCTTGCGTGCAGTAACCCATGGGA
ATGATGGATCCTTGCGTGCAGTAA
ATGGATCCTTGCGTGCAGTAA

Or directly with perl:

$ perl -lne 'print $1.$2 while /\G(?:...)*?(ATG)(?=((?:...)*?T(?:AA|AG|GA)))/g' <<< TTAATGATGGATCCTTGCGTGCAGTAACCCATGGGA
ATGATGGATCCTTGCGTGCAGTAA
ATGGATCCTTGCGTGCAGTAA

We're capturing the second to last codons in the match as part of a look ahead operator (captured in $2) so as to resume searching after the first codon to find potentially overlapping matches as seem to be the intention given your expected output. However note that in

ATG111ATG222TAA333TAA444ATG555TAA

While it will report ATG111ATG222TAA, ATG222TAA, and ATG555TAA, some of which overlap, it won't report ATG111ATG222TAA333TAA, ATG111ATG222TAA333TAA444ATG555TAA, ATG222TAA333TAA, ATG222TAA333TAA444ATG555TAA as it stops at the first occurrence of a TAA/TAG/TGA.

That \G(?:...)*? which matches on any amount (as few as possible) of character triplets after the end of the previous match (or the start of the line for the first match) is to make sure the ATG start is aligned on a 3 letter boundary, so that for instance it matches inside 123ATGxxxTGA but not 12ATGxxxTGA. If that's not what you want, just drop it.

3
  • Thanks Stephane, but I have to do it with the for loop
    – Natalia
    May 14 at 19:45
  • 2
    @Natalia why must you use a specific construct? Is this homework?
    – roaima
    May 15 at 10:10
  • @Natalia The loop is implicit in the tool.
    – Kusalananda
    May 15 at 10:26

Not the answer you're looking for? Browse other questions tagged or ask your own question.