EDITAR e uma solução
Como minha pergunta original foi mal formulada e eu estava tentando reinventar a roda, estou respondendo minha própria pergunta agora (talvez ajude outra pessoa):
gff2fasta é uma ferramenta que faz exatamente o que eu preciso, que é extrair um determinado pedaço de DNA de um genoma completo (um arquivo enorme em formato fasta chamado FULLGENOME.fasta).
Se eu souber onde está a peça que quero, posso fazer um arquivo chamado TEST.gff onde especifico o andaime (aqui: sca_5_chr5_3_0) e o início (aqui: 2390621) e fim (aqui: 2391041) da minha peça para gff2fasta, por exemplo:
sca_5_chr5_3_0 JGI CDS 2390621 2391041 . + 0 name "e_gw1.5.88.1"; proteinId 40463;
Então eu só preciso correr
gff2fasta.pl -gff TEST.gff -fasta FULLGENOME.fasta -out test
gff2fasta irá então obter as informações de TEST.gff e gerar o bit que desejo no arquivo chamado test, que ficará assim:
>sca_5_chr5_3_0_2390621_2391041_F
ATGACGACCCGCGCACCCAAAGACACATACGCTCAGCCCGACTATGAGGAGGCTCACCTAGCGACGTTTGCAGCCCCAAA
AGGCTACCCTATCGAGTCTATGCTACCCCCTAGCGTGAAGAGGGAGACCTTTGAACAGGCCCTAGCCGAGTTTACCGACG
CTATCGGCAAAGATTATGTCTTTATCGGCGATGCTCTCTCTCACTACATCGATCCCTACGACATCTATATCGATGATAGT
GAGAGAAGGAGGATGCCGAGCGCGGCTGTTTGGTACGTAACGCGGCACCCATCGAAACCTAGCAGCACTGACAAGTTTCC
GCGTCTAGTCCCTCTTCGCTCGAGGAAGCTAAGCAGGCTCTCAAGGTTGCGAACAAATACGGCATTCCGATCTGGGCATT
TTCCAGAGGCAAGAACCTGGG
Obrigado a terdon e a todos os outros pela ajuda!
-------------
Esta é a continuação desta pergunta, com mais informações e detalhesunix: obtém caracteres de 10 a 80 em um arquivo
Acho que estou quase lá, mas ainda preciso de ajuda.
Tentei explicar da forma mais clara possível, mas estou ciente de que é uma questão muito especializada, por isso, informe-me se puder esclarecer mais alguma coisa!
O que tenho são três arquivos:
nota para o administrador: é possível fazer upload de arquivos? Não tenho ideia de como...
Umarquivo(N_haematococca_1.fasta) do qual preciso extrair um nome:
head -1 N_haematococca_1.fasta | cut -f4 -d'|'
este nome neste caso é:
e_gw1.5.88.1
Problema 1: O código acima funciona bem, mas tenho problemas para salvar o nome (e_gw1.5.88.1) em uma variável que posso usar mais tarde...
Quero salvar esse nome em uma variável, vamos chamá-la:
firstline
Asegundoarquivo (Necha2_best_models.gff) onde desejo todas as linhas onde este nome ocorre:
grep -ir "e_gw1.5.88.1" --include="Necha2*.gff" > Necha2_in_genome.list
mas com a variável nomeada:
grep -ir $firstline --include="Necha2*.gff" > Necha2_in_genome.list
Isso funciona para o uso de "e_gw1.5.88.1"...
O arquivo que estou criando aqui me diz o nome do fragmento de DNA que quero cortar (sca_5_chr5_3_0) e qual parte dele eu preciso (do caractere número 2390621 ao caractere número 2392655). Preciso dessa informação para obter as partes certas do terceiro arquivo
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))
Em umterceiroarquivo do qual desejo cortar certas partes após uma palavra-chave (sca_5_chr5_3_0) neste caso:
Graças a Kamaraj e hschou, tenho uma solução parcial para isso agora:
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
No entanto, se eu depurar isso com números pequenos:
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
Eu recebo esta saída:
r1_1_0
CCTTATCCTAGCG
nmapped
CTTATATATTAT
nmapped
TAAAAGGAGTTA
unmapped
TCTTATATAAA
unmapped
AATCTTAAGAA
Parece que a opção RS é ignorada e imprime apenas os caracteres de 10 a 20 para determinadas linhas. Não tenho ideia de por que essas linhas foram selecionadas.
sca_5_chr5_3_0
ocorre apenas uma vez no arquivo.
outros nomes que existem são
>sca_66_unmapped
>sca_67_unmapped
etc.
Tenho que obter essas informações de 178 genomas, todos são arquivos enormes e pesquisar manualmente não é uma opção.
COMO ficam os arquivos:
N_haematococca_1.fasta (arquivo 1) é um arquivo fasta normal:
>jgi|Necha2|40463|e_gw1.5.88.1
MTTRAPKDTYAQPDYEEAHLATFAAPKGYPIESMLPPSVKRETFEQALAEFTDAIGKDYVFIGDALSHYI
DPYDIYIDDSERRRMPSAAVCPSSLEEAKQALKVANKYGIPIWAFSRGKNLGYGGPSARVNGSVAFDLHR
Necha2_best_models.gff (arquivo 2) tem esta aparência (apenas mais longa):
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 (arquivo 3) se parece um pouco com isto (bits GATC muito mais longos...):
>sca_8_chr1_1_0
CCTTATCCTAGCGAGGATTAAGGTCCTCGAAAAGAAAGGCAAAGATATAAAGGTAATATAATAGAAGATT
AAGGTATTCTAAGTAAAGGTTATAAAGAAATAAAATAAGAAGAATATTTATAGGCTAAGAAAGACCCCCC
TAAAGGTTAAGGACTTAATATTAAGATTTAATATTCCCTAATTAATTAATATATTAATAAAAATAAAGAT
>sca_5_chr5_3_0
ATGACCACTATATCCATCGGCACAACGGCGTTAATCACATTTGGGTCTGCAATTTTCTGTTTTTGCGGTT
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAGACTCGCTCTGCTTT
>sca_67_unmapped
CTACAGGAAGCCCTAGAGGCCCTGCAAGACCTTCCAAAGCAGCACTCTTTGCTTCTTCTTAGAGGCAGTA
TCCAGCTTCTTCTAAGGCACCTCCAGAGGCAGCTAGACCCAACCGGGCTAAAAGACCTTTGGGAAGAGGC
TAATACCCTTATAAGAGAGGCTATTATAGCCCTAGTGGCTAGAAGCCCTAGTGAGAGCCCTAAAGAGCCT
AATTCAAGCCTTATAGCCCTTCCAGTCAGAGAGGGAGGCCTAGGAATACCCTTACACAAGGACCTAGCCC
Resultado final esperado: é um único pedaço de texto após a palavra-chave >sca_5_chr5_3_0
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAG