shell script para executar o protocolo tophat

Question

Então, estou usando um protocolo chamado 'smoking' para análise de dados de sequenciamento de RNA. É mais uma espécie de questão técnica relacionada ao script de shell. Eu posso fazer isso na linha de comando e não tenho problemas como tal. Como estou fazendo isso em um cluster, gostaria de usar um script que possa automatizar minha tarefa.

Assim, o comando do protocolo é assim:

1. tophat
2. abotoaduras
3. cuffmerge

4. cuffdiff

   O primeiro comando faz todo o alinhamento que gera algum arquivo que eu tenho que usar para o próximo comando cufflinks , então cuffmerge e finalmente cuffdiff .

Alguém pode me ajudar a escrever um script de shell simples, que pode chamar cada um desses comandos e fazer a tarefa.

Qualquer ajuda seria muito apreciada.

argumentos

tophat -p 8 -G genes.gtf -o C1_R1_thout genome C1_R1_1.fq C1_R1_2.fq
cufflinks -p 8 -o C1_R1_clout C1_R1_thout/accepted_hits.bam
cuffmerge -g genes.gtf -s genome.fa -p 8 assemblies.txt 
cuffdiff -o diff_out -b genome.fa -p 8 –L C1,C2 -u merged_asm/merged.gtf \
./C1_R1_thout/accepted_hits.bam,./C1_R2_thout/accepted_hits.bam,./C1_R3_thout/accepted_hits.bam \
./C2_R1_thout/accepted_hits.bam,./C2_R3_thout/accepted_hits.bam,./C2_R2_thout/accepted_hits.bam

onde, 'p' corresponde ao número do processador o '-o' corresponde ao diretório de saída e repousa o '-g' corresponde ao arquivo de anotação que estou usando para anotar minhas leituras RAW que serão alinhadas.

    

Answer 1

A solução simples e frágil

Vamos escrever um script simples chamado hailmary.sh

#!/bin/bash
#The first line should always be just as it is above
#This script is called hailmary.sh
#because we run this script and we need to pray
#that all four commands will run correctly
#If one of them fail, you may not get the results

tophat -p 8 -G genes.gtf -o C1_R1_thout genome C1_R1_1.fq C1_R1_2.fq
cufflinks -p 8 -o C1_R1_clout C1_R1_thout/accepted_hits.bam
cuffmerge -g genes.gtf -s genome.fa -p 8 assemblies.txt 
cuffdiff -o diff_out -b genome.fa -p 8 –L C1,C2 -u merged_asm/merged.gtf ./C1_R1_thout/accepted_hits.bam,./C1_R2_thout/accepted_hits.bam,./C1_R3_thout/accepted_hits.bam ./C2_R1_thout/accepted_hits.bam,./C2_R3_thout/accepted_hits.bam,./C2_R2_thout/accepted_hits.bam

1. Copie e cole todas as linhas acima, incluindo as linhas que começam com "#" no gedit e salve como hailmary.sh.

2. No Nautilus, clique com o botão direito no arquivo que você acabou de criar e selecione %código%. Vá para a guia Properties e coloque uma marca de seleção ao lado Permitir executar o arquivo como programa .

   Como alternativa, em um terminal, digite:

   chmod + x hailmary.sh

3. Para executar o script em um terminal, digite:

   ./hailmary.sh

O Permissions antes do nome é necessário e assume que o arquivo está no local atual do diretório. Se você colocar o arquivo em uma pasta que esteja no caminho, como ./ , não precisará do /home/<userid>/bin . Se você colocar em outro lugar, precisará escrever o caminho inteiro, como:

/home/<userid>/myspecialfolder/hailmary.sh

Note que os quatro comandos e seus argumentos estão em quatro linhas separadas. Se você quiser colocá-los em uma única linha, você deve separá-los por ./ ou && . Não há necessidade de ; se estiverem em linhas separadas.

Em qualquer um desses casos, o segundo comando não será iniciado até que o primeiro seja concluído (ou falhe!).

O problema com essa abordagem é que ela não verifica se o primeiro comando foi executado com êxito antes de executar o segundo, e assim por diante. Portanto, se ; falhar por algum motivo, o script continuará com a sequência de cufflink, cuffmerge e cuffdiff. É por isso que eu chamo esse script de tophat .

Fonte: link

Script com verificação de saída de tophat

#!/bin/bash
#The first line should always be just as it is above
#This script is called hailmary2.sh
#This script runs tophat
#then checks for the existance of the output file
#If the output is found, it runs the rest

tophat -p 8 -G genes.gtf -o C1_R1_thout genome C1_R1_1.fq C1_R1_2.fq

if [[ -f "./C1_R1_thout/accepted_hits.bam" ]]; then
    echo "tophat finished. Proceeding with the rest"
    cufflinks -p 8 -o C1_R1_clout C1_R1_thout/accepted_hits.bam
    cuffmerge -g genes.gtf -s genome.fa -p 8 assemblies.txt 
    cuffdiff -o diff_out -b genome.fa -p 8 –L C1,C2 -u merged_asm/merged.gtf ./C1_R1_thout/accepted_hits.bam,./C1_R2_thout/accepted_hits.bam,./#C1_R3_thout/accepted_hits.bam ./C2_R1_thout/accepted_hits.bam,./C2_R3_thout/accepted_hits.bam,./C2_R2_thout/accepted_hits.bamfi
else echo "tophat did not complete"
fi

Espero que isso ajude até que alguém forneça uma resposta mais elegante.