Defining substring matches

[drhicks]

Hi, I am writing a script to search for motifs in DNA sequence. A user of this program would input at the command line, their fasta file containing the DNA sequence string, and a particular motif to search for. I wish to use the substr() function to search for the particular motif. So far, this is easy enough to do when using A = A, T = T, C = C, and G = G. However, I would like to add the rest of the IUPAC code so that users may define motifs such as ATCNRC, where N = A, T, C, or G and R = A or G? The full IUPAC nucleic acid code is below.

A = '[A]';
C = '[C]';
G = '[G]';
T = '[T]';
R = '[AG]';
Y = '[CT]';
M = '[AC]';
K = '[GT]';
W = '[AT]';
S = '[GC]';
B = '[CGT]';
D = '[AGT]';
H = '[ACT]';
V = '[ACG]';
N = '[ACGT]';

Thanks for any help!

[BrowserUk]

I wish to use the substr() function to search for the particular motif.

#

substr is *not* designed for (nor capable of) searching for anything; so why are you specifying that particular function?

You've defined your IUPAC codes in terms of regex character classes; so why are you eschewing the regex engine?

Given your table, it is trivial to convert IUPAC codes into a regex and use the regex engine to search your fasta file:

my %IUPAC = (
 A => '[A]',
 C => '[C]',
 G => '[G]',
 T => '[T]',
 R => '[AG]',
 Y => '[CT]',
 M => '[AC]',
 K => '[GT]',
 W => '[AT]',
 S => '[GC]',
 B => '[CGT]',
 D => '[AGT]',
 H => '[ACT]',
 V => '[ACG]',
 N => '[ACGT]',
);

my( $file, $motif ) = @ARGV;

my $re = join '', map $IUPAC{ $_ }, split '', $motif;

open FASTA, '<', $file or die $!;
getc( FASTA );                        ## discard first '>'
until( eof( FASTA ) ) {
    chomp( my $id = <FASTA> );                 ## read ident
    my $seq = do{ local $/ = '>'; <FASTA> };
    $seq =~ tr[\n>][]d;
    while( $seq =~ m[($re)]g ) {
        printf "Found: '$1' at '$id':%d\n", $-[0];
    }
}

NB: The above is untested code typed directly into my browser.

[jwkrahn]

until( eof( FASTA) {

You have two left parentheses but only one right parenthesis.

    $seq =~ tr[\n][];

You are replacing the newline character with a newline character?

[Anonymous Monk]

You are replacing the newline character with a newline character?

Indeed, to count newline characters, but not save the count? Hmmm:)

[BrowserUk]

See above:

NB: The above is untested code typed directly into my browser.

Corrected.

[AnomalousMonk]

while( $seq =~ m[($re)]g ) {
    printf "Found: '$1' at '$id':%d\n", $-[0];
}

A question of idle (and perhaps rather trivial) curiosity: In the quoted code, you use  $-[0] "offset of the start of the last successful match" (see  @- in perlvar). My reflexive choice would have been  $-[1] since capture group 1 is being matched. There's no difference in the behavior of the code since capture group 1 is all that's matched, but was there a particular reason you chose as you did?

[BrowserUk]

Um. I cannot remember ever using indexes above 0 with @- or @+. I'm sure I probably have at some point, but I don't remember doing so.

When I first typed it, I didn't use capture brackets and only printed the id and position. I added the capture brackets and '$1' as an afterthought.

[drhicks]

This code works perfect for the job. It is so much simpler and does more than the script I had written using substr(). Thanks Derrick