У меня есть файл последовательности белка, я хочу посчитать тримеры в нем, используя sed или grep

У меня есть файл последовательности белка в следующем формате

uniprotID\space\sequence

последовательность - строка любой длины, но только с 20 разрешенными буквами, т.е.

ARNDCQEGHILKMFPSTWYV

Пример 1 записи

Q5768D AKCCACAKCCAC

Я хочу создать CSV-файл в следующем формате

Q5768D     

12
ACA 1
AKC 2
CAC 2
CAK 1
CCA 2
KCC 2

Вот что я сейчас пытаюсь:

#!/bin/sh
while read ID SEQ # uniprot along with sequences
do
echo $SEQ | tr -d '[[:space:]]' | sed 's/./& /g'  > TEST_FILE
declare -a SSA=(`cat TEST_FILE`)
SQL=$(echo ${#SSA[@]})
  for (( X=0; X <= "$SQL"; X++ ))
      do
         Y=$(expr $X + 1)
         Z=$(expr $X + 2)
         echo ${SSA[X]} ${SSA[Y]} ${SSA[Z]}
     done  | awk '{if (NF == 3) print}' | tr -d ' ' > TEMPTRIMER
rm TEST_FILE # removing temporary sequence file
sort TEMPTRIMER|uniq -c > $ID.$SQL
done < $1

в этом коде я храню отдельные записи в другом файле, что не хорошо. Кроме того, программа работает очень медленно, за 12 часов из 12 миллионов записей получают доступ только к 12000 записей.

Источник

2 ответа

Решение

Этот Perl-скрипт обрабатывает приблизительно 550'000 "триммеров"/ сек. (случайные действительные тестовые последовательности длиной 0–8000 символов, записи по 100 КБ (~400 МБ) создают выходной CSV объемом 2 ГБ)

выход:

Q1024A;421;AAF=1;AAK=1;AFC=1;AFE=2;AGP=1;AHC=1;AHE=1;AIV=1;AKN=1;AMC=1;AQD=1;AQY=1;...
Q1074F;6753;AAA=1;AAD=1;AAE=1;AAF=2;AAN=2;AAP=2;AAT=1;ACA=1;ACC=1;ACD=1;ACE=3;ACF=2;...

код:

#!/usr/bin/perl
use strict;
$|=1;
my $c;

# process each line on input
while (readline STDIN) {
  $c++; chomp;
  # is it a valid line? has the format and a sequence to process
  if (m~^(\w+)\s+([ARNDCQEGHILKMFPSTWYV]+)\r?$~ and $2) {
    print join ";",($1,length($2));
    my %trimdb;
    my $seq=$2;
    #split the sequence into chars
    my @a=split //,$seq;
    my @trimmer;

    # while there are unprocessed chars in the sequence...
    while (scalar @a) {

      # fill up the buffer with a char from the top of the sequence
      push @trimmer, shift @a;

      # if the buffer is full (has 3 chars), increase the trimer frequency
      if (scalar @trimmer == 3 ) {
        $trimdb{(join "",@trimmer)}++;
        # drop the first letter from buffer, for next loop
        shift @trimmer;
      }
    }

    # we're done with the sequence - print the sorted list of trimers
    foreach (sort keys %trimdb) {

      #print in a csv (;) line
      print ";$_=$trimdb{$_}";
    }
    print"\n";
  }
  else {
    #the input line was not valid.
    print STDERR "input error: $_\n";
  }
  # just a progress counter
  printf STDERR "%8i\r",$c if not $c%100;
}
print STDERR "\n";

если у вас установлен Perl (большинство Linux-систем это делает, проверьте путь /usr/bin/perl или замените его своим), просто запустите: ./count_trimers.pl < your_input_file.txt > output.csv

Источник

Если это то, что вы хотите:

$ cat file
Q5768D AKCCACAKCCAC
OTHER FOOBARFOOBAR
$
$ awk -f tst.awk file
Q5768D  OTHER
12      12
AKC 2   FOO 2
KCC 2   OOB 2
CCA 2   OBA 2
CAC 2   BAR 2
ACA 1   ARF 1
CAK 1   RFO 1

Это сделает это:

$ cat tst.awk
BEGIN { OFS="\t" }
{
    colNr = NR
    rowNr = 0
    name[colNr] = $1
    lgth[colNr] = length($2)
    delete name2nr
    for (i=1;i<=(length($2)-2);i++) {
        trimer = substr($2,i,3)
        if ( !(trimer in name2nr) ) {
            name2nr[trimer] = ++rowNr
            nr2name[colNr,rowNr] = trimer
        }
        cnt[colNr,name2nr[trimer]]++
    }
    numCols = colNr
    numRows = (rowNr > numRows ? rowNr : numRows)
}
END {
    for (colNr=1;colNr<=numCols;colNr++) {
        printf "%s%s", name[colNr], (colNr<numCols?OFS:ORS)
    }
    for (colNr=1;colNr<=numCols;colNr++) {
        printf "%s%s", lgth[colNr], (colNr<numCols?OFS:ORS)
    }
    for (rowNr=1;rowNr<=numRows;rowNr++) {
        for (colNr=1;colNr<=numCols;colNr++) {
            printf "%s %s%s", nr2name[colNr,rowNr], cnt[colNr,rowNr], (colNr<numCols?OFS:ORS)
        }
    }
}

Если вместо этого вам нужен вывод, как в perl-ответе @rogerovo, это было бы намного проще, чем выше, и более эффективно и использовало бы намного меньше памяти:

$ cat tst2.awk
{
    delete cnt
    for (i=1;i<=(length($2)-2);i++) {
        cnt[substr($2,i,3)]++
    }
    printf "%s;%s", $1, length($2)
    for (trimer in cnt) {
        printf ";%s=%s", trimer, cnt[trimer]
    }
    print ""
}

$ awk -f tst2.awk file
Q5768D;12;ACA=1;KCC=2;CAK=1;CAC=2;CCA=2;AKC=2
OTHER;12;RFO=1;FOO=2;OBA=2;OOB=2;ARF=1;BAR=2
Источник
Другие вопросы по тегам