ajuda a entender uma fórmula do awk que descompacta arquivos fasta

3

Acabei de encontrar uma fórmula que pode ser usada para descompactar arquivos fasta. Antes de dar a fórmula, preciso explicar o que é desdobrar um arquivo fasta. Em suma, o formato fasta é assim:

>name_of_sequence$
xxxxxxxxxxxxxxxxxxxxxx$
>name_of_sequence_2$
xxxxxxxxxxxxxxxxxxxxxx$
>name_of_sequence_3$
xxxxxxxxxxxxxxxxxxxxxx$

Este seria um arquivo fasta normal, pois eu tenho apenas uma linha por sequência (xxxxxx ...). O cifrão são quebras de linha.

Por vezes, no entanto, pode encontrar ficheiros fasta wrapped como este:

>name_of_sequence$
xxxxxxxxx$
xxxxxxxxx$
xxxx$
>name_of_sequence_2$
xxxxxxxxx$
xxxxxxxxx$
xxxx$
>name_of_sequence_3$
xxxxxxxxx$
xxxxxxxxx$
xxxx$

Aqui, você ainda tem apenas três sequências, mas cada uma delas está dividida em três partes. Desempacotar um arquivo fasta significa converter o último formato para o primeiro (uma linha por sequência).

Para fazer isso, você precisa remover quebras de linha do último arquivo, mas não todas elas. Você precisaria manter a quebra de linha após o nome da sequência (por exemplo, aqui: > nome_de_sequência $) e no final da sequência (por exemplo, aqui: xxxx $).

Parece que essa fórmula faz isso:

cat infasta | awk '/^>/{print s? s"\n"$0:$0;s="";next}{s=s sprintf("%s",$0)}END{if(s)print s}' > outfasta

Minha pergunta é: alguém poderia me explicar como funciona?

    
por Agathe 19.02.2017 / 18:58

2 respostas

5

Este é o seu script awk :

/^>/ {
    print s ? s "\n" $0 : $0;
    s = "";
    next;
}

{
    s = s sprintf("%s", $0);
}

END {
    if (s)
      print s;
}

O primeiro bloco é acionado apenas para linhas que começam com > , ou seja, linhas de cabeçalho fasta.

No primeiro bloco, algo é impresso. Esse algo é s ? s "\n" $0 : $0 . Isso significa que "se s for diferente de zero (ou não definida), use s e adicione uma nova linha a ela seguida por toda a linha atual, caso contrário, use apenas toda a linha atual". Neste programa, s será uma sequência parcialmente lida pertencente à linha de cabeçalho processada mais recentemente, e quando o programa atingir uma linha de cabeçalho, essa instrução print exibirá a última sequência (que agora está completa), se houver foi qualquer, seguido pela linha de cabeçalho recém-encontrada em uma nova linha.

O bloco então define s para uma string vazia (ainda não lemos nenhuma sequência pertencente a este cabeçalho), e passamos para a próxima linha de entrada.

O próximo bloco é executado para todas as linhas de entrada (mas não para linhas de cabeçalho, pois elas serão ignoradas devido ao next no bloco anterior). Ele simplesmente anexa a linha atual a s . sprintf é usado, mas não sei bem por que ( s = s $0 provavelmente funcionaria também).

O último bloco será executado depois de ler todas as linhas de entrada. Ele imprimirá a seqüência que pertence à última linha de cabeçalho, se houver alguma.

Resumo:

O script awk concatena todas as linhas de sequência separadas, salvando-as em uma variável. Quando uma linha de cabeçalho é encontrada, ela exibe a sequência lida até o momento junto com o novo cabeçalho em uma linha própria. No final, a sequência pertencente ao último cabeçalho é exibida.

Script alternativo awk que não armazena sequência em uma variável (pode ser útil se você tiver genomas muito grandes em seus arquivos fasta):

/^>/ {
    if (NR == 1) {
        print;  # 1st header line, just print it.
    } else {
        # Print a newline for the prev. sequence, then the header line on its own line.
        printf("\n%s\n", $0);
    }
    next; # Skip to next input line.
}

{
    printf("%s", $0); # Print sequence without newline.
}

END {
    printf("\n"); # Add final newline to output.
}

Como "one-liner":

awk '/^>/{if(NR==1){print}else{printf("\n%s\n",$0)}next} {printf("%s",$0)} END{printf("\n")}' sequence.fasta
    
por 19.02.2017 / 19:19
0

O FWIW forneceu uma solução baseada em "sed" para arquivos fasta envoltos. O fluxo subjacente do método sed é encontrar a linha do nome da sequência, em primeiro lugar mostramos essa linha em uma linha por si só e então começamos a acumular no próprio espaço padrão as linhas de seqüência & também tirando a nova linha enquanto nós vamos. Esse fluxo quebra quando atingimos a próxima linha de nome de sequência ou eof.

sed -e '
  /^>/{                  # caught sequence name line
     n                   # print seq name, next line into pattern space
     :loop
        N                # read next line into PS, if not print PS/quit
        /\n>/!s/\n//     # join successive sequences
     /\n/!bloop          # go back for more seq if new seq name not got yet
     P;D                 # print the current seq then delete it, branch to the top with PS having new seq name
  }
' your_fasta_file
    
por 19.02.2017 / 20:31