如何统计字符串中n个长度的字符组合出现的次数

How to count the number of occurences of n-length combinations of characters in string

我正在使用下面的一行来列出 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

由于uniq -c计算一行出现的次数,所以不可能return0。请求的更改需要完全重写。

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"

您要执行的操作要复杂得多。要了解您没有看到的内容,您首先需要了解所有可能的字符组合,然后您可以根据这些组合进行筛选。

在这里,我使用 Perl 中的 substr 的滑动 window 方法来查找 As 字符串中的所有 "seen" ATCG 个字符,散列中的 TsCsGs(从 __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

显然会有一种 better/prettier 方法来计算字符串中字符的所有排列,而不是我所做的,但它确实有效。

正如 Borodin 所说,这主要打印出 "unseen" 字符串的变体。

最简单的方法是为每个模式构建一个出现次数的散列,然后打印所有可能模式的次数

此程序使用 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{} while /(?= ( [ATCG]{6} ) ) /gx;
    }
}

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

如果您有大量数据并且需要更快的速度,这里有一个 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;
}

处理一个 31MB 的 fastq 样本(碰巧包括所有 4096 个可能的六字符序列)只花了不到半秒的时间; Perl 解决方案分别用了 12 秒 (fugu) 和 18 秒 (ikegami/borodin)。