Trouver tous les répétées 4-mères dans une séquence d'ADN - Perl

Question

Bonjour,

J'essaie d'écrire un programme qui lit dans un fichier au format FASTA contenant des séquences d'ADN multiples, identifie tous les répétées 4-mères (c.-à-tout 4-mères qui se produisent plus d'une fois) dans une séquence, et imprime le 4-mer répété et l'en-tête de la séquence dans laquelle il a été trouvé. Un k-mer est simplement une séquence de k nucléotides (par exemple, "aaca", "gacg" et "tttt" sont des 4-mères).

Voici mon code:

use strict; 
use warnings; 

my $count = -1; 
my $file = "sequences.fa"; 
my $seq = ''; 
my @header =(); 
my @sequences =(); 
my $line = ''; 
open (READ, $file) || die "Cannot open $file: $!.\n"; 

while ($line = <READ>){ 
    chomp $line; 
    if ($line =~ /^>/){ 
     push @header, $line; 
     $count++; 
     unless ($seq eq ''){ 
      push @sequences, $seq; 
      $seq = ''; 
     } 
    } else { 
     $seq .= $line; 
    } 
} push @sequences, $line; 

for (my $i = 0; $i <= $#sequences+1; $i++){ 
    if ($sequences[$i] =~ /(....)(.)*\g{1}+/g){ 
     print $header[$i], "\n", $&, "\n"; 
    } 
} 

J'ai deux demandes: D'abord, je ne sais pas comment concevoir mon modèle regex pour obtenir la sortie désirée. Ensuite, et c'est moins important, je suis certain que mon code est très inefficace, donc s'il y a un moyen de le raccourcir, dites-le moi s'il vous plaît.

Merci d'avance!

Voici un exemple pour un fichier FASTA: (Notez qu'il ya une ligne supplémentaire entre les séquences, ce qui est le cas dans les fichiers FASTA d'origine)

> NC_001422.1 entérobactéries phage phiX174 sensu lato, complète génome GAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTttttttCGGATATTTCTGATGAGTCGAAAAAT CCCTTACTTGAGGATAtatataAATTATGTCTAATATTCAAACTGGCGCCGAGCGTATGCCGCATGACCT

> NC_001501.1 entérobactéries phage phiX184 sensu lato, génome complet AACGGCTGGTCAGTATTTAAGGTTAGTGCTGAGGTTGACTACATCTGTTTTTAGAGACCC AGACCTTTTA TCTCACTTCTGTTACTCCAGCTTCTTCGGCACCTGTTTTACAGACACCTAAAGCTACATCGTCAACGTTA TATTTTGATAGTTTGACGGTTAATGCTGGTAATGGTgagagagaGGTTTTCTTCATTGCATTCAGATGGA TCAACGCCGCTAATCAGGTTGTTTCTGTTGGTGCTGATATTGCTTTTGATGCCGACCCTAAATTTTTTGC CTGTTTGGTTCGCTTTGAGTCTTCTTCGGTTCCGACTACCCTCCCGACTGCCTATGATGTTTATCCTTTG

> NC_001622.5 entérobactéries phage phiX199 de Lato, génome complet de TTCGCTGAATCAGGTTATTAAAGAGTTGCCGAGATATTTATGTTGGTTTCATGCGGATTGGTCGTTTAAA TTGGACTTGGTGGCAAGTCTGCCGCTGATAAAGGAAAGGATAATGACCAAATCAAAGAACTCGTGATTAT CTTGCTGCTGCATTTCCTGAGCTTAATGCTTGGGAGCGTGCTGGTGCTGATGCTTCCTCTGCTGGTATGG TTGACGCCGGATTTGAGAATCAAAAATGTGAGAGAGCTTACTAA AATGCAACTGGACAATCAGAAAGAGA GATGCAAAATGAGACTCAAAAAGAGATTGCTGGCATTCAGTCGGCGACTTCACGCCAGAATACGAAAGAC CAGGTATATGCACAAAATGAGATGCTTGCTTATCAACAGAAGGAGTCTACTGCTCGCGTTGCGTCTATTA TGGAAAACACCAATCTTTCCAAGCAACAGCAGGTTTCCGAGATTATGCGCCAAATGCTTACTCAAGCTCA AACGGCTGGTCAGTATTTTACCAATGACCAAATCAAAGAAATGACTCGCAAGGTTAGTGCTGAGGTTGAC TTAGATGAGTGTTCATCAGCAAACGCAGAATCAGCGGTATGGCTCTTCTCATATTGGCGCTACTGCAAAG

Answer 1

Je serais probablement aborder votre problème un peu plus comme ceci:

#!/usr/bin/env perl 

use strict; 
use warnings; 

use Data::Dumper; 

#set paragraph mode. Iterate on blank lines. 
local $/ = ''; 

#read from STDIN or a file specified on command line, 
#e.g. cat filename_here | myscript.pl 
#or myscript.pl filename_here 
while (<>) { 
    #capture the header line, and then remove it from our data block 
    my ($header) = m/\>(.*)/; 
    s/>.*$//; 

    #remove linefeeds and whitespace. 
    s/\s*\n\s*//g; 
    #use lookahead pattern, so the data isn't 'consumed' by the regex. 
    my @sequences = m/(?=([atcg]{4}))/gi; 

    #increment a count for each sequence found. 
    my %count_of; 
    $count_of{$_}++ for @sequences; 

    #print output. (Modify according to specific needs. 
    print $header,"\n"; 

    print "Found sequences:\n"; 
    print Dumper \@sequences; 
    print "Count:\n"; 
    print Dumper \%count_of; 

    #note - ordered, but includes duplicates. 
    #you could just use keys %count_of, but that would be unordered. 
    foreach my $sequence (grep { $count_of{$_} > 1 } @sequences) { 
     print $sequence, " => ", $count_of{$sequence},"\n"; 
    } 
    print "\n"; 
} 

Nous parcourons re cordon par enregistrement, capturer et supprimer la ligne «en-tête», puis épisser ensemble le reste. Capturez ensuite chaque séquence (chevauchée) de 4 et comptez-les.

Ceci, pour vos données d'échantillon (première strophe par souci de concision):

NC_001422.1 Enterobacteria phage phiX174 sensu lato, complete genome 
Found sequences: 
    GAGT => 2 
    AGTT => 2 
    TTAT => 2 
    CATG => 2 
    ATGA => 3 
    TGAC => 2 
    CGCA => 2 
    AGTT => 2 
    ACTT => 2 
    tttt => 3 
    tttt => 3 
    tttt => 3 
    GGAT => 2 
    GATA => 2 
    ATAT => 2 
    TATT => 2 
    ATGA => 3 
    TGAG => 2 
    GAGT => 2 
    AAAA => 2 
    AAAA => 2 
    ACTT => 2 
    TGAG => 2 
    GGAT => 2 
    GATA => 2 
    tata => 2 
    tata => 2 
    TTAT => 2 
    TATG => 2 
    ATAT => 2 
    TATT => 2 
    GCCG => 2 
    TATG => 2 
    GCCG => 2 
    CGCA => 2 
    CATG => 2 
    ATGA => 3 
    TGAC => 2 

Remarque - parce qu'il est basé sur les séquences d'origine, il est basé sur la commande dans les données, et vous verrez TGAC là deux fois parce que ... c'est là deux fois.

Cependant vous pourriez plutôt:

foreach my $sequence (sort { $count_of{$b} <=> $count_of{$a} } 
          grep { $count_of{$_} > 1 } 
           keys %count_of) { 
     print $sequence, " => ", $count_of{$sequence},"\n"; 
    } 
    print "\n"; 

qui rejettera tout avec moins de 2 matches, et l'ordre par la fréquence.