Análise de Microarrays Usando o R
Data de criação: 2005-04-12.
Palavras-chave: microarrays, Acadêmico, R-Project.
1. Introdução
Descrevemos nessa página os detalhes da implementação computacional das ferramentas utilizadas para análise de microarrays de cDNA, desenvolvidas sob projeto de iniciação científica pela FAPESP.
As ferramentas foram desenvolvidas tendo em mente basicamente o problema prático estudado no Laboratório de Cardiologia Molecular do Instituto do Coração (InCor-USP), mas foram implementadas de forma a serem facilmente aplicadas à outros conjuntos de dados.
Todos os programas foram desenvolvidos baseando-se em um sistema operacional GNU/Linux, mas não há grandes dificuldades em portá-los para outro sistema.
2. Filtragem e armazenamento (banco de dados)
Decidimos em primeiro lugar, armazenar todos os dados iniciais do problema (antes da análise) em um banco de dados externo ao R. Tal decisão foi tomada levando em conta que o ponto forte do R é a análise e interpretação dos dados e seu desempenho não é particularmente bom na manipulação de grandes conjuntos de dados.
A estrutura do banco de dados foi adotada seguindo as linhas gerais descritas em [1]: mantemos uma tabela para cada hibridização (contendo os dados de intensidade para cada canal), e uma tabela para cada tipo de lâmina utilizada (contendo quais seqüências estão representadas em cada spot). Adicionalmente, colocamos uma tabela para guardar informações referentes às amostras de mRNA e às recíprocas das hibridizações.
Definida a estrutura do banco de dados, implementamos ele em MySQL [1], através de alguns scripts em Perl [2], desenvolvidos para esse fim. Apesar da escolha de um banco de dados específico (MySQL), devido aos pacotes DBI (Database Independent Interface) do Perl e do R, a adoção de outro banco de dados (como Postgres) é extremamente simples, sendo necessária a mudança de apenas algumas partes do código desenvolvido.
Criando o banco
Precisamos em primeiro lugar de um banco chamado 'microarray' dentro do MySQL. Isso é feito com o comando
A partir daí criamos algumas tabelas que vão ser utilizadas pelos nossos filtros para incluir as informações sobre as lâminas:
mysql> CREATE TABLE lamina (name VARCHAR(50) PRIMARY KEY,pixel_size INT(5),x_origin INT(10),y_origin INT(10),type VARCHAR(100),ratio_form VARCHAR(100));
mysql> CREATE TABLE add_info (name VARCHAR(50),reciproca VARCHAR(50),cy3_desc VARCHAR(200),cy5_desc VARCHAR(200),cy3_control INT(1),lamina VARCHAR(50));
Inserindo os dados
Dadas essas tabelas iniciais criadas, a inserção dos dados é feita de forma automatizada através de scripts em Perl. Para acrescentarmos um arquivo no formato GIPO (GAL), indicando as seqüências de cDNA contidas nos spots, usamos o script gipo_parser.txt. Vamos ilustrar sua sintaxe com um exemplo. Digamos que queiramos incluir o arquivo 11ka.gal e vamos usar como referência o nome 'noruega'. Basta rodar o script da seguinte maneira:
Para leitura dos dados de intensidade gerados pelo GenePix, temos o script gpr_parser.pl. Ilustramos seu uso novamente através de um exemplo. Digamos que queiramos inserir os dados de leitura de intensidade para a hibridização rat85 que estejam salvos no arquivo rat85.gpr:
De forma análoga, temos um script para leitura dos dados de intensidade gerados pelo Dapple, o dapple_parser.pl. Seu uso é idêntico ao do script anterior.
Algumas informações ainda podem ser completadas nos bancos, em especial na tabela add_info, para indicar informações como quais amostras de mRNA foram fixadas a quais corantes. Isso pode ser feito por exemplo, diretamente no MySQL:
Em particular é interessante acrescentar essas informações para aproveitar as funções que lidam diretamente com os spots, indicando suas recíprocas e os nome das seqüências neles contidas.
3. Interagindo com o R
Depois que as informações sobre os experimentos, incluindo os arquivos de intensidade produzidos pelo GenePix/Dapple, foram inseridas no banco de dados, temos uma série de classes e funções em R para manipulá-las.
As funções que fazem essa interação estão definidas no arquivo maexp.R. Para carregá-las, basta digitar dentro do R:
Trabalhando com lâminas inteiras
Uma das primeiras funções a serem utilizadas é a craw. Com elas criamos uma representação dos dados de uma hibridização no R. Por exemplo:
Com o comando acima criamos o objeto rat85.raw, que é uma representação das intensidades por spot para a hibridização rat85. A partir desse objeto, podemos fazer vários gráficos e análises descritivas, como por exemplo através do uso das funções mads,desvs,medias,medianas, criadas especialmente para atuar em objetos desse tipo.
Ao digitar o comando:
Será feito um gráfico MAPlot para as intensidades da hibridização rat85, como mostra a figura abaixo.
Podemos ainda gerar um gráfico de boxplots de M por subarray, com o comando:
Gerando a figura abaixo:
Trabalhando com spots individuais
É possível durante uma análise que tenhamos interesse em verificar como se comportam alguns poucos spots. Para facilitar essa tarefa criamos a classe spot e algumas funções que permitem trabalhar com objetos dessa classe.
A primeira delas é a cpot(lamina,block,row,column), que cria um objeto da classe spot que esteja na hibridização lamina, no subarray bloco, na coluna column e na linha row.
Com o comando abaixo criamos um objeto da classe spot com o nome de spot1, localizado na hibridização rat88, bloco 30, coluna 10 e linha 1.
Podemos então através da função summary(), obter informações sobre ele:
Spot summary information
array: rat88 subarray: 30 row: 10 column: 1
gene_id: Rn.9436
description: ESTs, Highly similar to (defline not available 6180013) [H.sapiens]
Corrected intensities (foreground - background):
Mean Median Description
Cy 3 29381 32901 SHR sal
Cy 5 19074 19074 2c sal
Ou ainda informações sobre sua réplica:
Spot summary information
array: rat88 subarray: 30 row: 9 column: 6
gene_id: Rn.9436
description: ESTs, Highly similar to (defline not available 6180013) [H.sapiens]
Corrected intensities (foreground - background):
Mean Median Description
Cy 3 29442 32776 SHR sal
Cy 5 19426 19426 2c sal
Podemos ainda ver a imagem do spot na lâmina, a fim de verificar se não houve um estouro na leitura ou problemas na hibridização:
Outras funções para trabalhar com spots individuais também estão disponíveis, como adv.replica ou reciproca.
Trabalhando com um conjunto de spots
Quando estamos fazendo uma análise exploratória, pode ser útil poder trabalhar com um conjunto de spots que desejamos estudar. Por exemplo no caso em que tenhamos um conjunto de spots controle e queiramos verificar como está a expressão desses controles entre as diferentes hibridizações.
Há várias maneiras de criar um conjunto de spots para trabalhar. Podemos usar o comando subset do R, para escolher spots que tenham certas medidas de M e A. De maneira análoga podemos selecionar spots que tenham certas seqüências de cDNA específicas, ou ainda, escolher spots a partir de um gráfico MAPlot.
Para facilitar essa última opção, temos a função setfinder(madat), que recebe como argumento um objeto de intensidades e transforma o cursor do mouse em um selecionador de spots. A partir daí, basta clicarmos com o botão esquerdo sobre os spots que estejamos interessados no MAPlot. Ao terminarmos a seleção, clicamos com o botão direito do mouse, sendo retornada então uma lista com os spots escolhidos.
Dada essa lista, temos um conjunto de funções análogas às criadas para spots individuais (internamente as funções para conjuntos de spots usam as funções para spots individuais, de forma a reutilizar o código), como setdesc(brc,madat), setview(brc,madat), setexpr(brc,madat).
Dado um conjunto de spots selecionados guardados no objeto rat85.brc.matched, ilustremos algumas das funções acima.
> rat85.brc.matched
block row column
412 1 15 6
436 1 16 1
6155 8 10 7
6179 8 11 2
7703 10 5 18
7727 10 6 13
7991 10 15 16
8015 10 16 11
24627 30 9 6
24651 30 10 1
Acima apenas listamos o conteúdo do conjunto de spots.
> setexpr(rat85.brc.matched,rat85.raw)
block row column m a flags
412 1 15 6 2.409810 13.65147 0
436 1 16 1 2.372897 13.55222 0
6155 8 10 7 1.980907 11.28392 0
6179 8 11 2 2.060793 11.41410 0
7703 10 5 18 2.328872 12.44463 0
7727 10 6 13 2.334672 12.56875 0
7991 10 15 16 2.426450 13.10539 0
8015 10 16 11 2.386922 13.07610 0
24627 30 9 6 2.252549 13.76584 0
24651 30 10 1 2.545238 14.22881 0
Acima temos os níveis de expressão para cada um deles na hibridização rat85.
> setexpr(rat85.brc.matched,rat89.raw)
block row column m a flags
412 1 15 6 -1.601081 13.63974 0
436 1 16 1 -1.600253 13.56189 0
6155 8 10 7 -1.125265 10.38281 0
6179 8 11 2 -1.094673 10.85568 0
7703 10 5 18 -1.354174 12.45187 0
7727 10 6 13 -1.316535 12.66915 0
7991 10 15 16 -1.354358 13.34761 0
8015 10 16 11 -1.328311 13.21466 0
24627 30 9 6 -1.336396 14.38395 0
24651 30 10 1 -1.404151 14.35570 0
Acima temos os níveis de expressão para cada um deles na hibridização rat89 (recíproca da rat85).
Há funções ainda para fazer a interseção entre duas listas de spots (intersectEXPR(brc1,brc2), útil na seleção de spots entre diferentes hibridizações), para transformar um conjunto de intensidades madat em um conjunto de spots (inten2set(madat)) ou ainda para transformar um objeto da classe spot num objeto do tipo lista de spots (spot2set(spot)).
A função maplot(madat,...) mencionada logo no começo da seção, pode ser ainda utilizada com um argumento opcional highlight=brc, de forma que os spots indicados pelo conjunto de spots brc são destacados no gráfico MAPLot. Assim, o comando:
Gera o gráfico:
onde os pontos em azul são os spots listados em rat85.brc.matched.
Normalizações
No arquivo norm.R implementamos alguns dos métodos de normalização encontrados na literatura ([3] e [4] por exemplo). Todos atuam em cima de objetos do tipo madat, gerados a partir do comando craw(lamina) ou a partir de outra normalização.
Normalizações globais
Para normalização global de um objeto de intensidades, temos a função global.norm(madat,flaglevel). Através de madat apontamos qual o objeto de intensidades que queremos usar, e através de flaglevel, apontamos (opcionalmente) se queremos deixar de fora spots que foram marcados como 'ruins' pelo software de leitura das imagens. Para normalizar as leituras de intensidade guardadas em rat88.raw, fazemos:
Em rat88.global ficam salvas portanto as intensidades normalizadas. Podemos então aplicar as funções normais para lidar com objetos de intensidade. Por exemplo, para ver o MAPlot, fazemos:
Obtendo:
Para obter uma normalização global por subarray o processo é o mesmo, mas através da função subarray.norm(madat,flaglevel).
Normalizações por lowess
Para normalização por lowess global/por subarray temos as funções global.lowess(madat,lowess.f,flaglevel) e subarray.lowess(madat,lowess.f,flaglevel) respectivamente. Elas funcionam de forma analoga às normalizações globais, com a diferença do parâmetro opcional lowess.f, que indica o parâmetro f a ser passado para o ajuste de lowess.
Ajuste de escala
Para fazer o ajuste de escala entre os subarrays, temos a função within.slide.scale(madat). Ela retorna outro madat com a normalização por escala utilizando o desvio médio absoluto, como descrito em [3].
Autonormalização
A função self.normalization(madat,madat.reverse,lowess.f,flaglevel) aplica a autonormalização na lâmina madat e na sua recíproca madat.reverse, com parâmetros opcionais flaglevel e lowess.f. Ela retorna a lâmina madat normalizada. Para obter a normalização da recíproca, basta inverter os argumentos madat e madat.reverse. Os comandos abaixo ilustram o uso dessa normalização, para a hibridização rat85 e sua recíproca rat89:
>rat89.self <- self.normalization(rat89.raw,rat85.raw)
Referências
[1] Yarger R. J., Reese G., King T. MySQL & mSQL. Cambridge: O'Reilly, 1999.[2] Wall L., Christiansen, T., Schwartz R. Programming Perl . Cambridge: O'Reilly, 1996.
[3] Yang Y.H., et al. Normalization for cDNA microarray data. Tech Report. January 2001.
[4] Finkelstein D.B., Gollub J, Cherry J.M. Normalization and systematic measurement error in cDNA microarray data. Journal of Computational and Graphical Statistics, v. 5, n.3, 2002.
*Projeto realizado com o apoio da Fundação de Amparo à Pesquisa do Estado de São Paulo.
