Bioinfo BR

Vaga de bioinformata na Life Technologies

2012-02-16T11:32:00.002-02:00

Oi Pessoal,

Estamos abrindo uma vaga para bioinformata aqui na Life Technologies
(ex Applied Biosystems). Essa vaga é para trabalhar com suporte de
bioinformática para os usuários dos sequenciadores de nova geração
SOLiD, 5500 e Ion Torrent. Prefere-se candiados com pós em
bioinformática e experiência em next generation sequencing analysis,
experiência com Linux e proficiência em inglês é mandatório.

Os interessados podem se inscrever no site:

http://www.lifetechnologies.com/br/en/home/about-us/careers.html

Basta ir no box de Job Search e colocar o código 9501BR

PS: Temos uma segunda vaga no México, quem for espano-hablante e
estiver interessado em se mudar pode se inscrever para a vaga 9372BR.

Ion Torrent: The Open Source Sequencer

2012-02-03T17:19:00.000-02:00

Apartir da versão 2.0 o PGM diponibiliza o código fonte da pipeline da análises, do mapeador, do Variant Caller e do utilitário TorrentScout no github:

Torrent Suite 2.0: http://github.com/iontorrent/TS

TMAP: http://github.com/iontorrent/TMAP

Torrent Scout: http://github.com/iontorrent/TorrentScout

Torrent Variant Caller Plugin: http://github.com/iontorrent/Ion-Variant-Hunter

Agora sim temos um Sequenciador Completamente open source :-)

5500 Wild Fire

2011-12-12T16:40:00.002-02:00

Vídeo sobre o 5500W (aka wild fire) que será lançado no próximo ano:

Simple SNP Calling and annotation in command line

2011-11-25T14:52:00.001-02:00

This is a simple recipe for doing SNP calling and annotation:

First, use samtools (http://samtools.sourceforge.net/) to call the SNP's:

samtools mpileup -uf ref.fa aln1.bam aln2.bam | bcftools view -vcg - > var.vcf

Them use snpEff (http://snpeff.sourceforge.net/) to annotate the vcf file:

java -jar snpEff.jar download hg19

java -jar snpEff.jar hg19 var.vcf > annotation.tsv

(If you are working with other genome change hg19 to the name of the genome, the list of available genomes are in the snpEff site).

The snpEff will inform the name of the SNP, the changed residue, codon, etc, very useful.

I tested this with an Ion Amplicon, I think should work with SOLiD too.

Hybrid Assembly Pipeline with SOLiD

2011-10-24T21:26:00.002-02:00

Foi disponiblizado hoje uma ferramenta para combinar os reads de uma biblioteca de mate-pair do SOLiD com contigs gerados por outra tecnologia (ou mesmo pelo próprio SOLID).

A idéia é mapear os pares nos contigs e depois montar um grafo de De Brujin com as possíveis ordenações dos contigs, corrigir o grafo e gerar um scaffold. Aparentemente esse programa tem um resultado muito mais acurado do que o SOPRA.

Para quem esta fazendo Denovo com o SOLiD é uma ferramenta muito útil.

Paper do Ion Torrent

2011-07-20T16:22:00.002-03:00

Saiu na Nature o paper descrevendo o Ion Torrent:

http://www.nature.com/nature/journal/v475/n7356/full/nature10242.html

Ion Torrent sequencia cepa de e. coli que está causando surto na Europa

2011-06-04T23:07:00.000-03:00

Utilizando o Ion Torrent pesquisadores identificaram que a e. coli que está causando o surto na europa é a compinação de duas outras cepas mortais:

http://www.courant.com/health/connecticut/hc-sequencing-e-coli-0603-20110602,0,2978920.story

Matéria no jornal da band:

http://videos.band.com.br/Exibir/Alemanha-diz-que-bacteria-nao-esta-em-pepino-espanhol/2c9f94b630297d9901304d534290161b?channel=587

De Novo assembly com o Ion Torrent

2011-05-26T11:03:00.000-03:00

Ótimo blog post sobre montagem de genomas utilizando dados do Ion Torrent:

http://pathogenomics.bham.ac.uk/blog/2011/05/ion-torrent-data-blog-post-a-week-is-a-long-time-in-genomics/

Notícias sobre o Lançamento do Ion Torrent

2011-04-20T11:51:00.000-03:00

No blog do Estadão:

http://blogs.estadao.com.br/link/projeto-genoma-para-todos/

No Valor Econômico:

Na Info:

http://info.abril.com.br/noticias/ciencia/novo-chip-sequenciador-de-dna-chega-ao-brasil-18042011-28.shl

Fotos do SOLiD 5500 e do Ion Torrent

2011-04-02T20:54:00.000-03:00

Algumas fotos que eu tirei no PAG 2011.

Fernando com o SOLiD 5500:

E o Ion Torrent:

Multithreading no Velvet 1.1

2011-04-01T23:48:00.000-03:00

A versão 1.1.01 do Velvet agora suporta multithreading via OpenMP, além disso segundo o anúncio houve uma melhora no uso de memória[1]

Logo em seguida foi lançado o 1.1.02, para corrigir alguns bugs da versão anterior [2]. Eu testei essa versão com a bibliteoca de DH10B feita com SOLiD, não medi o tempo, mas percebi que ocorre um real uso dos 8 cores da máquina de teste e que houve uma significativa redução no tempo velvetg.

[1] http://listserver.ebi.ac.uk/pipermail/velvet-users/2011-March/001303.html
[2] http://listserver.ebi.ac.uk/pipermail/velvet-users/2011-March/001311.html

Estatística de arquivo BAM

2011-03-10T11:29:00.000-03:00

Uma maneira fácil de extrair algumas estatísticas de um arquivo bam é utilizar o comando

samtools flagstat

Um outro programa que parece muito interessante é o samstat:

http://samstat.sourceforge.net/

Que gera um relatório HTML sobre o alinhamento.

Ion Torrent Sequencia o Genome de Gordon Moore da Intel

2011-02-25T08:18:00.002-03:00

A Ion Torrent anunciou o sequenciamento do genoma de Gordon Moore, fundador da Intel e criador da lei de Moore:

http://www.genomeweb.com/sequencing/ion-torrent-sequences-intel-co-founders-genome-finds-more-uniform-coverage

Visualizando os grafos do Velvet

2011-02-24T17:37:00.001-03:00

O velvet inclui um script para converter o arquivo LastGraph para o formato dot. O problema é que o visualizador do graphviz é bem tosco, porém é possível abrir um arquivo no formato dot no visualizador gratuíto gephi:

http://gephi.org/

O resultado da reenderização do grafo pode ser visto abaixo:

Chip de 1Gbp para Ion Torrent anunciado para a segunda metade de 2011

2011-02-23T21:04:00.000-03:00

Hoje foi anunciada uma nova versão do chip de sequenciamento do Ion Torrent, a 318, capaz de sequenciar até 1 Gbp em reads de pelo menos 300bp

Comentários adicionais sobre esse anuncio no blog Omics! Omics!

http://omicsomics.blogspot.com/2011/02/has-ion-torrent-taken-318-sized-lead.html

Ion Community

2011-02-17T14:31:00.001-02:00

Existe um site para os usuários de Ion Torrent chamado Ion Community:

http://ioncommunity.iontorrent.com/

Big Data, Big Problems

2011-02-08T14:01:00.000-02:00

Artigo muito interessante sobre estrutura computacional para analisar grandes volumes de dados:

http://cacm.acm.org/blogs/blog-cacm/103932-big-data-big-problems/fulltext

Genoma por US$ 10mil

2011-01-16T11:16:00.000-02:00

A empresa DNAVision em Bruxelas oferece comercialmente o sequenciamento de um genoma humano por US$ 10mil. A empresa recentemente adquiriu 2 SOLID's 5500xl e 2 SOLiD's 4 para essa tarefa.

Fonte: http://www.genomeweb.com/sequencing/dnavision-offer-10k-human-genome-sequencing-services-purchases-four-solids

SOLiD Paper: Identification of methylated regions with peak search based on Poisson model from massively parallel methylated DNA immunoprecipitation-sequencing data.

2010-12-16T18:09:00.000-02:00

Artigo com análise do metiloma humano utilizando uma estratégia de enriquecimento de regiões metiladas:

"Identification of methylated regions with peak search based on Poisson model from massively parallel methylated DNA immunoprecipitation-sequencing data"

Os autores verificarão que é necessário somente um quad do SOLiD para sequenciar todo o metiloma humano.

SOLiD Paper: De novo mutations of SETBP1 cause Schinzel-Giedion syndrome

2010-12-16T09:36:00.000-02:00

No artigo "De novo mutations of SETBP1 cause Schinzel-Giedion syndrome" os autores fazem o sequenciamento do exoma de 4 pacientes e descobrem 3 mutações no gene SETBP1 relacionadas com a doença.

SOLiD Paper: A small-cell lung cancer genome with complex signatures of tobacco exposure

2010-12-06T16:11:00.000-02:00

Eu vou postar aqui no blog alguns papers chave que utilizam o SOLiD em diferentes aplicações. Iniciando com o artigo

A small-cell lung cancer genome with complex signatures of tobacco exposure:
http://www.nature.com/nature/journal/v463/n7278/abs/nature08629.html

Que é um trabalho de resequenciamento de células de pulmão imortalizadas normais e tumorais. Comparando as mutações nas duas amostras os autores detecção as mutações somáticas adquiridas pelo indivíduo.

Servidor com 1 TB de RAM por US$ 100 mil

2010-11-30T16:39:00.000-02:00

O PowerEdge R910 da Dell é um servidor 4U que aceita até 1 TB de RAM, um valor que me parece muito adequado para fazer montagem de genomas. O custo básico do servidor com 128 GB de RAM é de US$ 27mil, para aumentar a RAM para 1 TB adiciona-se US$ 62mil. Aumentando também a quantidade de disco e adicionando processadores de 8C chega-se a uns US$ 100mil, o que me parece um valor viável para um servidor que em teoria pode montar o genoma de um mamífero.

SAM format specification 1.3

2010-11-24T16:44:00.000-02:00

Saiu a versão 1.3 da especificação do formato SAM/BAM:

http://samtools.sourceforge.net/SAM-1.3.pdf

A meu ver eles melhoraram a formatação do arquivo (aparentemente trocaram o Word por LaTeX) e mudaram a nomenclatura para lidar com leituras com mais de 2 tags (bibliotecas de fragmento tem 1 tag, de mate-pair ou pair-end tem 2 e as bibliotecas feitas com strobo lightining no PacBio tem diversas tags).

A de novo paradigm for mental retardation

2010-11-19T23:06:00.000-02:00

No artigo A de novo paradigm for mental retardation os autores utilzam dados de SOLiD para sequenciar o exoma de 10 trios (mãe, pai e filho/a) em que o filho/a sofrem de retardo mental. O artigo busca mutações não sinonimo não herdadas dos pais que poderiam estar ligadas ao retardo mental.

A referência é essa:

Referencia: Vissers LELM, Ligt J de, Gilissen C, et al. A de novo paradigm for mental retardation. Nature Genetics. 2010;(November).

Autor do Montador de Genomas Cortex vai estar em Curitiba em Dezembro

2010-11-19T18:36:00.000-02:00

Já faz um tempo eu ouço falar do montador Cortex, que seria um montador novo, sucessor do Velvet. As únicas que eu consegui foram comentários de pessoas que viram a palestra do autor do programa Mario Cacamo.

Felizmente dia 10 de Dezembro vai ter uma palestestra do Mario Cacamo em Curitiba:

http://www.bioinfo.ufpr.br/iwb/programme.php

Eu pretendo ir para descobrir mais detalhes diretamente.

Processing sequences with Hadoop

2010-10-27T11:04:00.000-02:00

The increasing volume of data generated by next-gen creates a huge pressure in the computer infra-structure. We need bigger and more expensive storage solutions to hold the results of the experiments and the intermediate data of the analysis, and because all the data are in a single system, we need more network bandwidth to deal with all-to-one pattern of network access (we can have a 10 GB connection to the storage, but what we do if this is not enought?).

Large internet sites like google, Facebook, etc have a similar problem with the data from their users (or from all the internet in the case of Google). To solve the storage and the processing problem they create new solutions to scale the storage and the processing to really big clusters (hundreds or thounsands of nodes). Google created the MapReduce Framework and the Google Filesystem (http://labs.google.com/papers/mapreduce.html and http://labs.google.com/papers/gfs.html), which are solutions for batch processing huge amount of data. Haddop is a open source Java Framework which implements both the Distributed filesystem and the MapReduce framework and have been used by big internet players like Yahoo and Microsoft.

In the Bifx space we have some projects using hadoop to process next-gen data:

Crossbow: http://bowtie-bio.sourceforge.net/crossbow/index.shtml

Myrna: http://bowtie-bio.sourceforge.net/myrna/index.shtml

Cloudburst: http://sourceforge.net/apps/mediawiki/cloudburst-bio/index.php?title=CloudBurst

Anyone can lear more about hadoop with the training videos from Cloudera (A hadoop consulting company):

http://www.cloudera.com/resources/?type=Training

I have been playing with Hadoop in my free type and I produced some code to deal with Biological data, more specificaly I created a InputFormat for Fasta files and a FastatoFastq converter using hadoop, the code is in github:

http://github.com/lvaruzza/bioseq-hadoop

This is more a self-study project, but anyone interested in Hadoop and distributed programming is invited to use this code.

Baylor College vai comprar 2dois servidores com 1TB de RAM cada

2010-10-22T11:27:00.002-02:00

Segundo matéria no GenomeWeb o Baylor College vai gastar US$ 262 mil na compra de dois servidores para montagem de genomas, cada um com 1 TB de RAM:

http://www.genomeweb.com/informatics/baylor-genome-sequencing-center-wins-262k-stimulus-funds-hpc-cluster

PS: Com 20 sequenciadores, o Baylor College é atualmente o maior site com SOLiD's instalados.

A mudança na geração de dados causada pelo NextGen

2010-09-10T19:21:00.000-03:00

Achei um comentário muito interessante em um artigo do Lincoln Stein:

"... the 1000 Genomes Project [25], which is cataloguing human genetic variation, deposited twice as much raw sequencing data into GenBank’s SRA division during the project’s first 6 months of operation as had been deposited into all of GenBank for the entire 30 years preceding (Paul Flicek, personal communication)."

Fonte: http://genomebiology.com/2010/11/5/207

SOLiD Probes

2010-07-19T16:33:00.000-03:00

O probe que faz a leitura da sequencia do SOLiD tem a seguinte estrutura:

XX NNN ZZZ

Onde XX é a dibase cujo o flurofito se refere, NNN são chamdas de bases degeneradas e ZZZ são bases universais (inosinas).

Eu pensava que as bases degeneradas eram algum tipo de base modificada mas na verdade são bases normais que são sintetizadas em todas as 1024 (4^5) combinações possíveis e portanto os probes anelam em qualquer sequencia, não por uma modificação na base mas sim pela forma como os probes são misturados,ou seja, a mistura que é degenerada. Com essa construção, apesar de lermos somente 2 bases por ciclo, a especificidade do probe é referente a cinco bases.

Mate Pair Paper

2010-07-19T14:22:00.000-03:00

Recentemente me perguntaram a referência da técnica de mate pair. Achei esse artigo que acredito ser a referência original da técnica:

Pairwise end sequencing: a unified approach to genomic mapping and sequencing.
Roach et al, 1995:
http://www.ncbi.nlm.nih.gov/pubmed/7601461

Introdução ao Python para Cientistas

2010-07-06T11:44:00.000-03:00

Coleção de tutoriais em python focados em aplicações para cientistas:

http://software-carpentry.org/

Nova versão da CLC Machine

2010-07-06T11:36:00.002-03:00

A CLC lançou uma nova configuração da CLC Machine:

http://www.clcmachine.com/system.php

Essa configuração me parece mais apropriada para análise de dados de SOLiD do que a configuração anterior: 12 cores, 48 GB de RAM e 8 TB de disco.

Valores de qualidade no SOLiD

2010-07-05T16:14:00.000-03:00

Saiu um application notes bem interessante sobre detecção de beads policlonais em uma corrida de SOLiD:

http://bioinformatics.oxfordjournals.org/cgi/content/abstract/26/6/849

Nos reads normais a qualidade cai na extremidade 3', porém nos reads policlonais a qualidade é ruim em todo o read, portanto se for feita uma filtragem pelos reads com baixa qualidade no incio, por exemplo nas 10 primeiras bases, é possível filtrar os reads policlonais.

Outro resultado interessante do artigo é que a maioria dos erros de mapeamento se concentra em bases com QV < 10.

OZZY no SOLiD

2010-06-29T22:41:00.003-03:00

E não é que vão fazer o Genoma do Ozzy Osbourne no SOLiD:

http://www.bizjournals.com/losangeles/stories/2010/06/28/daily5.html

PS: Na notícia eles só falam que vai ser feito em um sequenciador da Life Tech, mas acho seguro supor que eles estão falando do SOLiD.

Update: O artigo abaixo confirma que o sequenciamento vai ser feito no SOLiD 4:

http://www.bio-itworld.com/els/2010/06/30/ozzys-genome.html

Tutorial sobre o ABySS

2010-06-29T08:36:00.000-03:00

Achei um tutorial sobre a execução do montador distribuído ABySS:

http://ged.msu.edu/angus/tutorials/short-read-assembly.html

É muito mais fácil ter 512 GB ou 1 TB de memória distribuído em um cluster do que em uma única máquina, por isso acredito que montadores Eurelianos distribuídos são o caminho para a montagem De Novo de Eucaritos com short reads.

PS: O ABySS trabalha com dados do SOLiD em Color Space.

99.94% de Acurácia

2010-06-18T17:00:00.000-03:00

Uma acurácia de 99.94% significa que esperamos 600 diferenças entre os reads e a referência a cada 1 milhão de bases sequenciadas. Alguém pode argumentar que a diferença entre 99.90% e 99.94% é insignificante, mas em um sequenciamento de 100 GB isso significa que o SOLiD vai ter 40 milhões a menos de mismatchs, isso são quase 9 genomas de E. coli.

2010-06-18T09:16:00.000-03:00

Na lista de discussão do ABySS eu vi um post com o seguinte trecho:

...
The haploid genome size of the species whose genome I'm trying to
assemble is 3 billion bp. My colleague doesn't want me to name the
species so I need to honor his request. Let's just call it "rb".

The Illumina pe-libs are:

lib1: 10x @ 2x75bp
lib2: "20x" @ 2x150bp

lib1: 210 million read-pairs
lib2: 200 million read-pairs

The fragment sizes of both lib1 and lib2 are "500-600" bp.

Our experience with 2x150bp reads (an Illumina-supported trial) was
disappointing. Your mileage may vary. Thus, the lib2 spec "2x150bp"
is nominal. We don't expect much more than half that (2x75bp) to be
usable.
...

Esse comentário é bastante pertinente na questão do tamanho do read versus acurácia no sentido de que aumentar o tamanho sem aumentar a acurácia não faz sentido porque o ruído inserido torna a montagem inviável (no contexto de montagem Eureliana, nos algoritímos anteriores os problemas são outros).

Brasil no TOP500

2010-06-02T16:37:00.000-03:00

O Brasil está bem colocado na lista do TOP500 do mês de junho de 2010, 86a posição:

http://www.top500.org/lists/2010/06

Com um cluster da SUN instalado no NACAD da UFRJ:

http://www.nacad.ufrj.br/

Early History of GenBank

2010-04-07T11:50:00.000-03:00

Achei alguns textos interessantíssimos sobre o início do GenBank, quando ele ainda estava no Los Alamos National Laboratory e o NCBI ainda nem existia. Os textos são uma entrevista no site do LANL e um artigo escrito por Walter B. Goad.

Em 1979 eles já viam a importância da bioinformática e da análise de dados para a biologia molecular, eles até sugerem uma abordagem de "serviços remotos" para lidar com heterogeneidade dos computadores que existiam na época.

Gerando consenso a partir de um arquivo SAM

2010-03-31T16:16:00.000-03:00

Uma coisa muito legal do formato SAM é o utilitário samtools, ele permite fazer diversos processamentos muito úteis, sendo um deles gerar um consenso com base no alinhamento.

Primeiro é necessário converter o arquivo SAM em SORTED BAM, (uu já abordei esse assunto no post sobre Perm). Em seguida utiliza-se os comandos samtools pileup e o samtools.pl pileup2fq em conjunto:

samtools pileup -cf | samtools.pl pileup2fq -D100 > consensus.fastq

A opção -D100 limite a cobertura a no máxio 100 reads, ajuste esse parâmetro de acordo com a cobertura do seus dados.

Configurando o Maui, alternativa ao pbs_sched

2010-02-09T17:37:00.000-02:00

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

2010-02-05T08:31:00.006-02:00

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)

2010-02-03T10:48:00.002-02:00

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 resources_available.nodect = 20"

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 torque.

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

2010-02-02T19:31:00.006-02:00

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

2010-02-01T12:01:00.009-02:00

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 GB 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

2010-01-29T00:26:00.000-02:00

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

2010-01-28T13:57:00.002-02:00

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:

Convertar o arquivo de SAM para BAM:

samtools view -bt DH10B_WithDup_FinalEdit_validated.fasta mates.sam > mates.bam
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.
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.

De Novo Acessory Tools

2009-12-17T15:20:00.000-02:00

Foi lançada uma nova versão dos utilitários para fazer montagem De Novo com SOLiD reads. Essa nova versão está disponível no site:

http://solidsoftwaretools.com/gf/project/denovotools/

PS: A versão antiga ainda está disponível no site (http://solidsoftwaretools.com/gf/project/denovo/), portanto verifique qual versão você está baixando.

Visualização de Montagens

2009-12-17T12:04:00.000-02:00

Achei um visualizador de assemblies muito bom chamado Tablet. O programa é feito em java, testando em Mac, Win e Linux. A user interface é bastante agradável e achei a performance, para o exemplo pequeno disponível na página do programa, muito boa.

A página do programa é a seguinte:

http://bioinf.scri.ac.uk/tablet/index.shtml

O programa esta descrito no artigo Tablet – Next Generation Sequence Assembly Visualization.

Uso de memória pelo Velvet

2009-12-16T15:40:00.004-02:00

O fator limitante para a execução do Velvet é a quantidade de memória disponível no computador. Segundo essa
thread no fórum SEQAnswers forum a fórmula para calcular o uso da memória em kb é:

Ram required for velvetg = -109635 + 18977*ReadSize 
    + 86326*GenomeSize + 233353*NumReads - 51092*K

Onde:

ReadSize está em bases
GenomeSize está em megabases
NumReads está em milhões de reads
E k é o tamanho do k-tamero utilizado pelo Velvet

Para facilitar a vida dos usuários do Velvet (e para aprender a usar o Google App Engine), eu criei uma aplicaçãozinha que faz essa conta. O endereço é o seguinte:

http://denovoutils.appspot.com/

Next-gen Map

2009-12-07T11:51:00.000-02:00

om a lista Achei esse site com um mapa bem legal com a lista dos sequenciadores de próxima geração:

http://pathogenomics.bham.ac.uk/hts/

A lista está bastante desatualizada para o Brasil, mas os próprios usuários podem atualiza-la clicando diretamente no mapa.

Hadoop for bioinformatics

2009-12-02T17:22:00.003-02:00

Hadoop for Bioinfomatics - Deepak Singh from Cloudera on Vimeo.

Tutorial: Mapeando Reads utilizando Hadoop e EC2

2009-11-09T11:10:00.000-02:00

Hoje foi postado um tutorial no blog Biodivertido sobre a instalação e uso do programa de mapeamento CloudBurst em uma instância do Amazon EC2 (Elastic Cloud Computing).

Hadoop, Cloud Computing e Amazon EC2 são três coisas que eu acredito terão muita importância para os usuários do sequênciadores de nova geração.

De Novo Assembly

2009-10-29T12:02:00.002-02:00

O trabalho An Eulerian path approach to DNA fragment assembly propôs um nova abordagem para a montagem de sequencias que permite a montagem De Novo de genoma pequenos utilizando os reads de SOLiD ou Solexa.

Segue abaixo a lista de programas que implementam esse metodo:

Velvet	Site	Paper
Edena	Site	Paper
EULER-SR	Site	Paper
SSAKE	Site	Paper
ADiR	Site
ALLPATHS2	Site	Paper
ALLPATHS	Site	Paper
QSRA	Site	Paper
VCAKE	Site	Paper
MIRA	Site
Abyss	Site	Paper1, Paper2

Hadoop na bioinformática

2009-10-27T12:37:00.000-02:00

O Hadoop é a combinação de um framework open-source em Java para processamento distribuído e um filesystem também distribuído que foi desenvolvido nos moldes do MapReduce e GFS do Google. Essas tecnologias permitem que se desenvolvam programas que podem ser escalados para a execução em clusters com dezenas de milhares de computadores. Atualmente o Hadoop é a arquitetura por trás do sistema de busca do Yahoo.

O aumento brutal no volume de dados gerado pelos sequenciadores de nova geração cria novos desafios para a bioinformática. O hadoop é uma solução bastante adequada para muitos dos problemas ligados com a análise desses novos dados de sequenciamento. No blog da Cloudera foi colocado um artigo sobre o Crossbow, que é um pipeline para resequenciamento.

Além de permitir uma grande escalabilidade, o Hadoop é suportado pelo sistema de Cloud Computing Amazon EC2. Com o cloud computing é possível executar análises pesadas sem investir um grande volume de dinheiro na compra de computadores, pois cria-se um cluster on-demand no data center da Amazon e paga-se por hora de CPU usada.

Acredito que no futuro devem surgir mais pipelines que utilizam o Hadoop.

Lista de programas para mapear short sequencies

2009-10-27T08:35:00.011-02:00

A nova geração de sequênciadores de alto desempenho geram uma quantidade enorme de reads muito menores do que o método de Sanger. Alinhar centenas de milhões de seqüências curtas é um grande desafio computacional, por isso estão sendo criados diversos programas para essa tarefa.

Abaixo está uma lista de programas que eu levantei:

Nome	Site	Paper
SOAP	http://soap.genomics.org.cn/	SOAP2: an improved ultrafast tool for short read alignment
MapReads	http://solidsoftwaretools.com/		cs
SHRiMP	http://compbio.cs.toronto.edu/shrimp/	SHRiMP: Accurate Mapping of Short Color-space Reads	bs, cs
perM	http://code.google.com/p/perm/	PerM: efficient mapping of short sequencing reads with periodic full sensitive spaced seeds	bs,cs
Mosaik	http://code.google.com/p/mosaik-aligner/		cs,bs
bwa	http://bio-bwa.sourceforge.net/	Fast and accurate short read alignment with Burrows-Wheeler Transform	cs,bs
maq	http://maq.sourceforge.net/index.shtml
bowtie	http://bowtie-bio.sourceforge.net/index.shtml	Ultrafast and memory-efficient alignment of short DNA sequences to the human genome	bs,cs
SOCS	http://socs.biology.gatech.edu/	Efficient mapping of Applied Biosystems SOLiD sequence data to a reference genome for functional genomic applications	cs
bfast	https://secure.genome.ucla.edu/index.php/BFAST	BFAST: An Alignment Tool for Large Scale Genome Resequencing	cs,bs
RMAP	http://rulai.cshl.edu/rmap/	Using quality scores and longer reads improves accuracy of Solexa read mapping
CloudBurst	http://sourceforge.net/apps/mediawiki/cloudburst-bio/index.php?title=CloudBurst

Extração de tags de SAGE

2005-08-04T15:11:00.000-03:00

Estou analisando alguns programas que estraem as tags de SAGE das seqüências. O processo consiste em extrair as ditags que estão entre os sítios CATG. Já analisei dois programas, o SAGEparser e o tagcalling. O SAGEparser tem um sistema simples de corrigir o número de ditags apartir da freqüência de susas monotags. Já o tagcalling possui um sistema de EM que corrige a contagem das tags considerando a possibilidade de mutações. Esse segundo programa me parece melhor.

Cuidado com o CAP3

2005-07-15T10:22:00.000-03:00

Quem for utilizar o cap3, tome cuidado com a versão do programa. Digite cap3 e veja se ele imprime a seguinte linha:

VersionDate: 04/15/05

A versão antiga do cap3 não indica qualquer informação sobre a versão. Essa versão antiga contém um bug em linux, não sei se ele se repete em outras plataformas. Se vc for montar uma seqüência e resolver mascarar as pontos com N, o cap3 antigo vai gerar um número muito grande de singlets, a versão nova não faz isso.