Как подсчитать количество вхождений комбинаций символов n-длины в строке

Question

Я использую следующий один лайнер, чтобы перечислять вхождения комбинаций ATCG, образуя цепочку длиной 6. Он отлично работает в стороне от того, чтобы не печатать событие из 0 матчей. Есть ли способ изменить регулярное выражение или другую часть на то, где он будет печатать что-то вроде «0 ATTTAG»?

#!/bin/bash 
for file in e_coli.fa 
do 
    base=$(basename $file .fa) 
    cat $file | perl -nE 'say for /(?<=([ATCG]{6}))/g' \ 
     | sort | uniq -c >> ${base}_hexhits_6mer.txt 
done 

stdout: 
    465 AAAAAA 
    607 AAAAAC 
    661 AAAAAG 
    581 AAAAAT 
    563 AAAACA 
    807 AAAACC 
    770 AAAACG 
    373 AAAACT 
    663 AAAAGA 
    1213 AAAAGC 

Answer 1

С uniq -c подсчитывает количество раз происходит линия, она не может вернуться 0. Требуемое изменение требует полного переписывания.

perl -e' 
    while (<>) { 
     ++$counts{$_} for /(?=([ATCG]{6}))/g; 
    } 

    for my $seq (glob("{A,C,G,T}" x 6)) { 
     printf("%7d %s\n", $counts{$seq}, $seq); 
    } 
' "$file" >"${base}_hexhits_6mer.txt" 

Answer 2

То, что вы просите, намного сложнее. Чтобы узнать, что вы не видели, вам необходимо знать все возможные комбинации символов, которые затем можно фильтровать.

Здесь я использую подход скользящего окна с помощью substr в Perl, чтобы найти число всех «видели» ATCG символов в строке As, Ts, Cs и Gs (как прочитанный из __DATA__) в хэш. Затем они сортируются, так что наиболее часто наблюдаемый 6-мер сначала отображается и распечатывается.

use strict; 
use warnings; 

my @bases = qw/ A G C T /; 

my %data; 

for my $a1(@bases){ 
    for my $a2(@bases){ 
     for my $a3(@bases){ 
      for my $a4(@bases){ 
       for my $a5(@bases){ 
        for my $a6(@bases){ 
         $data{"$a1$a2$a3$a4$a5$a6"} = 0; 
        } 
       } 
      } 
     } 
    } 
} 

my $nucs = <DATA>; 

my $len = length($nucs); 

for (my $i = 0; $i <= $len - 6; $i++) { 
    my $kmer = substr($nucs, $i, 6); 
    next if $kmer =~ tr/ACGT//c; 
    $data{$kmer}++; # populate hash with "seen" 6-mers 
} 

# print out sorted hash 
foreach my $seq (sort { $data{$b} <=> $data{$a} } keys %data){ 
    print "$seq,$data{$seq}\n"; 
} 

__DATA__ 
ATGCCCGTCGTAGTCATGCATGCATCGATCGATGCATGCTACGTGTTGT 

Там явно будет лучше/красивее способ вычисления всех перестановок символов в строке, чем то, что я сделал, но это работает.

Как говорит Бородин, это печатает главным образом вариации «невидимых» строк.

Answer 3

Самый простой способ построить хэш подсчетов встречаемости для каждого шаблона, а затем напечатать количество всех возможных моделей

Эта программа использует glob трюк, чтобы создать список всех возможных строк шесть-символов, сформированных из A, T, C и G

use strict; 
use warnings 'all'; 

my @files = qw/ e_coli.fa /; 

my %counts; 

for my $file (@files) { 

    open my $fh, '<', $file or die qq{Unable to open "$file" for input: $!}; 

    while (<$fh>) { 
     ++$counts{$1} while /(?= ([ATCG]{6})) /gx; 
    } 
} 

for my $pattern (glob '{A,T,C,G}' x 6) { 
    printf "%4d %s\n", $counts{$pattern} // 0, $pattern; 
} 

Answer 4

в случае у вас есть много данных, и вам нужно что-то немного быстрее, вот C раствор:

#include <stdio.h> 
#include <stdlib.h> 
#include <string.h> 

void reader(FILE* in, unsigned long hist[4096]) { 
    for (unsigned long key=0, count=0;;) { 
    switch(getc(in)) { 
     case EOF:      return; 
     case 'A': key <<= 2;   break; 
     case 'C': key <<= 2; key += 1; break; 
     case 'G': key <<= 2; key += 2; break; 
     case 'T': key <<= 2; key += 3; break; 
     default: count=0;    continue; 
    } 
    if (count == 5) ++hist[key & 0xFFF]; 
    else ++count; 
    } 
} 

int putkey(FILE* out, unsigned long key) { 
    char s[6]; 
    for (int j=6; j--; key >>= 2) s[j] = "ACGT"[key&3]; 
    return fprintf(out, "%.6s", s); 
} 

void writer(FILE* out, unsigned long hist[4096]) { 
    for (unsigned long key = 0; key < 4096; ++key) { 
    fprintf(stdout, "%7lu ", hist[key]); 
    putkey(out, key); 
    putchar('\n'); 
    } 
} 

int main(int argc, char** argv) { 
    FILE* in = stdin; 
    if (argc > 1) in = fopen(argv[1], "r"); 
    if (!in) { perror(argv[1]); exit(1); } 
    unsigned long hist[4096] = {0}; 
    reader(in, hist); 
    writer(stdout, hist); 
    return 0; 
} 

Потребовалось чуть меньше половины секунды, чтобы обработать образец fastq 31 МБ (который, как это бывает, включает в себя все 4096 возможных шестисимвольных последовательностей); решения Perl занимали 12 секунд (фугу) и 18 секунд (икегами/бородин) соответственно.