1

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.

  • 1
    Unix text processing tools are not the best match for this, as they expect to work on "lines" of text, rather than one massive single line. Adding a newline after every TGA or TAA or TAG reduces the problem to "find the first ATG on each line, discarding any line without ATG". – icarus Oct 14 '19 at 8:40
  • 1
    Sigh. My condolences, your professor is not helping you. I've been a bioinformatician for 15 years now and I have never done sequence manipulation in bash. But are you sure they really expect you to code this in bash? Perhaps they just meant you need to use the command line? – terdon Oct 14 '19 at 9:31
  • 1
    Wait, what? grep awk and sed are not bash! Please clarify your requirements. You can do this in bash, I just wrote a script for it, but it's needlessly complicated. Are you now saying we can use non-bash things like awk or sed? – terdon Oct 14 '19 at 9:45
  • 1
    The language is bash. Tools like sed, awk and grep have nothing to do with bash. Awk itself is a full programming language, so allowing you to use awk and not perl is absurd. They are the same class of thing. Even sed is Turing-complete! This sounds like you need to ask your professor to explain what they need. No sane bioinformatician would ever do this in pure bash. But, if you must, I have posted one approach below. It only deals with the forward strand, you can adapt it to deal with the reverse. – terdon Oct 14 '19 at 9:55
  • 1
    bash is just a (powerful) shell; it executes other powerful programs that are built to handle regular expressions and text processing. – Jeff Schaller Oct 14 '19 at 12:32
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
  • BTW, if you don't have grep -P or pcregrep (or perl), see this answer in How to reduce the greediness of a regular expression in AWK? – cas Oct 15 '19 at 3:41
  • Surprisingly, the ? also works with BSD grep (macOS), although it doesn't support PCRE, e.g. grep -Eo 'ATG.*?(TGA|TA[AG]). – Freddy Oct 15 '19 at 4:15
  • @Freddy good to know. i didn't even try testing it on my freebsd vm - though it wouldn't work. – cas Oct 15 '19 at 13:27
0

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
  • Yes I want the output to be the sequences that match the criteria – Mridula M Menon Oct 14 '19 at 8:48
  • I just realized I sound very stupid, can you show me the script you used to frame the sequences with the x and y? – Mridula M Menon Oct 14 '19 at 8:57
  • Okay I think I didn't explain my question well, I want the code to be getting any random sequence, and then outputting all the subsequences that start with ATG and end with TGA or TAA or TAG. I think your code only finds the answer to my specific string – Mridula M Menon Oct 14 '19 at 8:59
  • @MridulaMMenon I have written an echo command but the perl command does the right thing for a file, too. I will modify the answer to make this more clear. Your question did not specify where data comes from or where the filtered data goes so I could only guess. – Ned64 Oct 14 '19 at 9:03
  • I used your cat statement, and it said head: illegal line count -- -1 – Mridula M Menon Oct 14 '19 at 9:14
-1

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
  • You are a God send. I figured out how to do it for the reverse strand too. Now I just have to hope my prof doesn't look up what I turn in because then he will find this here. – Mridula M Menon Oct 14 '19 at 10:35
  • @MridulaMMenon if your prof isn't an idiot they will find this. And they will also tell you this is a stupid way of doing it. You are there to learn, just blindly copying answers you don't understand, especially bad answers like this one, will only harm you. Don't copy this, learn from it and adapt it. – terdon Oct 14 '19 at 10:40
  • My prof is a biochemist who doesn't know too much about programming. However, I wouldn't put it past them to look this up. I'll adapt it for sure, but the logic and the syntax will largely be the same. I'm worried that they will see the date this question was posted and connect the dots. Also, there isn't a question on the internet with this specific requirement so... – Mridula M Menon Oct 14 '19 at 10:53
  • @MridulaMMenon well, you should have thought of that before you decided to ask strangers to do your work for you. It hadn't really occurred to me that you would just copy/paste this. If I had your professor's email, I would inform them myself. You should take this opportunity to learn, not just get a grade. – terdon Oct 14 '19 at 10:58
  • I can save you the trouble and tell the professor myself, I feel like I know how to get this done in other programming languages like C++ so I have achieved the objective, if my prof wants me to do it some particular way, while giving me bare minimal tools to do it, then it's an exercise in futility. The other classmates of mine just wrote a tediously long code without using loops so that beats the point of using bash but he has to be okay with that because he hasn't taught us enough syntax to write something like what you wrote which I would argue is more elegant – Mridula M Menon Oct 14 '19 at 11:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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