Convertendo o fastq em fasta com sed/awk

StackOverflow https://stackoverflow.com/questions/1542306

  •  20-09-2019
  •  | 
  •  

Pergunta

Eu tenho dados que sempre vem no bloco de quatro no formato a seguir (chamado fastq):

@SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
+SRR018006.2016 GA2:6:1:20:650 length=36
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+!
@SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+SRR018006.19405469 GA2:6:100:1793:611 length=36
7);;).;);;/;*.2>/@@7;@77<..;)58)5/>/

Existe uma maneira simples de sed/awk/bash de convertê -los neste formato (chamado fasta):

>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

Em princípio, queremos extrair as duas primeiras linhas em cada bloco de 4 e substituir @ com >.

Foi útil?

Solução

Esta é uma pergunta antiga, e houve muitas soluções diferentes oferecidas. Como a resposta aceita usa sed, mas tem um problema gritante (que é que ele substituirá @ por> :

sed -n '1~4s/^@/>/p;2~4p'

A única suposição feita é que cada leitura ocupa exatamente 4 linhas no arquivo FastQ, mas isso parece bastante seguro, na minha experiência.

O script fastQ_TO_FASTA no kit de ferramentas Fastx também funciona. (Vale a pena mencionar que você precisa especificar a opção -q33 para acomodar as codificações agora comuns Phred+33 Qual. O que é engraçado, já que está jogando fora os dados de qualidade de qualquer maneira!)

Outras dicas

sed não está morto. Se estamos jogando golfe:

sed '/^@/!d;s//>/;N'

Ou, emulando http://www.ringtail.tsl.ac.uk/david-studholme/scripts/fastq2fastha.pl Postado por Pierre, que apenas imprime a primeira palavra (o ID) da primeira linha e faz (alguns) erros de manuseio:

#!/usr/bin/sed -f
# Read a total of four lines
$b error
N;$b error
N;$b error
N
# Parse the lines
/^@\(\([^ ]*\).*\)\(\n[ACGTN]*\)\n+\1\n.*$/{
  # Output id and sequence for FASTA format.
  s//>\2\3/
  b
}
:error
i\
Error parsing input:
q

Parece haver muitas ferramentas existentes para converter esses formatos; Você provavelmente deve usá -los em vez de qualquer coisa postada aqui (incluindo o acima).

Conforme detalhado em Cock, et al (2009) NAR, muitas dessas soluções estão incorretas, pois "o caractere de marcador"@"(ASCII 64) pode ocorrer em qualquer lugar da sequência de qualidade. Isso significa que qualquer analisador não deve tratar uma linha começando por '@' como indicando o início do próximo registro, sem verificar adicionalmente o comprimento da sequência de qualidade até agora corresponde ao comprimento da sequência. "

Ver http://ukpmc.ac.uk/articlerender.cgi?accid=pmc2847217 para detalhes.

Apenas estranho, não precisa de outras ferramentas

# awk '/^@SR/{gsub(/^@/,">",$1);print;getline;print}' file
>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

Eu escrevia

awk '
    NR%4 == 1 {print ">" substr($0, 2)}
    NR%4 == 2 {print}
' fastq > fasta

Este é o mais rápido que eu tenho, e eu o prendi no meu arquivo .bashrc:

alias fq2fa="awk '{print \">\" substr(\$0,2);getline;print;getline;getline}'"

Ele não falha nas linhas de qualidade pouco frequentes, mas não impossíveis, que começam com @... mas falham no fastq embrulhado, se isso é legal (porém, existe).

Aqui está a solução para a parte "pular todas as outras linhas" do problema que acabei de aprender com:

while read line
do
    # print two lines
    echo "$line"
    read line_to_print
    echo "$line_to_print"

    # and skip two lines
    read line_to_skip
    read line_to_skip
done

Se tudo isso precisa ser feito é mudar um @ para >, então eu acho

while read line
do
    echo "$line" | sed 's/@/>/'
    read line
    echo "$line"

    read line_to_skip
    read line_to_skip
done

vai fazer o trabalho.

Algo como:

awk 'BEGIN{a=0}{if(a==1){print;a=0}}/^@/{print;a=1}' myFastqFile | sed 's/^@/>/'

Deveria trabalhar.

Eu acho que, com GNU Grep, isso pode ser feito com isso:

grep -A 1 "^@" t.txt | grep -v "^--" | sed -e "s/^@/\>/"

awk 'BEGIN{P=1}{if(P==1||P==2){gsub(/^[@]/,">");print}; if(P==4)P=0; P++}' data

>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

abaixo de

awk '{gsub(/^[@]/,">"); print}' data

onde dados são o seu arquivo de dados. Eu recebi:

>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
+SRR018006.2016 GA2:6:1:20:650 length=36
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+!
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+SRR018006.19405469 GA2:6:100:1793:611 length=36
7);;).;);;/;*.2>/@@7;@77<..;)58)5/>/

Eu sei que estou no futuro, mas para o benefício dos Googlers:

Você pode querer usar fastq_to_fasta do kit de ferramentas Fastx. Ele manterá o sinal @, no entanto. Ele também removerá linhas com ns, a menos que você diga não.

Você pode estar interessado no Bioawk, é uma versão adaptada do AWK, que é sintonizada para processar arquivos de fasta

bioawk -c fastx '{ print ">"$name ORS $seq }' file.fastq

Observação: Bioawk é baseado em Brian Kernighan's Awk que está documentado em "The Awk Programming Language", de Al Aho, Brian Kernighan e Peter Weinberger (Addison-Wesley, 1988, ISBN 0-201-07981-X). Não tenho certeza se esta versão é compatível com Posix.

Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top