La conversión de FASTQ a FASTA con SED / AWK

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

  •  20-09-2019
  •  | 
  •  

Pregunta

Tengo un dato en que siempre viene en bloque de cuatro en el siguiente formato (llamado 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/>/

¿Hay una manera sencilla de sed / awk / bash para convertirlos en este formato (llamado FASTA): 

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

En principio, queremos extraer las dos primeras líneas de cada bloque-de-4 y reemplazar @ con >.

¿Fue útil?

Solución

Esta es una vieja cuestión, y ha habido muchas soluciones diferentes que se ofrecen. Dado que la respuesta aceptada utiliza sed, pero tiene un problema evidente (que es que va a sustituir a @ con> cuando el signo @ aparece como la primera letra de la línea de calidad), me siento obligado a ofrecer una solución a base de sed simple que realmente funciona : 

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

El único supuesto hecho es que cada lectura ocupa exactamente 4 líneas en el archivo FASTQ, pero que parece bastante seguro, en mi experiencia.

El script fastq_to_fasta en el kit de herramientas FASTX también funciona. (Vale la pena mencionar que es necesario especificar la opción -Q33 para dar cabida a las codificaciones + 33 Qual ahora comunes Phred. Lo que es divertido, ya que es tirar los datos de calidad de todos modos!)

Otros consejos

sed no está muerto. Si estamos golf: 

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

O, emulando http: //www.ringtail.tsl .ac.uk / David-Studholme / scripts / fastq2fasta.pl publicado por Pierre, que sólo imprime la primera palabra (el identificador) de la primera línea y lo hace (a) el tratamiento de errores: 

#!/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 que hay un montón de herramientas existentes para convertir estos formatos; que es mejor usar estos en lugar de cualquier información incluida aquí (incluyendo los anteriores).

Como se detalla en Cock, et al (2009) NAR, muchas de estas soluciones son incorrectas, ya que "el carácter '@' marcador (ASCII 64) puede ocurrir en cualquier parte de la cadena de calidad. Esto significa que cualquier programa de análisis no debe tratar una línea que comienza con '@' como indica el inicio del siguiente registro, sin comprobar, además, la longitud de la cadena de calidad hasta la fecha coincide con la longitud de la secuencia. "

http://ukpmc.ac.uk/articlerender.cgi?accid=PMC2847217 para más detalles.

acaba de awk, sin necesidad de otras herramientas 

# 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

Me gustaría escribir 

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

Este es el más rápido que tengo, y me lo metió en mi archivo .bashrc: 

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

no falla en las líneas poco frecuentes pero no imposibles de calidad que comienzan con @ ... pero que fallan en FASTQ envuelto, si eso es incluso legal (aunque existe).

Aquí está la solución a la parte de "saltar cada dos líneas" del problema que acabo de aprender de SO: 

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

Si todo lo que hay que hacer es cambiar una @ a >, entonces calculo 

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

    read line_to_skip
    read line_to_skip
done

hará el trabajo.

Algo así como: 

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

debería funcionar.

Creo que, con grep de GNU esto se podría hacer con esto: 

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

abajo 

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

donde los datos es el archivo de datos. He recibido: 

>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/>/

Sé que estoy lejos en el futuro, pero en beneficio de los empleados de Google:

Es posible que desee utilizar fastq_to_fasta de la caja de herramientas FASTX . Asimismo, mantendrá el signo @, sin embargo. Además, eliminará las líneas con Ns menos que le diga que no.

Usted podría estar interesado en bioawk, es una versión adaptada del awk que está sintonizado para procesar archivos FASTA 

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

Nota: BioAwk se basa en awk de Brian Kernighan que se documenta en el "El lenguaje de programación AWK", por Al Aho, Brian Kernighan, y Peter Weinberger (Addison-Wesley, 1988, ISBN 0-201-07981-X) . No estoy seguro de si esta versión es compatible con POSIX .

Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top