EDITAR y una solución
Debido a que mi pregunta original estaba mal redactada y estaba tratando de reinventar la rueda, ahora estoy respondiendo mi propia pregunta (tal vez ayude a alguien más):
gff2fasta es una herramienta que hace exactamente lo que necesito, que es extraer un fragmento determinado de ADN de un genoma completo (un archivo enorme en formato fasta llamado FULLGENOME.fasta).
Si sé dónde está la pieza que quiero, puedo crear un archivo llamado TEST.gff donde especifico el andamio (aquí: sca_5_chr5_3_0) y el principio (aquí: 2390621) y el final (aquí: 2391041) de mi pieza para gff2fasta, Por ejemplo:
sca_5_chr5_3_0 JGI CDS 2390621 2391041 . + 0 name "e_gw1.5.88.1"; proteinId 40463;
Entonces solo necesito correr
gff2fasta.pl -gff TEST.gff -fasta FULLGENOME.fasta -out test
gff2fasta luego obtendrá la información de TEST.gff y generará el bit que quiero en el archivo llamado prueba, que se verá así:
>sca_5_chr5_3_0_2390621_2391041_F
ATGACGACCCGCGCACCCAAAGACACATACGCTCAGCCCGACTATGAGGAGGCTCACCTAGCGACGTTTGCAGCCCCAAA
AGGCTACCCTATCGAGTCTATGCTACCCCCTAGCGTGAAGAGGGAGACCTTTGAACAGGCCCTAGCCGAGTTTACCGACG
CTATCGGCAAAGATTATGTCTTTATCGGCGATGCTCTCTCTCACTACATCGATCCCTACGACATCTATATCGATGATAGT
GAGAGAAGGAGGATGCCGAGCGCGGCTGTTTGGTACGTAACGCGGCACCCATCGAAACCTAGCAGCACTGACAAGTTTCC
GCGTCTAGTCCCTCTTCGCTCGAGGAAGCTAAGCAGGCTCTCAAGGTTGCGAACAAATACGGCATTCCGATCTGGGCATT
TTCCAGAGGCAAGAACCTGGG
¡Gracias a terdon y a todos los demás por la ayuda!
-------------
Esta es la continuación de esta pregunta, con más información y detalles.Unix: obtiene caracteres del 10 al 80 en un archivo
Creo que ya casi he llegado, pero todavía necesito ayuda.
He intentado explicarlo lo más claramente posible, pero soy consciente de que es una pregunta muy especializada, así que ¡avíseme si puedo aclarar algo más!
Lo que tengo son tres archivos:
nota para el administrador: ¿es posible cargar archivos? No tengo ni idea de cómo...
Unoarchivo(N_haematococca_1.fasta) del cual necesito extraer un nombre:
head -1 N_haematococca_1.fasta | cut -f4 -d'|'
este nombre en este caso es:
e_gw1.5.88.1
Problema 1: El código anterior funciona bien pero tengo problemas para guardar el nombre (e_gw1.5.88.1) en una variable que puedo usar para más adelante...
Quiero guardar ese nombre en una variable, llamémosla:
firstline
Asegundoarchivo (Necha2_best_models.gff) donde quiero todas las líneas donde aparece este nombre:
grep -ir "e_gw1.5.88.1" --include="Necha2*.gff" > Necha2_in_genome.list
pero con la variable nombrada:
grep -ir $firstline --include="Necha2*.gff" > Necha2_in_genome.list
Esto funciona para el uso de "e_gw1.5.88.1"...
El archivo que estoy creando aquí me dice el nombre del fragmento de ADN que quiero cortar (sca_5_chr5_3_0) y qué parte necesito (desde el carácter número 2390621 hasta el carácter número 2392655). Necesito esa información para obtener los bits correctos del tercer archivo.
a=(`awk -F '\t' '{print $4}' Necha2_in_genome.list`)
startDNA=${a[1]}
endDNA=${a[${#a[@]}-1]}
#add 1000 or other number, depending on the problems with the gene
correctedstartDNA=$(($startDNA-1000))
correctedendDNA=$(($endDNA-1000))
en unterceroarchivo del cual quiero cortar ciertas partes después de una palabra clave (sca_5_chr5_3_0) en este caso:
Gracias a Kamaraj y hschou, ahora tengo una solución parcial para esto:
cat Necha2_scaffolds.fasta | sed -n -e '/sca_5_chr5_3_0/,$p' | grep -v '^>' | tr -d '\n'|awk -v start="${correctedstartDNA}" -v end="${correctedendDNA}" '{print substr($0,start,end)}' RS= Necha2_scaffolds.fasta
Sin embargo, si depuro esto con números pequeños:
cat Necha2_scaffolds.fasta | sed -n -e '/sca_5_chr5_3_0/,$p' | grep -v '^>' | tr -d '\n'|awk -v start=10 -v end=20 '{print substr($0,start,end)}' RS= Necha2_scaffolds.fasta
Obtengo este resultado:
r1_1_0
CCTTATCCTAGCG
nmapped
CTTATATATTAT
nmapped
TAAAAGGAGTTA
unmapped
TCTTATATAAA
unmapped
AATCTTAAGAA
Parece que la opción RS se ignora y solo imprime los caracteres del 10 al 20 para determinadas líneas. No tengo idea de por qué se seleccionan estas líneas.
sca_5_chr5_3_0
solo ocurre una vez en el archivo.
otros nombres que hay son
>sca_66_unmapped
>sca_67_unmapped
etc.
Tengo que obtener esta información de 178 genomas, todos son archivos enormes y buscar manualmente no es una opción.
CÓMO se ven los archivos:
N_haematococca_1.fasta (archivo 1) es un archivo fasta normal:
>jgi|Necha2|40463|e_gw1.5.88.1
MTTRAPKDTYAQPDYEEAHLATFAAPKGYPIESMLPPSVKRETFEQALAEFTDAIGKDYVFIGDALSHYI
DPYDIYIDDSERRRMPSAAVCPSSLEEAKQALKVANKYGIPIWAFSRGKNLGYGGPSARVNGSVAFDLHR
Necha2_best_models.gff (archivo 2) se ve así (un poco más largo):
sca_5_chr5_3_0 JGI exon 2390621 2390892 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2390621 2390892 . + 0 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 1
sca_5_chr5_3_0 JGI start_codon 2390621 2390623 . + 0 name "e_gw1.5.88.1"
sca_5_chr5_3_0 JGI exon 2390949 2391041 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2390949 2391041 . + 2 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 2
sca_5_chr5_3_0 JGI exon 2391104 2391311 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2391104 2391311 . + 2 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 3
sca_5_chr5_3_0 JGI exon 2391380 2392367 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2391380 2392367 . + 0 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 4
sca_5_chr5_3_0 JGI exon 2392421 2392485 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2392421 2392485 . + 1 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 5
sca_5_chr5_3_0 JGI exon 2392541 2392657 . + . name "e_gw1.5.88.1"; transcriptId 40463
sca_5_chr5_3_0 JGI CDS 2392541 2392657 . + 0 name "e_gw1.5.88.1"; proteinId 40463; exonNumber 6
sca_5_chr5_3_0 JGI stop_codon 2392655 2392657 . + 0 name "e_gw1.5.88.1"
sca_5_chr5_3_0 JGI exon 2396205 2396730 . + . name "e_gw1.5.227.1"; transcriptId 41333
sca_5_chr5_3_0 JGI CDS 2396205 2396730 . + 0 name "e_gw1.5.227.1"; proteinId 41333; exonNumber 1
sca_5_chr5_3_0 JGI start_codon 2396205 2396207 . + 0 name "e_gw1.5.227.1"
Necha2_scaffolds.fasta (archivo 3) se parece un poco a esto (bits GATC mucho más largos...):
>sca_8_chr1_1_0
CCTTATCCTAGCGAGGATTAAGGTCCTCGAAAAGAAAGGCAAAGATATAAAGGTAATATAATAGAAGATT
AAGGTATTCTAAGTAAAGGTTATAAAGAAATAAAATAAGAAGAATATTTATAGGCTAAGAAAGACCCCCC
TAAAGGTTAAGGACTTAATATTAAGATTTAATATTCCCTAATTAATTAATATATTAATAAAAATAAAGAT
>sca_5_chr5_3_0
ATGACCACTATATCCATCGGCACAACGGCGTTAATCACATTTGGGTCTGCAATTTTCTGTTTTTGCGGTT
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAGACTCGCTCTGCTTT
>sca_67_unmapped
CTACAGGAAGCCCTAGAGGCCCTGCAAGACCTTCCAAAGCAGCACTCTTTGCTTCTTCTTAGAGGCAGTA
TCCAGCTTCTTCTAAGGCACCTCCAGAGGCAGCTAGACCCAACCGGGCTAAAAGACCTTTGGGAAGAGGC
TAATACCCTTATAAGAGAGGCTATTATAGCCCTAGTGGCTAGAAGCCCTAGTGAGAGCCCTAAAGAGCCT
AATTCAAGCCTTATAGCCCTTCCAGTCAGAGAGGGAGGCCTAGGAATACCCTTACACAAGGACCTAGCCC
Producto final esperado: es un solo fragmento de texto después de la palabra clave >sca_5_chr5_3_0
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAG