Configurando o Maui, alternativa ao pbs_sched

O Torque é distribuído com um Scheduler muito simples, o pbs_sched. O Maui é um scheduler muito mais potente, feito pelos mesmos desenvolvedores do Torque (eles também tem uma solução mais completa chamada Moab que é paga).

Para instalar o Maui é preciso primeiro baixar o código fonte nesse endereço (você vai precisar criar uma conta no site andes de baixar o programa). Em seguida descompacte o arquivo em algum diretório, por exemplo em /usr/src. Dentro do diretório do código fonte digite (como root):


./configure
make
make install


Em seguida adicione o diretório /usr/local/maui/bin no path, no debian/ubuntu basta editar o arquivo /etc/environment para fazer isso.

Antes de iniciar o Maui é preciso finalizar o pbs_sched e remove-lo da sequencia de boot, para isso (ainda como root):


/etc/init.d/pbs_sched stop
update_rc.d pbs_sched disable


Edite o arquivo de configuração do maui em /usr/local/maui:

1) Troque:

RMCFG[TANGERINE] TYPE=PBS@RMNMHOST@

Por

RMCFG[base] TYPE=PBS

Adicione a linha
ADMIN3 all

E adicione o seu usuário como ADMIN1, por exemplo:

ADMIN1 root varuzza

Depois de editar o arquivo é preciso configurar a inicialização do maui no boot. Como o maui não vem com um script de iniciação pronto coloque esse arquivo no diretório /etc/init.d que eu criei a partir dos scripts do pbs.


Finalmente digite


/etc/init.d/maui start


Se todo ocorrer corretamente ao digitar o comando showq você deverá obter um output similar à esse:

ACTIVE JOBS--------------------
JOBNAME            USERNAME      STATE  PROC   REMAINING            STARTTIME


     0 Active Jobs       0 of    2 Processors Active (0.00%)
                         0 of    1 Nodes Active      (0.00%)

IDLE JOBS----------------------
JOBNAME            USERNAME      STATE  PROC     WCLIMIT            QUEUETIME


0 Idle Jobs

BLOCKED JOBS----------------
JOBNAME            USERNAME      STATE  PROC     WCLIMIT            QUEUETIME


Total Jobs: 0   Active Jobs: 0   Idle Jobs: 0   Blocked Jobs: 0

Convertendo de FastQ para Fasta

Conveter de FastQ para csfasta é fácil, basta utilizar uma linha de shell script:

zcat SRR034220.fastq.gz
       | awk 'BEGIN{a=0}{if(a==1){print;a=0}}/^@/{print;a=1}'
       | sed 's/^@/>/' | gzip > SRR034220.csfasta.gz


(Fonte: StackOverflow)

Já a conversão de FastQ para Qual é mais complicada porque é preciso converter os valores de qualidade do formato texto para o formato numérico, e para isso é melhor utilizar uma linguagem de programação. Python é uma linguagem de programação elegante e simples de programar, a biblioteca BioPython é bastante sofisticada e tem uma função pronta para fazer essa conversão. O pequeno trecho de código abaixo faz a conversão que queremos:


from Bio import SeqIO
print "Converting files..."
count = SeqIO.convert("SRR034220.fastq", "fastq", "SRR034220.py.qual", "qual")
print "Converted %i records" % count


Extremamente simples, porém, leeento. O script gastou 3 horas de CPU em sua execução. Isso não é muito prático.

Procurando por uma alternativa encontrei a bibiloteca BioHaskell. Haskell é uma linguagem funcional compilada fortemente tipada extremamente elegante e eficiente, com performance próxima do C. O código em haskell é quase tão compacto quanto o código em python:


module Main
where

import Bio.Sequence
import System

main = do
  [inputFile,outputQual] <- getArgs 
  seqs <- (readFastQ inputFile) 
  writeQual outputQual seqs

Porém o tempo de execução do programa foi de 34 minutos, mais de 5 vezes mais rápido do que a versão em Python.

Performance é algo relativo, muitos programas passam mais tempo esperando do efetivamente usando a CPU, como user interfaces por exemplo, mas nesse caso a diferença de performance é importante, pois a versão em haskell pode ser rodada durante o almoço ou uma reunião, enquanto que a versão em python vai gastar toda a manhã para rodar.

Instalando o PBS (torque)

O PBS (Portable Batch System) é um sistema de execução batch de jobs desenvolvido originalmente pela NASA. Existem três versões do PBS: PBS Pro, openPBS e o Torque. Esqueça as outras versões, a unica versão gratuita que funciona é o Torque (Também não use a versão que do Debian, o empacotamento está cheio de bugs).

Baixe o código fonte do torque aqui. Descompacte em um lugar apropriado (por exemplo, /usr/src) e execute o seguinte comando:

./configure --with-rcp=scp --disable-gcc-warnings

No Debian o script de configuração não consegue encontrar a Tk, por isso se você quiser compilar o Torque as ferramentas gráficas use o seguinte comando:

./configure --with-rcp=scp --with-tk=/usr/lib/tk8.4/ --disable-gcc-warnings

Depois de configurado é só digitar:

make && make install

Depois de instalar o sistema é preciso colocar os scripts de inicialização no /etc/init.d para que o torque seja iniciado no boot. Dentro do diretório fonte do torque digite:

cd contrib/init.d

Em seguida é preciso copiar os scripts xxx.pbs_mom, xxx.pbs_sched e xxx.pbs_server para /etc/init.d, onde xxx pode ser debian, suze ou nada para a versão redhat do script. Abaixo está um comando que copia e renomeia os arquivos na versão para o debian/ubuntu:

for i in debian.pbs_*; do cp $i /etc/init.d/`echo $i | sed 's/debian.//'`; done

Por fim:

update-rc.d pbs_server defaults

update-rc.d pbs_mom defaults

update-rc.d pbs_sched defaults

Isso finaliza a instalação  do Torque, o passo seguinte é configura-lo. Você precisa editar dois arquivos, primeiro o arquivo /var/spool/torque/server_priv/nodes, nesse arquivo é incluída a lista de nós no cluster, como eu tenho somente um nó, eu coloquei a seguinte linha:

/var/spool/torque/server_priv/nodes:

tangerine  np=2

Na primeira coluna têm-se o nome do nó (é importante que o arquivo hosts em todo o cluster contenham esse nome e que esse nome seja o mesmo do hostname.
Esse arquivo só precisa ser configurado no headnode, como tenho somente um computador, o meu headnode é o mesmo que o compute node.

Em seguida configura-se o pbs_mom em cada um dos compute nodes, de novo, no meu caso é no mesmo computador. Edite o seguinte arquivo:

/var/spool/torque/mom_priv/config:

$pbsserver      tangerine
$logevent       255

O pbsserver é o nome do headnode e o valor de logevent eu peguei da documentação do torque. Feita a configuração os deamons podem ser iniciados:

/etc/init.d/pbs_server start

/etc/init.d/pbs_mom start

/etc/init.d/pbs_sched start

(Você śo precisa iniciar os deamons essa primeira vez pois já configuramos para inicia-los no momento do boot)

Verifique se os daemons estão rodando utilizando o comand:

ps awx | grep pbs

Se você não ver os três daemons rodando é porque deu problema, olhe nos logs dos daemons em /var/spool/torque/server_logs e /var/spool/torque/mom_logs.

Com os daemons rodando podemos fazer o último passo que é criar as filas de execução, abaixo está a receita para criar uma fila chamada secondary:

qmgr -c "set server scheduling=true"
qmgr -c "create queue secondary queue_type=execution"
qmgr -c "set queue secondary started=true"
qmgr -c "set queue secondary enabled=true"
qmgr -c "set queue secondary resources_default.nodes=1"
qmgr -c "set queue secondary resources_default.walltime=1000:00:00"
qmgr -c "set server default_queue=secondary"
qmgr -c "set server scheduling=true"

As pipelines do Corona utilizam como duas files, a primeira chamada secondary que acabamos de criar e uma segunda chamada tracking que é criada com a seguinte receita:

qmgr -c "create queue tracking queue_type=execution"
qmgr -c "set queue tracking started=true"
qmgr -c "set queue tracking enabled=true"
qmgr -c "set queue tracking resources_default.nodes=1"
qmgr -c "set queue tracking resources_default.walltime=1000:00:00"

Pronto!!!!

Podemos submeter jobs ao Torque, como por exemplo esse aqui:

echo echo echo | qsub

A execução do job pode ser acompanhada pelo comando qstat, porém nesse caso geralmente job é executado instantaneamente. Quando esse job terminar serão criados dois arquivos no diretório atual, STDIN.o??? e STDIN.e??? onde ??? é o id do job. Se esses arquivos não forem criados, existe um problema no toruqe.

Um dos problemas mais comuns na execução de jobs no torque é a falta de acesso sem senha via ssh, é preciso que o ssh faça um login automático para que os resultados sejam copiados de volta, mas isso é assunto para um próximo post.

Utilizando O MosaikAligner

O Mosaik é um programa de mapeameamento com uma interface muito boa. Uma característica do Mosaik é trabalhar com arquivos binários em formato próprio, por isso é preciso pré-processar tanto a seqüencia de referência quanto os reads.

Existe um bug no processamento de seqüências em color space, por isso é necessário aplicar o patch descrito nesse bug report antes de utilizar o programa.

Para processar a referência utilize o seguinte comando:

MosaikBuild -fr fruitfly.fa -oa fruitfly.ref.bs.dat

Que vai gerar a referência em Base-Space. É preciso também gerar a referência em Color-Space, para isso deve-se utilizar o comando:

MosaikBuild -fr fruitfly.fa -cs -oa fruitfly.ref.cs.dat

Por fim, os reads são processados com o comando:

MosaikBuild -q SRR034220.fastq.gz -out SRR034220.dat -st solid

Existe mais um pre-processamento opcional pois o Mosaik pode fazer a busca das sementes utilizando uma hash-table ou então o que eles chamam de Jump Database. Para gerar o Jump Database utilize o seguinte comando:

MosaikJump -ia fruitfly.ref.cs.dat -out fruitfly.ref.cs_15.dat  -hs 15 -mem 6

A opção -mem é muito importante, pois caso ele não seja ajustada o programa vai utilizar toda a memória do sistema, no caso, eu coloquei 6 GB para o programa. Agora podemos ir apara o alinhamento:

MosaikAligner -in SRR034220.dat -out SRR034220_align.dat -ibs fruitfly.ref.bs.dat -ia fruitfly.ref.cs.dat -hs 15 -mm 4 -j fruitfly.ref.cs_15.dat -p 2

Note que estamos utilizando tanto a referência em color-space e em base-space, além disso o tamanho do hash (-hs 15) tem que bater com o valor passado para o MosaikJump, além disso estamos aceitando até 4 mismatchs (-mm 4). Estou também utilizando os cores que eu tenho. O problema é que com esse comando eu consigo alinhar somente 90 reads/s, e portanto demoraria mais de 200 horas para alinhar todo o dataset, por isso acabei utilizando a opção -bw 19 para utilizar o algorithmo de Smith-Waterman com bandas, o manual sugere utilizar 13 para um reads 36, eu cheguei ao valor de 19 através de uma regra de 3. Eu também utilizei o parâmetro -mhp 100 para limitar o número de posições por seed. Portanto o comando final foi:

MosaikAligner -in SRR034220.dat -out SRR034220_align.dat -ibs fruitfly.ref.bs.dat -ia fruitfly.ref.cs.dat -hs 15 -mm 4 -j fruitfly.ref.cs_15.dat -p 2  -bw 19 -mhp 100

Com esse comando eu consigo alinhar em torno de 260 reads/s, uma boa melhora. Porém  depois de um tempo o  valro cai para 40 reads/s. Acredito que o fato de ter somente 8GB de RAM está limitando o desempenho do programa.

Mapeando reads de Transcriptoma com Bowtie

Já existem diversos datasets de SOLiD disponíveis no banco de dados do SRA (Sequence Read Archive) do NCBI. Para testar o bowtie eu peguei um dataset de transcriptoma de D. Melanogaster SRX015641. Esse arquivo contém o resultado de somente uma corrida, mesmo assim são 4.2 Gb (O site do NCBI permite utilizar um plugin chamado Aspera Connect para facilitar o download).

O genoma da Drosófila pode ser baixado no BDGP. O genoma vem em diversos arquivos, antes de fazer o alinhamento é preciso juntar todas as partes em um único arquivo fasta utilizando o velho e bom cat:


cat na*.RELEASE5 > fruitfly.fasta

Em seguida é feita a indexação do do genoma utilizando o comando bowtie-build, supondo que o bowtie está no path, digite:

bowtie-build -C fruitfly.fasta  fruitfly

A opção -C é muito importante, ela diz para gera um índice em color-space. Esse processo demora alguns minutos. Por fim, é só executar o bowtie, como o download já está no formato fastaq não é preciso fazer nenhuma conversão. Como o arquivo está compactado, podemos fazer a descompactação em paralelo ao alinhamento com o seguinte comando:


zcat SRR034220.fastq.gz |  
    bowtie -S -C -p 2 fruitfly - | 
    gzip > SRR034220.sam.gz

O comando deve ser digitado em uma única linha. A opção -S pede para o programa gerar o alinhamento em formato SAM, a opção -C diz que estamos utilizando color-spae e -p 2 pede para o programa utilizar duas threads (que é o número de cores que eu tenho). Nessa forma são utilizadas duas pipes para descomprimir os reads e comprimir o arquivo SAM, o uso do gzip reduz muito o uso de disco e o IO, note o - após o nome da referência, ele diz para o bowtie ler os reads do STDIN e portanto o resultado do do zcat.

Alternativamento pode-se decompactar o arquivo fastq (você vai precisar de 16 de disco para fazer isso) e depois fazer o mapeamento com esse comando:


bowtie -S -C -p 2 fruitfly SRR034220.fastq |
         gzip > SRR034220.sam.gz

O bowtie possui a opção --mm que faz o programa utilizar o função memmap para ler os reads, porém, para o tamanho de input que o SOLiD gera, essa opção supera o limite de 4GB e faz com que o bowtie tenho um Segmentation Fail.

O mapeamento demorou 1h45 (wall time) e mapeou 30% dos reads.

Mapeando reads do SOLiD com mais programas

Consegui mapear os reads do dataset de E. coli com o Mosaik, com o bowtie e o com o perm (detalhes no último post).  O bowtie funcionou sem problemas e o Mosaik funcionou após aplicar o patch indicado nesse bug report. Um problema em todos os programas foi a identificação de mate-pairs, todos eles tiveram problemas em identificar os mates, o mapeamento funcionou mas na hora de fazer o match dos mates os programas tiveram problemas, preciso descobrir o que há de errado. 

Tentei usar o bfast, mas achei ele meio sacal de usar, é preciso especificar manualmente quais os padrões de discontinuous words que se deseja usar, um por um, talvez depois eu de mais uma olhada nele.

Utilizando o PerM

PerM é um alinhador de sequencias curtas bastante rápido e que alinha sequencias tanto color-space quando base-space. Dentre as vantagens desse programa estão o output em formato SAM e o uso de multiplos processadores.

Para mapear um conjunto de dados de mate-pair, por exemplo o conjunto de exemplo de E. coli disponível no site do solid tools: ecoli2x50,  basta utilizar o seguinte comando:

PerM DH10B_WithDup_FinalEdit_validated.fasta 
          -1 Rosalind_20080729_2_Chris5_F3.csfasta 
          -2 Rosalind_20080729_2_Chris5_R3.csfasta 
          -o mates.sam 

Estou suponto que o PerM está no path e que o diretório atual é o que contém os reads (as quebras de linha não devem ser utilizadas, foram colocadas somente para aumentar a legibilidade).

Para fazer o alinhamento é somente isso. Para visualizar o resultado o IGV é um ótima opção, mas para abrir o arquivo nesse programa é preciso primeiro convertar o arquivo SAM em BAM indexado. Para isso utiliza-a o samtools e as seguintes etapas:

1. Convertar o arquivo de SAM para BAM:
   
   samtools view -bt DH10B_WithDup_FinalEdit_validated.fasta mates.sam > mates.bam
   
2. Ordenar o arquivo BAM:
   
   samtools sort [-m MEM] mates.bam mates.sorted
   
   É importante ajustar o valor de MEM para que a você não tenha que esperar algumas horas pela ordenação dos reads. Esse valor deve ser informado em bytes.
   
3. Indexar o arquivo BAM:
   
   samtools  index  mates.sorted.bam 
   
   Essa última etapa vai gerar o arquivo mates.sorted.bam.bai.
   

Depois de executada a conversão e a indexação, abra o IGV. Na primeira vez que você for visualizar o genoma da E. Coli será preciso importa-lo, vá no menu File/Import Genome...

Depois de importar o genoma vá em File/Load From File... e carregue o arquivo mates.sorted.bam.

Pronto! Se você aumentar o nível de zoom vai ser possível ver todos os alinhamentos, como na figura abaixo:

O alinhamento demora em torno de 15minutos no meu Atlhom X2 3000+. A conversão de SAM para BAM e a indexação demoram mais tempo que o alinhamento.

O único problema é na resolução dos mate-pairs, em pouquíssimas sequencias o PerM consegue fazer o pareamento dos mates.