Perl-программа для поиска k-mer с определенной последовательностью

Я пытаюсь улучшить программу perl, которую я ранее написал, чтобы она распознала верхнюю 1000-ю длину 23 k-mers, которая заканчивается GG, и распечатывать k-mers, которые только появляется один раз в последовательности. Однако независимо от того, где я добавляю reg exp, я не могу получить ожидаемый результат.Perl-программа для поиска k-mer с определенной последовательностью

код у меня есть:

#!/usr/bin/perl 
use strict; 
use warnings; 

my $k   = 23; 
my $input  = 'Fasta.fasta'; 
my $output  = 'Fasta2.fasta'; 
my $match_count = 0; 

#Open File 
unless (open(FASTA, "<", $input)) { 
    die "Unable to open fasta file", $!; 
} 

#Unwraps the FASTA format file 
$/ = ">"; 

#Separate header and sequence 
#Remove spaces 
unless (open(OUTPUT, ">", $output)) { 
    die "Unable to open file", $!; 
} 

<FASTA>; # discard 'first' 'empty' record 

my %seen; 
while (my $line = <FASTA>) { 
    chomp $line; 
    my ($header, @seq) = split(/\n/, $line); 
    my $sequence = join '', @seq; 

    for (length($sequence) >= $k) { 
     $sequence =~ m/([ACTG]{21}[G]{2})/g; 

     for my $i (0 .. length($sequence) - $k) { 
      my $kmer = substr($sequence, $i, $k); 

      ##while ($kmer =~ m/([ACTG]{21}[G]{2})/g){ 
      $match_count = $match_count + 1; 
      print OUTPUT ">crispr_$match_count", "\n", "$kmer", "\n" unless $seen{$kmer}++; 
     } 
    } 
}

Входной файл FASTA выглядит следующим образом:

> >2L type=chromosome_arm; loc=2L:1..23011544; ID=2L; dbxref=REFSEQ:NT_033779,GB:AE014134; MD5=bfdfb99d39fa5174dae1e2ecd8a231cd; length=23011544; release=r5.54; species=Dmel; 
CGACAATGCACGACAGAGGAAGCAGAACAGATATTTAGATTGCCTCTCAT 
TTTCTCTCCCATATTATAGGGAGAAATATGATCGCGTATGCGAGAGTAGT 
GCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATGGTGGCGGA 
TGAACGAGATGATAATATATTCAAGTTGCCGCTAATCAGAAATAAATTCA 
TTGCAACGTTAAATACAGCACAATATATGATCGCGTATGCGAGAGTAGTG 
CCAACATATTGTGCTAATGAGTGCCTCTCGTTCTCTGTCTTATATTACCG 
CAAACCCAAAAAGACAATACACGACAGAGAGAGAGAGCAGCGGAGATATT 
TAGATTGCCTATTAAATATGATCGCGTATGCGAGAGTAGTGCCAACATAT 
TGTGCTCTCTATATAATGACTGCCTCTCATTCTGTCTTATTTTACCGCAA 
ACCCAAATCGACAATGCACGACAGAGGAAGCAGAACAGATATTTAGATTG 
CCTCTCATTTTCTCTCCCATATTATAGGGAGAAATATGATCGCGTATGCG 
AGAGTAGTGCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATG 
GTGGCGGATGAACGAGATGATAATATATTCAAGTTGCCGCTAATCAGAAA 
TAAATTCATTGCAACGTTAAATACAGCACAATATATGATCGCGTATGCGA 
GAGTAGTGCCAACATATTGTGCTAATGAGTGCCTCTCGTTCTCTGTCTTA 
TATTACCGCAAACCCAAAAAGACAATACACGACAGAGAGAGAGAGCAGCG 
GAGATATTTAGATTGCCTATTAAATATGATCGCGTATGCGAGAGTAGTGC 
CAACATATTGTGCTCTCTATATAATGACTGCCTCTCATTCTGTCTTATTT 
TACCGCAAACCCAAATCGACAATGCACGACAGAGGAAGCAGAACAGATAT

и так далее ...

ожидаемый результат (распечатать 23k- с окончанием GG, появляются только после в последовательности) Я надеюсь получить:

>crispr_1 
GGGTGGAGCTCCCGAAATGCAGG 
>crispr_2 
TTAATAAATATTGACACAGCGGG 
>crispr_3 
ATCGTGGGGCGTTTTGTGAAAGG 
>crispr_4 
AGTTTTTCACATAATCAGACAGG 
>crispr_5 
GTGTTGGATGAGTGTCCTCTGGG 
>crispr_6 
ATAGGTTGGTTGTTTTAAAAGGG 
>crispr_7 
AAATTTTTGTTGCCACTGAATGG 
>crispr_8 
AAGTTTCGAACTACGATGGTTGG 
>crispr_9 
CATGCTTTGTGGAAATAAGTCGG 
>crispr_10 
CACAGTGGGTGTTTGCACCTCGG 
.... and so on

Текущий код, который я сделал создать файл Fasta со следующим:

>crispr_1 
CGACAATGCACGACAGAGGAAGC 
>crispr_2 
GACAATGCACGACAGAGGAAGCA 
>crispr_3 
ACAATGCACGACAGAGGAAGCAG 
>crispr_4 
CAATGCACGACAGAGGAAGCAGA 
>crispr_5 
AATGCACGACAGAGGAAGCAGAA 
>crispr_6 
ATGCACGACAGAGGAAGCAGAAC 
>crispr_7 
TGCACGACAGAGGAAGCAGAACA 
>crispr_8 
GCACGACAGAGGAAGCAGAACAG 
>crispr_9 
CACGACAGAGGAAGCAGAACAGA 
>crispr_10 
ACGACAGAGGAAGCAGAACAGAT 
.... and so on

, а если я удалить

for (length($sequence) >=$k){ 
$sequence =~m/([ACTG]{21}[G]{2})/g;

и добавьте ## в то время как ($ kmer = ~ м/([АКТГ] {21} [G] {2})/г) {

while ($kmer =~ m/([ACTG]{21}[G]{2})/g){

Я получаю FASTA файл (с результатами, которые не пронумерованы правильно и unabl е различать дублированные и уникальные последовательности):

>crispr_1 
CATTTTCTCTCCCATATTATAGG 
>crispr_2 
ATTTTCTCTCCCATATTATAGGG 
>crispr_3 
TATTGTGCTCTTTGATTTTTTGG 
>crispr_4 
GATTTTTTGGCAACCCAAAATGG 
>crispr_5 
TTTTTGGCAACCCAAAATGGTGG 
>crispr_6 
TTGGCAACCCAAAATGGTGGCGG 
>crispr_7 
ACGACAGAGAGAGAGAGCAGCGG 
>crispr_8 
AAATCGACAATGCACGACAGAGG 
>crispr_91 
TATTGTGATCTTCGATTTTTTGG 
>crispr_93 
TTTTTGGCAACCCAAAATGGAGG 
.... and so on

Я попытался переместить регулярное выражение вокруг моего кода, но ни один из них не генерируются ожидаемый результат. Я не знаю, что я сделал не так. Я не добавляю выход программы, когда счетчик достигает 1000 в код.

Заранее благодарен!

источник

2016-11-04 Kaiser

Можете ли вы предоставить нам (вырезать при необходимости) ввод вместе с выводом, который вы хотите от этого входа. Очень сложно помочь, когда вы не видите данные, с которыми имеете дело. –

@DaveCross Я предоставил входной файл; это просто общий файл fasta, содержащий генетическую последовательность. Благодарю. – Kaiser

Я в замешательстве: первая ожидаемая строка (одна для 'crispr_1'):' CATTTTCTCTCCCATATTATAGG', но во входном файле, который вы показываете, нет последовательной последовательности. Как вы пришли к этой конкретной последовательности? –

Я не уверен, что полностью понимаю ваш вопрос, но может быть следующее, что вам нужно.

<FASTA>; # discard 'first' 'empty' record 

my %data; 
while (my $line = <FASTA>){ 
    chomp $line; 
    my($header, @seq) = split(/\n/, $line); 
    my $sequence = join '', @seq; 

    for my $i (0 .. length($sequence) - $k) { 
     my $kmer = substr($sequence, $i, $k); 

     $data{$kmer}++ if $kmer =~ /GG$/; 
    } 
} 
my $i = 0; 
for my $kmer (sort {$data{$b} <=> $data{$a}} keys %data) { 
    printf "crispr_%d\n%s appears %d times\n", ++$i, $kmer, $data{$kmer}; 
    last if $i == 1000; 
}

Некоторые выход на файл у меня есть,:

crispr_1 
ggttttccggcacccgggcctgg appears 4 times 
crispr_2 
ccgagctgggcgagaagtagggg appears 4 times 
crispr_3 
gccgagctgggcgagaagtaggg appears 4 times 
crispr_4 
gcacccgggcctgggtggcaggg appears 4 times 
crispr_5 
agcagcgggatcgggttttccgg appears 4 times 
crispr_6 
gctgggcgagaagtaggggaggg appears 4 times 
crispr_7 
cccttctgcttcagtgtgaaagg appears 4 times 
crispr_8 
gtggcagggaagaatgtgccggg appears 4 times 
crispr_9 
gatcgggttttccggcacccggg appears 4 times 
crispr_10 
tgagggaaagtgctgctgctggg appears 4 times 
crispr_11 
agctgggcgagaagtaggggagg appears 4 times 

. . . . 

ggcacccgggcctgggtggcagg appears 4 times 
crispr_50 
gaatctctttactgcctggctgg appears 4 times 
crispr_51 
accacaacattgacagttggtgg appears 2 times 
crispr_52 
caacattgacagttggtggaggg appears 2 times 
crispr_53 
catgctcatcgtatctgtgttgg appears 2 times 
crispr_54 
gattaatgaagtggttattttgg appears 2 times 
crispr_55 
gaaaccacaacattgacagttgg appears 2 times 
crispr_56 
aacattgacagttggtggagggg appears 2 times 
crispr_57 
gacttgatcgattaatgaagtgg appears 2 times 
crispr_58 
acaacattgacagttggtggagg appears 2 times 
crispr_59 
gaaccatatattgttatcactgg appears 2 times 
crispr_60 
ccacagcgcccacttcaaggtgg appears 1 times 
crispr_61 
ctgctcctgggtgtgagcagagg appears 1 times 
crispr_62 
ccatatattatctgtggtttcgg appears 1 times 

. . . .

Update Чтобы получить результаты, которые вы сослались в своем комментарии (ниже), замените выходной код:

my $i = 1; 

while (my ($kmer, $count) = each %data) { 
    next unless $count == 1; 
    print "crispr_$i\n$kmer\n"; 
    last if $i++ == 1000; 
}

Чтобы ответить, мой собственный комментарий, чтобы получить first 1000.

<FASTA>; # discard 'first' 'empty' record 

my %seen; 
my @kmers; 
while (my $line = <FASTA>){ 
    chomp $line; 
    my($header, @seq) = split(/\n/, $line); 
    my $sequence = join '', @seq; 

    for my $i (0 .. length($sequence) - $k) { 
     my $kmer = substr($sequence, $i, $k); 

     if ($kmer =~ /GG$/) { 
      push @kmers, $kmer unless $seen{$kmer}++; 
     } 
    } 
} 

my $i = 1; 
for my $kmer (@kmers) { 
    next unless $seen{$kmer} == 1; 
    print "crispr_$i\n$kmer\n"; 
    last if $i++ == 1000; 
}

Ответ Для проверки уникальности последних 12 символов, заканчивающихся в GG, код ниже достигает этого.

 if ($kmer =~ /(.{10}GG)$/) { 
      my $substr = $1; 
      push @kmers, $kmer unless $seen{$substr}++; 
     } 

my $i = 1; 
for my $kmer (@kmers) { 
    my $substr = substr $kmer, -12; 
    next unless $seen{$substr} == 1; 
    print "crispr_$i\n$kmer\n"; 
    last if $i++ == 1000; 
}

источник

2016-11-04 22:46:38

То, что я пытаюсь достичь, - это распечатать первые 1000 23-километровых последовательностей с помощью GG, который появляется только один раз. Благодаря! – Kaiser

@Sunny Я обдумал это, и мое решение не будет печатать * first * kmers, заканчивающиеся на GG. Он печатает любой из сортировщиков, которые появляются в файле в случайном порядке. Если вам нужно * first * 1000, тогда потребуется другое решение. –

Я проверил результат, полученный из вашего обновленного кода, а головка -100 и хвост -100 совпадают с ожидаемым результатом вывода. И это, похоже, не производит в случайном порядке.Поэтому я не думаю, что с этим решением что-то не так. Большое спасибо за проверку! – Kaiser

На самом деле, эта строка кода

$sequence =~m/([ACTG]{21}[G]{2})/g;

эта линия только на матч регулярного выражения, если вы пытаетесь напечатать эту $sequence, он, несомненно, распечатать первоначальный результат.

пожалуйста, добавьте код segement как этот

if($sequence =~/([ACTG]{21}[G]{2}$)/g) 
{ 


}#please remember to match the end with $.

КСТАТИ Похоже, мультипликатор для цикла разбора эти данные не очень разумно, скорость синтаксического анализа без лучшей эффективности.

источник

2016-11-04 03:55:37

Эрик, эта методология не работает, она распечатает пустой файл fasta. – Kaiser

Perl-программа для поиска k-mer с определенной последовательностью

ответ

Смежные вопросы