Recuperando sequências fasta usando informações de arquivo de leito do arquivo instalado localmente

Question

Eu tenho um arquivo .bed contendo cerca de 30.000 linhas para as quais tenho as seqüências recuperadas usando fetch-sequences module da ferramenta rsat ( link ). [Nota: Esta ferramenta conecta-se ao servidor a cada vez para recuperar as seqüências]

Agora, tenho cerca de 10000 subconjuntos do mesmo arquivo de leito classificados aleatoriamente para os quais gostaria de recuperar as sequências. Usando o módulo fetch-sequences , leva 10 horas para conseguir isso. Mas eu gostaria de fazê-lo rapidamente primeiro recuperando as seqüências do arquivo original. Usando este arquivo local, gostaria de recuperar as seqüências para o subconjunto.

Existe uma maneira de fazer isso usando comandos linux ou perl? Não sei como fazer isso com o arquivo instalado localmente. Por favor ajude.

Aqui está o arquivo da minha cama de amostra:

chr10   91061477    91062132    peak4069    41  .   134.220 -1  -1
chr12   77456450    77457116    peak7216    97  .   313.820 -1  -1
chr7    150754939   150755706   peak30140   87  .   281.000 -1  -1
chr3    11643031    11643625    peak20536   135 .   435.740 -1  -1
chr19   6521662 6522869 peak14719   264 .   851.950 -1  -1
chr14   35008667    35009076    peak9034    40  .   131.480 -1  -1
chr6    88851148    88852148    peak27572   55  .   178.560 -1  -1
chr6    20212045    20212627    peak26579   44  .   144.630 -1  -1
chr2    136422022   136422722   peak17485   83  .   270.330 -1  -1
chr11   45951365    45952105    peak4995    284 .   915.840 -1  -1
chr19   50181053    50181876    peak15733   101 .   328.650 -1  -1
chr9    36140202    36140695    peak32014   42  .   137.080 -1  -1
chr4    40992483    40993431    peak23120   40  .   129.060 -1  -1
chr14   50528290    50529818    peak9133    256 .   826.310 -1  -1
chr18   57542948    57543911    peak14298   244 .   789.750 -1  -1
chr21   44528945    44529572    peak19741   45  .   148.260 -1  -1
chr6    16763102    16763743    peak26552   84  .   272.680 -1  -1
chr1    44678710    44679433    peak803 97  .   312.860 -1  -1
chr12   11323475    11324633    peak6358    123 .   396.790 -1  -1
chr9    98271450    98271859    peak32325   37  .   120.160 -1  -1
chr19   2167913 2169475 peak14575   455 .   1470.150    -1  -1
chr16   80613819    80615001    peak12054   261 .   843.100 -1  -1
chr12   118574314   118574830   peak7774    43  .   141.040 -1  -1
chr19   59083917    59085058    peak15917   129 .   418.440 -1  -1
chr2    191184311   191184984   peak17942   103 .   332.220 -1  -1
chr15   85956548    85957254    peak10816   179 .   578.110 -1  -1
chr4    84031272    84032016    peak23411   60  .   195.570 -1  -1
chr12   32012425    32013045    peak6654    45  .   148.210 -1  -1
chr6    70575973    70577060    peak27441   52  .   168.800 -1  -1
chr12   57852728    57853291    peak7023    59  .   192.900 -1  -1
chr10   120879718   120880402   peak4449    209 .   675.760 -1  -1
chr1    28833424    28834023    peak526 35  .   114.020 -1  -1
chr8    60963955    60965013    peak30803   329 .   1062.570    -1  -1
chr7    77326420    77326889    peak29382   32  .   103.320 -1  -1
chr5    133476115   133476527   peak25468   37  .   121.150 -1  -1
chr19   45909117    45910074    peak15572   69  .   222.490 -1  -1
chr5    16467271    16468036    peak24373   117 .   380.290 -1  -1
chr15   66682042    66683502    peak10489   182 .   589.480 -1  -1
chr9    35603806    35604394    peak31993   71  .   229.000 -1  -1
chr4    48249067    48249653    peak23178   50  .   162.070 -1  -1
chr3    112775853   112776466   peak21577   62  .   202.250 -1  -1
chr12   12867020    12867982    peak6428    33  .   106.930 -1  -1
chr6    35655359    35655985    peak27066   53  .   171.010 -1  -1
chr6    74171305    74172116    peak27466   161 .   521.390 -1  -1
chr11   12195905    12196539    peak4741    256 .   826.330 -1  -1
chr2    55386393    55386871    peak16583   40  .   131.810 -1  -1
chr9    37030245    37030957    peak32041   89  .   290.090 -1  -1
chr4    30431566    30431997    peak22948   66  .   214.250 -1  -1
chr10   16612633    16613149    peak3304    49  .   160.900 -1  -1

Aqui está a amostra de minhas seqüências fasta buscadas (para as primeiras 3 linhas no arquivo de exemplo acima):

>hg19_chr10_91061478_91062132_+ range=chr10:91061478-91062132 5'pad=0 3'pad=0 strand=+ repeatMasking=none
AATTGTATTACAAGTCCCCAACTTAATCTTTTCGAATATGAAATAAGAGAGGGACAGTGCACAAGAGCAATGTCCCCAGACCCATCTTTAAGTGAAGCACCAGGCCGATGAAACATCCCTCTCTGCTGCCTTCTTTCTCTGATCACAACTCAGCTCCGGAGGAAAAAGAGTCCTCTAAAGTATAATAAAAAGAAAAAAAGAAAAAGAGTCCTGCCAATTTCACTTTCTAGTTTCACTTTCCCTTTTGTAACGTCAGCTGAAGGGAAACAAACAAAAAGGAACCAGAGGCCACTTGTATATATAGGTCTCTTCAGCATTTATTGGTGGCAGAAGAGGAAGATTTCTGAAGAGTGCAGCTGCCTGAACCGAGCCCTGCCGAACAGCTGAGAATTGCACTGCAACCATGAGGTAAATATTTTCCCTTCGTATTCGGTAGTGCTGTTGAGTCATCTTGTCCAATGCAAATCCTGAGAAGCTATGTTCCCAAAGAGGGCCAGCTCCATTTTAGTGTTTGTTTATAGCCTTACTATGCCTCTACCTCTGTTGGTTGTAAATCTGTCTTACCAATGGTGGTTTGTTCCCTCCTGAACAATTTTCTGCTTCACACTGGCAAGCTTCCTAAATTCATCTCCAGAACTGCACGCCTGGGGAGTTG
>hg19_chr12_77456451_77457116_+ range=chr12:77456451-77457116 5'pad=0 3'pad=0 strand=+ repeatMasking=none
GTCTTTGTTGAAGGTACTTGATAAAAGTCATTGACCCAGAGGAAGAGAAGTAAAGCACTGACTTAGACGTTATATAAATGTATGGATGTGTATTTTTTCAAGGCTGAACCATCCAAATTGGAAAGGAAAACAAAGTTTTGCTCTAAAACTCTCAAAGCCAAAACTCTGAATATATACTTTAAGTCTGGGCATTTCCACCCTCATGACTTAGATAATTAAAAAAAAAAAAAAAAGGCCACTTTAAATAATCTTCACTTTATCTGTGGTTTCACTTTCAGTGGCCAACTGCGGTCCAAAAATATCACATGGAAAATTCCAGAAATAAACAATTCATGAGTTTTAGATTGTGTGCAGTTCTGTGTAATGAAATCTCACGTCATCCTGCTCCGTCCTGCTTCGGATGTGACTCACCCCTTTGTCCAGCGTATTTGCACGGTAGATACTACCTGCTCGAGCAGCCACTGTGTTTTCAGGCTGGCTGTCACGGTATTGCAGTGCTCATGTTCGAGTAACTCTTATTTGACTTCATAATGGCTCCAAAGCACAAGAGTAGTGATGCTGGCAATTTGGATATGCCAAAGGGAAGCCATAAAGTGCTTCTTTTAAGTGAAAAGGTGAACGTTCTTGACTTAAGGAAAGAAAATCGTACGCCAAGGTTGCTAAGAT
>hg19_chr7_150754940_150755706_+    range=chr7:150754940-150755706 5'pad=0 3'pad=0 strand=+ repeatMasking=none
GCGGCCGCGGGGACCCCTGCGGGCCCTCGGTTTTAAGACTCTGGCCCCGGCGTTGCAGAGGAGGCGGCACCTGAGCCTCCAGCCCCGCCCCGTCTGCCCGCGCAGGGCCTTCTGGGGCTTGTAGTCCTAAAACGACTAGGTCTCCCCGCAATGCCCGAGACCGAGGACAAGAAGACTACAACCCCCAGCCGTCTGCGCTCACGCTCTCTGCAGCACTCGACTCCAGGGCTCCGTTTCTACCGCGGAGGCAAACCTTGGACTTCAAGTCCCAGGAAGCAACGCGTGTCCGTTTTCCGGGCGCCCCGCCGCGGCGCCGTGGTTCCCAGTGTGCCCCGCTGTTGTTATTCCTTTTATGTGCTCCCAGCCCTCTTGAAAAGGGCCGCTCCGGGACTACGCGTTCCAGAATGCAGCGGAAATGGGGGCGGAGCGCTCTCGGTTAGGGGTTTGGGGTTTGGCGGCCTAGATCCCGGGCACTGGCGGCCCAGCGCTGACCTGGTTGGTGGCATTGTGTTCCCAACGGCCTCTTGACGACCTCAGCACGGGTTTCCACCTCTCCCCAAGCCACCTAGTGACCCCAGAATTGACTGGGGAATGCCTGTGAGCGATGATGACCTCACAGGGAACAGCTGACCGCAGGGCTGGGAGAACAGCTGTGCCCCTTCGAGGCTGGATTTTAGTGGAGGGACACACGCCAAAGACCCCCTCTCTGCTGAGCCCCGTTTGTTGTCTCGGAGCCCACCCGACTCTAGCCGCTGAACTCTGACATG

    

Answer 1

Isso deve fazer o que você precisa:

for i in $(awk '{print $1"_"$2+1"_"$3}' foo.bed); do grep -A 1 $i seqs.fa ; done

Explicação

• awk '{print $1"_"$2+1"_"$3}' foo.bed : isso imprimirá o cromossomo e as coordenadas de início e fim de cada linha do arquivo .bed. Observe o $2+1 , pois as coordenadas no seu arquivo são diferentes das do .bed . Por exemplo, para as três primeiras linhas do seu arquivo:

  $ awk '{print $1"_"$2+1"_"$3}' foo.bed 
  chr10_91061478_91062132
  chr12_77456451_77457116
  chr7_150754940_150755706
  

• for i in $(command); do ...; done : isso salvará cada valor retornado por command como $i e, em seguida, fará algo com ele.

• grep -A 1 $i seqs.fa : esse é o "algo". Ele irá grep para cada resultado do comando awk e imprimirá a linha correspondente, bem como a próxima ( -A 1 ).

Se você precisar fazer isso novamente, não use o RSAT! ¹ Há maneiras muito mais simples de fazer isso. Em vez disso, baixe as seqüências FASTA de cada cromossomo e, em seguida, use as ferramentas do pacote exonerate (instalável em sistemas baseados no Debian com sudo apt-get install exonerate ). O procedimento (supondo que você tenha um arquivo FASTA para cada cromossomo chamado chrN.fa ) seria:

awk '{print $1,$2,($3-$2)}' foo.bed | while read chr start length; do
    fastasubseq /path/to/$chr.fa $start $length >> subseqs.fa
done

O comando acima irá extrair as subseqüências de interesse (assumindo que as coordenadas estão sempre com relação à faixa mais como deveriam) e deve levar apenas alguns segundos.

¹ _{Não me entenda mal, a suíte RSAT é uma ótima ferramenta e eu tenho muito respeito por seu autor com quem trabalhei por alguns anos, é apenas não necessariamente a melhor ferramenta para trabalhos de grande escala.}