Perl Script To Get Upstream And Downstream Regions Of Start Codons

Question

[eddie.im]

Hi,

I'm trying to get upstream and downstream regions of start codons in a genome (20nt upstream and 4nr downstream) but my script gets only the first start codon even if I use the "g" modifier on the regex. How can I get it to read all start codons (ATG)?

use strict;
use warnings;
my @regions;
my $term="ATG";
my $seqsample="CCCCATAGAGATAGAGATAGAGAACCCCGCGCGCTCGCATGGGGATGCATGATTCGG";
while ( $seq =~ m/(\S{20})$term(\S{4})/g ) {
        my $xx = $1.$term.$2;
    push (@regions, $xx);
    }
print "@regions\n"

Here's the script that I wrote.

Thanks in advance.

[harpreetmanku04]

Hello I am also trying this, but instead of ATG I need to extract region from sequence position suppose in long sequence of 5000 nucl I want to extract +/- region from position

1------------------------------------------2946_GENE_SEQUENCE_2995-------------------------------------5000
                                                  200 residues here                                                200 residues here

[harpreetmanku04]

>gi|546867469
TAGCTACATCAGTTGAAGATGGGGTTGCTGGGTAGCTGTCAGGCATAGGAGTCCCATCTTTCTCGGAAAG
AATGACAATAAGGTTGTTAAAGGGCAATCAGTCAAGGACAGGGGGAGCCACTCAAATCAGCACATGTGGC
ATGACAGTGATGCAGCATTCCCCTCAACCACCAATTGCTCCCAGGAATAGCACCTGGCTGGTCTGGCATC
TCGGGATGAGGGCTAAGGCACAGTGCCATTAGTGCCAAGAACCAAAGATCTGAGTCTTTGCAGACCCTAA
AAATTAAATCCTTTCTGTGCCCCGGGTGGAGGGGGAGAGTTCAAGTCACTTTATAATATTTGCCTCATAT
ACATATCACCTGTCCCCTTCCACCTCCTTTTTCTTGCAGTCCTAATTATCCAAGTAGCATCTACAACTCT

region to extract from

150 to 160 20-20 up/down stream

Kindly help me to extract this region with perl script or awk command.

Thank you

Answer 1

[swbarnes2]

First

while ( $seqsample =~ m/(\S{20})$term(\S{4})/g ) {

Second, this construction of the variable flanked by other stuff just isn't working. I'd use pos to get the positions of ATG in the string, then use substr to pull out the flanking sequence.

[Kenosis]

The OP only needs to use zero-width positive lookbehind and positive lookahead assertions, and capture their contents:

while ( $seqsample =~ /(?<=(\S{20}))$term(?=(\S{4}))/g ) { ...

Answer 2

[JC]

The errror is because you are extending the ATG 20 bases before and 4 after, the regex is matching non-overlapping hits. Its better to find the ATG positions and get the string with substr as swbarnes2 suggested:

use strict;
use warnings;
my @regions;
my $term = "ATG";
my $len  = length $term;
my $seq  = "CCCCATAGAGATAGAGATAGAGAACCCCGCGCGCTCGCATGGGGATGCATGATTCGG";
while ( $seq =~ m/$term/g ) {
    my $pos = $+[0];
    my $xx = substr ($seq, $pos - 20 - $len, 20 + $len + 4);
    push (@regions, $xx);
}
print join "\n", @regions, "";

Output:

$ perl atg.pl 
AGAGAACCCCGCGCGCTCGCATGGGGA
CCCCGCGCGCTCGCATGGGGATGCATG
GCGCGCTCGCATGGGGATGCATGATTC

[Kenosis]

The OP could just capture the contents of term-flanking zero-width assertions.

Answer 3

[Kenosis]

You're so very close! Just capture the contents of positive lookbehind and positive lookahead, zero-width assertions, so the search position isn't impacted by them:

use strict;
use warnings;

my ( $term, @regions ) = "ATG";
my $seq = "CCCCATAGAGATAGAGATAGAGAACCCCGCGCGCTCGCATGGGGATGCATGATTCGG";

while ( $seq =~ /(?<=(.{20}))$term(?=(.{4}))/g ) {
    push @regions, $1 . $term . $2;
}

print "$_\n" for @regions;

Output:

AGAGAACCCCGCGCGCTCGCATGGGGA
CCCCGCGCGCTCGCATGGGGATGCATG
GCGCGCTCGCATGGGGATGCATGATTC

Hope this helps!

[eddie.im]

Thanks for your help! it worked very well!

[Kenosis]

You're most welcome!

Answer 4

[cyril-cros]

Another solution, without Perl:

grab the official annotation in gtf format, it should contains some lines with start_codon. Use grep to extract them, gtf2bed (from BEDOPS) and bedtools slop -l 20 -r 4 -s (from bedtools).