use gff2fasta em vez de um script bash para obter partes de sequências de DNA de um genoma completo

2

EDIT 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 alguém):

O gff2fasta é uma ferramenta que faz exatamente o que eu preciso, que é extrair um dado pedaço de DNA de um genoma completo (um arquivo enorme no formato fasta chamado FULLGENOME.fasta).

Se eu sei onde a peça que eu quero é, eu posso fazer um arquivo chamado TEST.gff onde eu especifico o scaffold (aqui: sca_5_chr5_3_0) e o começo (aqui: 2390621) e end (aqui: 2391041) do meu pedaço 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 mostrar o bit que eu quero 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 pela ajuda!

-------------

Esta é a continuação desta questão, com mais informações e detalhes unix: pega os caracteres de 10 a 80 em um arquivo

Acho que estou quase lá, mas ainda preciso de ajuda.

Eu tentei explicar o mais claramente possível, mas estou ciente de que é uma pergunta muito especializada, então, por favor, deixe-me saber se eu posso esclarecer algo mais!

O que eu tenho são três arquivos:

nota para admin: é possível fazer upload de arquivos? Eu não tenho ideia de como ...

Um arquivo (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 eu tenho problemas para obter o nome (e_gw1.5.88.1) salvo em uma variável que eu possa usar para mais tarde ...

Eu quero salvar esse nome em uma variável, vamos chamá-lo:

 firstline

Um arquivo segundo (Necha2_best_models.gff) onde eu quero 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 eu quero cortar (sca_5_chr5_3_0) e qual bit dele eu preciso (do número de caractere 2390621 para o número de caractere 2392655). Preciso dessa informação para obter os bits certos 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 um terceiro arquivo do qual eu quero cortar certas partes depois de uma palavra-chave (sca_5_chr5_3_0) neste caso:

Graças a Kamaraj e hschou eu 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 10 a 20 para determinadas linhas. Eu não tenho ideia de por que essas linhas são selecionadas.

 sca_5_chr5_3_0

ocorre apenas uma vez no arquivo.

outros nomes que existem

 >sca_66_unmapped 
 >sca_67_unmapped

etc

Eu tenho que obter essa informação de 178 genomas, todos eles são arquivos enormes e pesquisar manualmente não é uma opção.

COMO os arquivos parecem:

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) se parece com isso (só mais):

 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) parece um pouco com isso (muito mais bits GATC ...):

 >sca_8_chr1_1_0
 CCTTATCCTAGCGAGGATTAAGGTCCTCGAAAAGAAAGGCAAAGATATAAAGGTAATATAATAGAAGATT
 AAGGTATTCTAAGTAAAGGTTATAAAGAAATAAAATAAGAAGAATATTTATAGGCTAAGAAAGACCCCCC
 TAAAGGTTAAGGACTTAATATTAAGATTTAATATTCCCTAATTAATTAATATATTAATAAAAATAAAGAT
 >sca_5_chr5_3_0
 ATGACCACTATATCCATCGGCACAACGGCGTTAATCACATTTGGGTCTGCAATTTTCTGTTTTTGCGGTT
 TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
 TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
 CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAGACTCGCTCTGCTTT
 >sca_67_unmapped
 CTACAGGAAGCCCTAGAGGCCCTGCAAGACCTTCCAAAGCAGCACTCTTTGCTTCTTCTTAGAGGCAGTA
 TCCAGCTTCTTCTAAGGCACCTCCAGAGGCAGCTAGACCCAACCGGGCTAAAAGACCTTTGGGAAGAGGC
 TAATACCCTTATAAGAGAGGCTATTATAGCCCTAGTGGCTAGAAGCCCTAGTGAGAGCCCTAAAGAGCCT
 AATTCAAGCCTTATAGCCCTTCCAGTCAGAGAGGGAGGCCTAGGAATACCCTTACACAAGGACCTAGCCC

Saída final esperada : é um único bit de texto após a palavra-chave > sca_5_chr5_3_0

 TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
 TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
 CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAG
    
por gugy 07.04.2017 / 10:22

0 respostas