Análise de Dados Biológicos em R I. Operações básicas em R 1) Crie uma variável chamada x e atribua-­‐lhe o valor 3 (digitando x=3 ou x <-­‐3). Para consultar o valor de x, simplesmente digite x. Crie uma segunda variável chamada y e atribua-­‐lhe o valor 5. Realize operações aritméticas simples sobre x e y (adição, subtração, multiplicação, divisão, exponenciação -­‐ símbolo ^, divisão inteira %/%, resto de divisão inteira %%). Verifique que x e y são vectores de tamanho 1. 2) Defina o valor r = π/6 (o valor de π é obtido pela constante pi). Calcule as funções trigonométricas de r -­‐ coseno (função cos), seno (sin), tangente (tan). 3) Verifique a lista de variáveis criadas na sessão usando a função objects() ou ls(). 4) Crie uma variável do tipo string fazendo cor = “azul”. Verifique o seu valor. Crie uma outra variável cor2, atribuindo-­‐lhe o valor “amarelo”. a) Teste a função paste para criar uma variável que seja a junção de ambas as variáveis anteriores. Teste o efeito do parâmetro sep. b) Teste as funções pré-­‐definidas nchar e substr. Use a instrução help para investigar o que cada uma das funções faz. 5) Crie um vector v com os valores 3, 6, 9, 12, 15. Use o construtor v = c(3, 6, 9, 12, 15). Aceda ao 3º valor de v. a) Teste a função typeof. Faça typeof(v) para testar o tipo dos elementos de v. b) Defina o vector d com os valores de v divididos por 3. c) Defina o vector rq que contém as raízes quadradas dos valores de v (use a função sqrt). d) Crie o vector com todos os elementos de v menos o 3º. 6) Defina um vector p com os 20 primeiros números naturais (usando 1:20). a) Calcule o vector inverso de p (função rev). b) Calcule o vector resultante de somar a constante 3 a cada elemento de p. c) Calcule o vector 1/p. 7) Usando a função seq crie um vector com todos os números pares menores do que 20. Use examples(seq) para ver exemplos de utilização desta função. Selecione de p os elementos com índice par. 8) Crie um vector r de 100 números aleatórios de uma distribuição uniforme entre 0 e 1, usando a função runif. a) Experimente outros tipos de distribuições estatísticas (e.g. normal – função rnorm). b) Ordene o vector r com a função sort. Repita a operação colocando os valores por ordem decrescente. c) Verifique que r é um vector usando a função is.vector. d) Calcule o tamanho do vector r com a função length. e) Calcule algumas medidas estatísticas sobre r, tais como a média (função mean), a soma (sum), o maior valor (max), o menor valor (min), a mediana (median), a variância (var). 9) Um vector pode conter valores lógicos (booleanos). a) Crie o vector p com os números inteiros de 1 a 20. b) Defina vl = p > 12. Veja o conteúdo de vl e teste o seu tipo. c) Crie um vector q que tenha apenas os elementos de p maiores do que 12. 10) Os vectores podem ter nomes associados a cada posição (vectores com strings). Crie um vector que represente a tabela: João Pedro Rui Joana Manuel 12 8 15 11 a) Modifique a nota da Joana para 18. Arrays/ matrizes 1. Usando um comando matrix crie uma matriz mat com 2 linhas e 3 colunas, com os valores dos números inteiros de 1 a 6. Faça help(matrix) para saber como usar. a. Identifique a linha 2 da matriz (mat[2,]) e a coluna 3 (mat[,3]). b. Verifique as dimensões da matriz com a função dim. Se fizer v = dim(mat), de que tipo é o objecto v ? c. Calcule a soma dos elementos da matriz mat. d. Calcule a soma dos elementos da 1ª linha. e. Crie um vector cujos elementos sejam as médias dos valores de cada coluna. 2. Crie um array tri-­‐dimensional com as dimensões 2,5,3 e os trinta primeiros números inteiros. Use a função array. a. Calcule a soma dos elementos do array cujo valor na segunda dimensão é 3. b. Calcule a média dos valores que têm nas primeiras duas dimensões 1. c. Calcule o array que se obtém multiplicando todos os valores do array anterior por 2 e somando a constante 3 (as operações aritméticas funcionam nos arrays tal como nos vectores). 3. Usando as funções cbind e rbind (que permitem criar matrizes a partir da junção de colunas ou linhas provenientes de vectores ou outras matrizes), crie duas matrizes: a matriz A com dimensões 2 x 3 e com o valor 2 em todas as células; a matriz B com dimensões 2 x 3 e os valores entre 1 e 6. a. Calcule A+ B (com o operador aritmético +); b. Altere as dimensões da matrix B para que fique com 3x2. c. Calcule a multiplicação de A por B (usando o operador %*%); Listas, factores e data frames 1. Os objetos do tipo lista permitem definir estruturas de dados mais flexíveis pois permitem juntar vários objetos de vários tipos. Cada um destes objetos pode ser designado por campo e é acedido pelo símbolo $. a. Defina uma variável l como uma lista com dois campos: nome = “João” e idade = 30. Use a função list. b. Veja o conteúdo de l, aceda ao campo nome (l$nome) e ao campo idade. Verifique o valor de l[[1]] e l[[2]]. Tente perceber o comportamento do operador [[ ]]. Objetos do tipo factor São casos especiais de vectores usados para especificar classes (ou grupos) para os componentes de um vector. Vamos ver um exemplo. 2. Defina um vector onde estão representadas o grupo sanguíneo de 10 dadores, fazendo dadores = c(“A”,“B”,“O”,“AB”,“A”,“O”,“A”,“A”,“AB”,“A”). a. Crie a partir de dadores um factor usando a função factor: faça fdadores = factor(dadores) e consulte o conteúdo do objecto fdadores. b. Consulte a lista de classes com a função levels. c. Defina um vector p com os valores 72, 55, 75, 80, 60, 58, 53, 65, 72, 84 que se referem aos pesos dos dadores da amostra. d. Usando a função tapply(p, fra, mean) calcule a média dos pesos para cada uma das classes. Experimente diferentes funções como max, range ou sd). II. Data frames e pré-­‐processamento 1. Neste exercício será utilizado o data frame “CO2”. Este data frame regista o consumo de CO2 em plantas de 2 regiões distintas (Mississippi e Quebec) quando expostas a diferentes concentrações de CO2 e quando são sujeitas, ou não, a arrefecimento antes do início da experiência. a. Utilize o comando help para analisar a estrutura do data frame. Veja o conteúdo de cada um dos seus campos com o operador $ e verifique as dimensões do objecto. b. Qual o número de plantas sujeita a tratamento (campo “Treatment”) de arrefecimento (valor “chilled”)? c. Crie um sub-­‐conjunto dos dados que contenha apenas os dados referentes à região “Quebec”. d. Crie uma matriz com os dados de concentração de CO2 e de uptake. Atribua nomes às linhas e colunas, onde para as linhas é usada a identificação da planta e para as colunas as strings “concentr” e “consumo”. e. Considerando a matriz anterior, selecione todos os valores de concentração e consumo referentes à planta com identificação “Mc3”. f.
Crie um vector com os valores de consumo, a partir da matriz criada em d). g. Se tivesse um campo com 1000 plantas do tipo “QC3” sujeitas a uma concentração constante de 675 mL/L, qual o consumo de CO2? h. Em que condições o consumo de CO2 é o mais elevado? Como faria para obter a mesma informação mas para cada uma das regiões? i.
Qual a média de consumo de CO2 realizado pelas plantas da região “Mississipi” sujeitas a tratamento? 2. Considere a seguinte informação relacionada com as classificações de alunos numa disciplina do curso de Bioinformática: Nome TP1 TP2 Exame Extra 12345 Cátia Azevedo 15 17 17 1 32345 Paulo Cardoso 12 14 12 0 78532 Sandrina Alves 14 11 15 0,5 21456 Tiago Peixoto 18 19 17 0,5 a) Crie um data frame para representar a informação atrás apresentada. Note que os campos a negrito correspondem ao nome das linhas e colunas. b) Acrescente um novo atributo chamado “Class. Final”, onde a classificação é calculado utilizando a fórmula CF = 0,25*tp1 +0,25*tp2 +0,5 exame + extra
c) Crie uma função que dado um data frame , um valor de corte e campo (“tp” ou “ex”) retorne um novo data frame com a informação referente aos alunos que têm classificação superior ao valor de corte no exame (se o campo for “ex”) ou pelo menos num dos trabalhos práticos (se o campo for “tp”)
3. Considere a informação constante no ficheiro dados.csv e importe-­‐os para um data frame. a. Explore brevemente o data frame importado. b. Crie um novo data frame que contenha apenas os dados posteriores a 1990 referentes aos estados “Alabama”, “California”, “Florida” e “New York”. c. Referente aos dados obtidos na alínea b) crie um sumário dos dados onde apresente, para casa estado, o número total de crimes por ano. 4. Considere o data frame airquality, que contém dados sobre a qualidade do ar em New York entre Maio e Setembro de 1973. a. Visualize a estrutura do data frame utilizando o comando str. b. Verifique quantos valores dos campos Ozone e Solar.R não estão disponíveis. c. Calcule o valor médio da radiação solar durante o mês de Julho. d. Quais os índices das linhas que têm missing values (NA) no atributo Ozone? e. Crie um novo data frame onde todos os atributos têm valores válidos. f.
Para cada um dos meses, calcule a o máximo, mínimo, media, mediana, desvio padrão e variância dos valores de radiação solar. III. Gráficos 1. Tendo em conta o dataset CO2 atrás utilizado, crie um histograma (hist) para o campo uptake considerando intervalos de 5 umol/m^2 sec. Para calcular o limite superior considere limite = (round(max(CO2$uptake)/5)+1)*5. 2. Tendo em consideração os dado airquality, crie um diagrama de caixa (boxplot) para analisar a variação da Radiação solar entre os vários meses. Titulo: “Diagrama de caixas -­‐ Airquality”; Cor das caixas : verde; Texto dos eixos de X e Y : "Meses" e "Rad. Solar (inL angleys)"; Alterar nomes de meses para "Maio","Jun","Jul", "Agos", "Set". 3. Crie um gráficos de dispersão que represente a relação entre as variáveis Ozone e Solar.R. Deve filtrar os dados em falta. Titulo: “Ozono vs Radiação Solar”; Eixos: cor azul com números a laranja; Texto dos eixos : modificar cor e tamanho; Usar para o desenho dos pontos “pch =19” explorar outros símbolos;
Insira uma linha no ponto médio dos valores de Radiação Solar usando a função abline; 4. Crie dois vectores com 100 valores aleatórios de uma distribuição uniforme (runif). Utilizando a função qqplot, avalise se os vectores seguem a mesma distribuição. 5. Faça download da informação referente à classificação de 200 estudantes nas áreas: leitura, escrita e matemática hsb2 <-­‐ read.csv("http://www.ats.ucla.edu/stat/data/hsb2.csv") Crie uma matriz onde valor a_i,j indica se o aluno i tem aproveitamento superior a 60 na disciplina j. 6. Coloque na mesma janela, lado a lado, dois gráficos que representem as funções cos e sin no intervalo [–2 , 2]. IV. Inferência 1. Considere os dados do data frame “wines” do package kohonen. Este dataset consiste na análise química a 177 amostras de vinho cultivados na mesma região, provenientes de 3 castas. install.packages("kohonen")
data(wines, package = "kohonen")
wines
vintages
wine.classes
a. Importe os dados do dataset wines do package kohonen e explore o conteúdo das variáveis wine, vintages e wine.classes. Verifique o seu tipo e sumarie os seus dados. 2. Implemente em R testes estatísticos que permitam realizar as seguintes tarefas: a. Comparar os valores da coluna “alcohol” da matriz wines, entre os grupos “Barolo” e “Barbera”. Será que podemos concluir que as duas castas têm teores alcoólicos distintos com base neste teste ? b. Comparar os valores da coluna “alcohol” da matriz wines, entre os três grupos de castas. Realize a análise de variância deste teste. O que conclui ? c. Complete a análise da alínea anterior realizando testes entre cada par de grupos de castas. Compare os resultados que obteve com o boxplot. d. Verifique de que forma a variável o grau alcoólico (“alcohol”) pode variar em função da variabilidade das variáveis “malic acid”, “flavonoids” e “tot. phenols”. V. PCA e SVD 1. Neste exercício será utilizado de novo o data frame “wines” do package kohonen. a. Utilizando a função prcomp, efectue a análise dos principais componentes dos dados (normalizados). b. Apresente em forma de gráfico o resultado da função anterior e analise quantos componentes serão necessários considerar para explicar 90 % da variabilidade dos dados. c. Proceda ao mesmo tipo de análise utilizando a função princomp. Verifique as diferenças ao nível da estrutura de resultados e ao nível dos valores numéricos calculados (e.g. para os campos rotation / loadings) 2. Utilizando SVD calcule as matrizes de scores e loadings para o dataset anterior. (Obs: SVD consta da fatorização de uma matriz M (de dimensões n x m) em M = UDVT. A matrix de loadings corresponde à matriz V e a matriz de scores à matriz U*D. 3. Calcule a importância de cada um dos componentes sabendo que é dada pela fracção da variância explicada, isto é, FV = wines.svd$d^2 / sum(wines.svd$d^2), sendo wines.svd o resultado da aplicação da função svd aos dados. 4. Apresente três gráficos que representem a variância (wines.svd$d^2 / (nrow(wines) -­‐ 1)), variância relativa (FV ) e a percentagem acumulada da variância total usando a função cumsum. 5. Sabendo que a matriz de scores mostra a posição de cada amostra relativamente a cada PC, produza um gráfico que apresente a distribuição dos exemplos considerando os 2 principais componentes PC1 e PC2. Coloque símbolos e cores diferentes para cada classe. 6. Usando a função pairs, apresente a posição das amostras considerando os 4 principais componentes (use os atributos pch e col para atribuir símbolos e cores diferentes para cada classe). VI. Clustering 1. Neste grupo de exercícios vamos considerar um subgrupo das amostras presentes do dataset wines. Retire 30 amostras de forma aleatória. 2. Calcula as matrizes de distâncias, entre as várias amostras, considerando vários tipos de distância: euclidiana, manhattan e máxima (use o help da função dist). 3. Utilizando clustering hierárquico agrupe as amostras considerando as distâncias anteriormente calculadas. Teste a função com vários métodos ( "single", "complete" e "average"). Visualize os resultados de cada opção. Tente verificar até que ponto o clustering realizado está de acordo com as classes presentes na variável vintages. 4. Considere a função kmeans do R para agrupar os dados do dataset wines em 3 grupos. Verifique quantas amostras foram colocadas num grupo diferente do real (assumindo a classificação dada pela variável vintages). VII. Análise preditiva 1. Nesta ficha vamos usar dois conjuntos de dados, um para um problema de classificação e um outro para um problema de regressão. a. Carregue os dois conjuntos de dados: cats do package MASS e bodyfat do package mboost. b. Explore os conjuntos de dados realizando as seguintes tarefas: verifique se estes têm dados omissos (NAs), verifique os tipos de dados de cada campo, sumarie a distribuição de valores de cada variável. c. Visualize os dados verificando a distribuição de cada variável, visualizando as suas correlações. Use funções para criar boxplots, gráficos de dispersão, histogramas, onde for apropriado. Identifique possíveis outliers. 2. Para cada um dos conjuntos de dados, crie um conjunto de dados para treino (com 75% dos exemplos) e um conjunto de dados para teste (com os restantes 25%). Verifique a distribuição de exemplos de cada uma das classes com o comando table (no caso da classificação) e a distribuição da variável a prever no caso da regressão. Repita os comandos se a distribuição for muito desbalanceada. 3. Crie e visualize modelos baseados em árvores de decisão/ regressão para cada um dos conjuntos de dados anteriores, usando os dados de treino. Avalie o modelo nos exemplos de teste calculando métricas de erro adequadas em cada caso. Compare as árvores que obtém usando os packages party e rpart. 4. Use o método dos k vizinhos mais próximos no conjunto de dados cats. Calcule as medidas apropriadas para diferentes valores de k. 5. Use redes neuronais e support vector machines para criar modelos de previsão para o conjuntos de dados cats. 6. Crie um modelo linear para o conjunto de dados bodyfat. Calcule as medidas de erro adequadas. Visualize os valores previstos pelo modelo vs os valores reais. 7. Usando o método de naive bayes para criar modelos, faça a sua validação cruzada com 5 folds calculando a percentagem de exemplos corretamente classificados. Use o conjunto de dados cats. VIII. Análise de dados de expressão genética Esta sessão será dedicada a uma exploração de algumas das potencialidades dos packages integrados no Bioconductor, neste caso centrada no uso do package limma e na sua aplicação análise de dados de microarrays. 1. Carregar Dados Procure no site do ArrayExpress (http://www.ebi.ac.uk/arrayexpress/) a experiência “E-­‐GEOD-­‐25200” e faça download da raw data. Faça download dos ficheiros targets.txt e spottypes.txt do site: http://darwin.di.uminho.pt/aula_bioc Se não conseguir carregar os dados do Array Express tem os mesmos dados disponíveis no mesmo site (ficheiro gse25200.zip) a) Carregue o package Limma fazendo: library(limma)
b) Com o comando targets carregue a informação referente aos ficheiros Agilent Feature Extraction raw data e a informação da condição a que cada uma das amostras pertence. targets <- readTargets("targets.txt")
c) Usando a função read.maimages carregue os seus dados para um objecto RGList. Esta função permite carregar informação de microarrays dual channel, pelo que neste caso o campo R estará sem informação. RG <- read.maimages(targets, columns = list(G = "gMedianSignal", Gb = "gBGMedianSignal",
R = "gProcessedSignal", Rb = "gIsPosAndSignif"), annotation =
c("ProbeName","SystematicName"))
d) Substraia o background e aplique a normalização entre arrays. RG <- backgroundCorrect(RG, method="normexp", offset=1)
RG$G <- normalizeBetweenArrays(RG$G, method="quantile") e) Aplique o logaritmo às intensidades do canal Green. RG$G <- log2(RG$G)
2. Filtragem a) Leia a informação com os tipos de controlo do array. (Necessita do ficheiro spottypes.txt) spottypes = readSpotTypes(…) b) Adicione a informação do tipo a cada um dos spots do microarray. RG$genes$Status = controlStatus(spottypes,RG) c) Filtre os dados para apenas os que estão num spot do tipo “Gene”. i <- RG$genes$Status=="Gene"
RG.nocontrol = RG[i,]
d) Converta o objecto RG para um MAList. E <- new("MAList", list(targets=RG.nocontrol$targets, genes=RG.nocontrol$genes,
source=RG.nocontrol$source, M=RG.nocontrol$Gb, A=RG.nocontrol$G))
e) Use a função avereps para calcular a média dos spots correspondentes a um mesmo gene. E.avg <- avereps(E, ID=E$genes$ProbeName)
3. Carregar a anotação a) Faça load da biblioteca hgug4112a.db com anotação deste microarray. library(hgug4112a.db)
b) Associe a cada Probe o respectivo GENENAME e SYMBOL. IDs <- E.avg$genes$ProbeName
listaG <- mget(IDs, hgug4112aGENENAME)
listaS <- mget(IDs, hgug4112aSYMBOL)
E.avg$genes$Symbol <- unlist(listaS)
E.avg$genes$Description <- unlist(listaG)
4. Modelo linear a) Construa o design do modelo linear para as condições carregadas. f <- factor(targets$Condition, levels = unique(targets$Condition))
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)
b) Aplique a função lmfit às intensidades. fit <- lmFit(E.avg$A, design) 5. Comparações a) Crie a matriz de contrastes para as comparações que pretende fazer. contrast.matrix <- makeContrasts("CsA-Control","TGF-Control","TGF_CsA-Control",
levels=design)
b) Aplique a matriz de contrastes ao modelo de dados e calcule as estatísticas para estes dados. fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
6. Filtar resultados a) Para cada uma das condições e seguindo o exemplo crie uma tabela com os 20 genes mais significativos para cada uma das comparações. output <- topTable(fit2, adjust="BH", coef="CsA-Control", genelist=E.avg$genes,
number=20)
Expressão diferencial – miocroarray de 2 cores Background: “This experiment was conducted by the Erle Lab in UC San Francisco. This experiment aims to study the cell adhesion molecule integrin alpha4/beta7 which assists in directing the migration of blood lymphocytes to the intestine and associated lymphoid tissues. The goal of the study is to identify differentially expressed genes between the alpha4/beta7+ and alpha4/beta7-­‐ memory T helper cells. The study hypothesizes that differentially expressed genes may play a role in the adhesion or migration of T cells.” Após o download do ficheiro beta7.zip disponibilizado (no site referido acima): a) Carregue o ficheiro que contem os targets da experiência. > library(limma)
> targets <- readTargets("TargetBeta7.txt")
b) Ler os dados dos ficheiros, a função f filtra os spots que são consideradas como ‘bad’. >
>
>
>
f <- function(x) as.numeric(x$Flags > -75)
RG <- read.maimages(targets$FileName, source="genepix", wt.fun=f)
RG$printer <- getLayout(RG$genes)
RG$printer
c) Qual o tipo de dados de RG? Analise o conteúdo do objeto. d) Faça a correção de background, explorando 3 alternativas para o método de correção “substract” , “none” e “normexp” (com offset=25). Utilizando a função plotMA apresente os gráficos, de cada uma das amostras, para a correção de background “normexp”. >
>
>
#
>
>
>
RGsu <- …
RGno <- …
RGne <- …
Cria gráficos que relacionam as matrizes M e A para cada uma das correções
plotMA3by2(RGsu, prefix="MAsu")
plotMA3by2(RGno, prefix="MAno")
plotMA3by2(RGne, prefix="MAne")
e) Faça a normalização dos dados dentro do array e entre os arrays. MA <- normalizeWithinArrays(RGne)
MA <- normalizeBetweenArrays(MA)
f)
Compare a relação entre as matrizes M e A após a normalização, através da função plotMA, com os obtidos na alínea d. g) Construa a matriz modelo da experiência a utilizar na linearização do modelo. design <- modelMatrix(targets, ref="b7 -")
h) Utilizando as funções lmFit, eBayes e topTable apresente as 10 probes com maior expressão diferencial. Na função topTable usa o método de ajuste “fdr”. i)
Qual o nome dos genes que são diferencialmente expressos?