contando múltiplos padrões em uma única passagem com grep?

Question

Eu escrevi um loop de grep para contar os trinucleotídeos de DNA em um arquivo fasta de DNA compactado contendo seqüências de DNA, por exemplo

declare -a tri=(AAA AAC AAG AAT CAA .. etc)

for i in ${tri[@]}
do
   gzip -cd gencode.v18.pc_transcripts.fa.gz | grep -v "^>" | grep -o $i | wc -l
done

Onde o arquivo fasta está nesse formato (embora muito maior)

head test.fa
>id1
TTTTTAAAAA
>id2
GGGGGCCCCC
etc..

Embora isto funcione (isto é, contagens ocorrências de cada trinucleótido) é, na minha opinião, bastante ineficiente, uma vez que tem de passar através dos dados 64 vezes (uma vez para cada trinucleótido possível).

Minha pergunta é como usar bash ou grep existe uma maneira de contar cada trinucleotídeo em uma única passagem pelo arquivo (como os arquivos são bem grandes)?

thx

    

Answer 1

IFS=$'\n'
gzip -dc file.gz | grep -v '^>' | grep -Foe "${tri[*]}" | sort | uniq -c

Mas, a propósito, AAAC corresponde a AAA e AAC , mas grep -o produzirá apenas uma delas. É isso que você quer? Além disso, quantas ocorrências de AAA em AAAAAA ? 2 ou 4 ( [AAA]AAA , A[AAA]AA , AA[AAA]A , AAA[AAA] )?

Talvez você queira:

gzip -dc file.gz | grep -v '^>' | fold -w3 | grep -Fxe "${tri[*]}" | sort | uniq -c

Isso é dividir as linhas em grupos de 3 caracteres e contar as ocorrências como linhas inteiras (encontraria 0 ocorrência de AAA em ACAAATTCG (como ACA AAT TCG )).

Ou, por outro lado:

gzip -dc file.gz | awk '
  BEGIN{n=ARGC;ARGC=0}
  !/^>/ {l = length - 2; for (i = 1; i <= l; i++) a[substr($0,i,3)]++}
  END{for (i=1;i<n;i++) printf "%s: %d\n", ARGV[i], a[ARGV[i]]}' "${tri[@]}"

(encontraria 4 ocorrências de AAA em AAAAAA ).

    

Answer 2

@ O segundo exemplo do stéphane-chazelas é ótimo, mas o comando sort pode realmente desacelerar quando os dados ficarem maiores.

Podemos supor que os caracteres em linhas sem cabeçalho são nucleotídeos válidos? Isso removeria a tri correspondente.

gzip -dc file.gz | grep -v '^>' | fold -w3 | awk '{a[$0]++} END{for(codon in a) {printf "%s: %d\n", codon, a[codon]}'

Em uma nota relacionada, converter os tokens de 3-nucleotídeos em octal de 2 bytes (ou binários de lá) reduziria o tamanho do fluxo de dados em 2/3 ou 1/3.

gzip -dc file.gz | grep -v '^>' | tr ACGT 0123 | fold -w3 | cat <(echo "obase=8; ibase=4;" ) - | bc | xargs printf "%02d\n" | tee foo.octal | xxd -r -p > foo.bin

De volta aos nucleotídeos

cat foo.bin | xxd -p | fold -w2 | cat <(echo "obase=4; ibase=8;") - | bc | xargs printf "%03d" | tr 0123 ACGT

xargs é provavelmente o passo limitante aqui, talvez analisar N linhas por vez iria corrigir isso ou usar gnu paralelo.