# IBI5086 - Lista 02
[Índice: Estatística IBI5086](/8z457TaXRmKsDwO5iHG6-w)
Feito pra disciplina IBI-5086, em Agosto/2020 (lista 01).
Alunos:
- Maria Fernanda Zaneli Campanari (nº USP: 12099257)
- Rafael Correia da Silva (nº USP: 9786869)
- Marília Viana Albuquerque de Almeida (nº USP: 12084806)
- Leonardo Yoshida (nº USP: 10399769)
- Gabriella Berwig Möller (nº USP: 11387062)
- Isabella Blagiem de Campos (nº USP: 12099240)
[toc]
## Problema 01.
### Respostas finais.
- Item a):
- P1: Não há evidência de interação, efeito principal tanto de A quanto de B, pois não há desvio do paralelismo
- P2: Não há evidência de interação, efeito principal tanto de A quanto de B, pois não há desvio do paralelismo
- P3: Há evidência de interação pois ocorre desvio do paralelismo;
- P4: Há evidência de interação pois ocorre desvio do paralelismo;
- P5: Não há evidência de interação, efeito principal tanto de A quanto de B, pois não há desvio do paralelismo
- P6: Não há evidência de interação, efeito principal tanto de A quanto de B, pois não há desvio do paralelismo
- P7: Há evidência de interação pois ocorre desvio do paralelismo;
### Item "a)" - perfis de médias
Primeiro, é necessário conceitualizar como será o gráfico de médias. Tomando como exemplo a tabela "P1", é visível que conforme a presença do gene dominante aumenta para qualquer um dos genes, também diminui o índice numérico, tanto para o gene A quanto para o gene B. Assim, o gráfico será similar ao seguinte, só que decrescente (Slide da Aula 3, p. 9):
### Gerando o gráfico para P1 no R.
Para gerar o gráfico de médias considerando os dois alelos, serão necessárias as informações gráficas:
- No eixo X; se o gene é "aa", "aA" ou "AA"
- No eixo Y; o valor atribuído;
- Colorir a reta em uma de três cores, considerando as opções "bb", "Bb" e "BB".
Assim, criei um arquivo `.txt` simples para dar entrada no R dos dados de P1:
```dataset
GenotypeSNP_A SNP_B Mean
aabb aa bb 20
aaBb aa Bb 18
aaBB aa BB 16
Aabb Aa bb 17
AaBb Aa Bb 15
AaBB Aa BB 13
AAbb AA bb 14
AABb AA Bb 12
AABB AA BB 10
```
Para construir o gráfico usando `ggplot2`:
```r=
# Bibliotecas
## GGplot: pretty plots
library(ggplot2)
## Problema 1.
### Gerando o gráfico para P1.
### Leitura da matriz
P1.Plot01.data <- read.delim(file = "./data/L02/P01_plots/P1.txt")
### Testes de plot
ggplot(data = P1.Plot01.data,
aes(x = SNP_A, y = Mean, color = SNP_B)) +
geom_point() + geom_line(aes(group = SNP_B), size=2)
```
E o resultado:
Considerando o paralelismo entre as médias por alelo, concluímos que não há evidência de interação para P1.
### Gerando todos os outros gráficos (P2 a P7).
Basicamente é a repetição dos gráficos acima, trocando apenas a matriz de entrada. Vou omitir o código, colocando apenas os gráficos e interpretação.
Em alguns dos gráficos, foi adicionado um modesto jitter (desvio) para que as linhas não se sobrepusessem.
Forma genérica do código R:
```r
### Método genérico
P1g.data <- read.delim(file = "./data/L02/P01_plots/P7.txt")
### Método genérico funcional pra plot com jitter
ggplot(data = P1g.data,
aes(x = SNP_A, y = Mean, color = SNP_B)) +
geom_point() + geom_line(aes(group = SNP_B),
position=position_jitter(w=0.1, h=0.1),
size=2)
### Método genérico funcional pra plot sem jitter
ggplot(data = P1g.data,
aes(x = SNP_A, y = Mean, color = SNP_B)) +
geom_point() + geom_line(aes(group = SNP_B),
position=position_jitter(w=0, h=0),
size=2)
```
### Gráfico P2
- Só há efeito (de ambas) quando em dupla homozigose;
### Gráfico P3
### Gráfico P4
### Gráfico P5
### Gráfico P6
### Gráfico P7
### P1. Item "b)" - gerando dados em DCA e interpretando
### Prêambulo
- É necessário usar os dados do enunciado para gerar os dados de entrada deste problema, com os critérios:
- Fatorial 3x3 (tal qual o item A);
- $\sigma^2$ = 1
- Usando como base as médias do item a)
- $r = 25$ (supondo $r$ da função `rnorm()`)
Dos parâmetros acima, entendo que para cada **valor** dos itens de P1 a P7 é necessário gerar 25 réplicas de mesma média, com variância 1. Um conjunto simples pode ser gerado com `rnorm()`, da seguinte forma:
```r
# Fórmula genérica:
# Para média 20, variância 1, 25 repetições:
rnorm(25, 20, 1)
```
Assim, parti iterativamente de cada tabela, gerando matrizes 9x28 para cada um dos genótipos, que foram testados todos contra todos. A seguir, as tabelas de anova.
### Interpretação das ANOVA
- P1, ação gênica intermediária: todos os alelos apresentaram efeito, quanto maior a presença do dominante menor a média da variável resposta;
- P2, dominância completa: só foram significativos os seguintes tratamentos: aaBB, AaBB, AAbb, AABb e AABB, sendo que o maior impacto fenotípico ($\beta = -10$) foi no genótipo AABB.
- P3, ação gênica complementar: houve significância exatamente na presença de qualquer um duplo dominante (AA ou BB ou ambos), com efeito muito similar (approx. $-9$ pra todos)
- P4, epistasia complexa: todos foram significativos com efeito muito similar, com o maior $\beta$ para AaBB e AABB
- P5, dominância parcial: todos foram significativos, de efeitos maiores nos casos de dominante ou duplo dominante (AaBB, AABB, AAbb e AABb)
- P6, superdominância: assim como acima, todos foram significativos, de efeitos maiores nos casos de dominante ou duplo dominante (AABB, AAbb e AABb)
- P7, ação gênica duplicada: somente houve efeito para o genótipo AABB, de $\beta = 10.3$
### Código fonte
```r
## Gerar com rnorm
rnorm(25, 20, 1)
set.seed(123)
# Todas as tabelas
list.files("./data/L02/P01_plots/")
# Entrada nas tabelas
P1.data <- read.delim(file = "./data/L02/P01_plots/P1.txt")
P2.data <- read.delim(file = "./data/L02/P01_plots/P2.txt")
P3.data <- read.delim(file = "./data/L02/P01_plots/P3.txt")
P4.data <- read.delim(file = "./data/L02/P01_plots/P4.txt")
P5.data <- read.delim(file = "./data/L02/P01_plots/P5.txt")
P6.data <- read.delim(file = "./data/L02/P01_plots/P6.txt")
P7.data <- read.delim(file = "./data/L02/P01_plots/P7.txt")
# Base data frame
P1.rnorm.mat <- data.frame(matrix(NA, nrow = 9, ncol = 28))
P2.rnorm.mat <- data.frame(matrix(NA, nrow = 9, ncol = 28))
P3.rnorm.mat <- data.frame(matrix(NA, nrow = 9, ncol = 28))
P4.rnorm.mat <- data.frame(matrix(NA, nrow = 9, ncol = 28))
P5.rnorm.mat <- data.frame(matrix(NA, nrow = 9, ncol = 28))
P6.rnorm.mat <- data.frame(matrix(NA, nrow = 9, ncol = 28))
P7.rnorm.mat <- data.frame(matrix(NA, nrow = 9, ncol = 28))
# Gera os três genótipos
P1.rnorm.mat[,c(1:3)] <- P1.data[,c(1:3)]
P2.rnorm.mat[,c(1:3)] <- P1.data[,c(1:3)]
P3.rnorm.mat[,c(1:3)] <- P1.data[,c(1:3)]
P4.rnorm.mat[,c(1:3)] <- P1.data[,c(1:3)]
P5.rnorm.mat[,c(1:3)] <- P1.data[,c(1:3)]
P6.rnorm.mat[,c(1:3)] <- P1.data[,c(1:3)]
P7.rnorm.mat[,c(1:3)] <- P1.data[,c(1:3)]
# Preenche as matrizes com as médias
for (genotype in P7.data[,1]) {
# Pega o número da linha
print(genotype)
genotype_no = grep(pattern = genotype, x = P7.data$Genotype)
## Gera os valores pra esse genótipo
new_values = rnorm(n = 25, mean = P7.data[genotype_no, 4], sd = 1)
## Guarda os novos valores
P7.rnorm.mat[genotype_no, 4:28] = new_values
}
# Nomes de colunas
colnames(P1.rnorm.mat) = c(colnames(P1.data)[1:3], paste("Mean", 1:25, sep="_"))
colnames(P2.rnorm.mat) = c(colnames(P1.data)[1:3], paste("Mean", 1:25, sep="_"))
colnames(P3.rnorm.mat) = c(colnames(P1.data)[1:3], paste("Mean", 1:25, sep="_"))
colnames(P4.rnorm.mat) = c(colnames(P1.data)[1:3], paste("Mean", 1:25, sep="_"))
colnames(P5.rnorm.mat) = c(colnames(P1.data)[1:3], paste("Mean", 1:25, sep="_"))
colnames(P6.rnorm.mat) = c(colnames(P1.data)[1:3], paste("Mean", 1:25, sep="_"))
colnames(P7.rnorm.mat) = c(colnames(P1.data)[1:3], paste("Mean", 1:25, sep="_"))
# Reformata os dados
P1.new <- melt(P1.rnorm.mat)
P2.new <- melt(P2.rnorm.mat)
P3.new <- melt(P3.rnorm.mat)
P4.new <- melt(P4.rnorm.mat)
P5.new <- melt(P5.rnorm.mat)
P6.new <- melt(P6.rnorm.mat)
P7.new <- melt(P7.rnorm.mat)
# Executa as anovas
P1.aov <- summary(lm(value ~ as.factor(Genotype), data = P1.new))
P2.aov <- summary(lm(value ~ as.factor(Genotype), data = P2.new))
P3.aov <- summary(lm(value ~ as.factor(Genotype), data = P3.new))
P4.aov <- summary(lm(value ~ as.factor(Genotype), data = P4.new))
P5.aov <- summary(lm(value ~ as.factor(Genotype), data = P5.new))
P6.aov <- summary(lm(value ~ as.factor(Genotype), data = P6.new))
P7.aov <- summary(lm(value ~ as.factor(Genotype), data = P7.new))
```
### Anova P1 - Ação gênica intermediária
### Anova P2 - Dominância completa
### Anova P3 - Ação gênica complementar
### Anova P4 - Epistasia complexa
### Anova P5 - Dominância parcial
### Anova P6 - Superdominância
### Anova P7 - Ação gênica duplicada
### Gerando os dados por efeito "P" com base nas médias das tabelas
## Problema 02.
### Item "a)" - Avaliar o efeito do SNP
Plano:
1. [x] Importar a matriz no R e entender os dados
2. [x] Escolher uma variável (escolhida: PESO)
3. [x] Testar a variável para as premissas da anova, avaliar se é necessária a transformação
1. [x] Shapiro-Wilk (distribuição normal)
2. [x] Homoscedasticidade (distribuição igual das variâncias)
Para a variável peso, foi impossível atender as condições para ANOVA paramétrica, sendo necessário lançar mão de uma análise não-paramétrica.
Primeiro, testamos a distribuição normal e homogeneidade da variâncias:
```r=
## Ler os dados
lista02_dat <- read.delim(
file = "./data/L02/map.txt",
header = TRUE,
sep = "\t"
)
## Teste de modelo para o segundo SNP (rs7289830_T)
as.factor(lista02_dat$rs7289830_T)
test.fitSNP <- lm(Peso ~ as.factor(rs7289830_T),
data = lista02_dat)
aov.test.fitSNP <- anova(test.fitSNP)
summary(test.fitSNP)
colnames(lista02_dat)
## Passo 3. Teste da variável
### Teste global
### Aqui é testada toda a população
shapiro.test(lista02_dat$Peso)
### Porém, quero ver só uma parte, se algum subconjunto é normal
test.Peso.SNP1 <- c(lista02_dat[, "rs7289830_T"] == 0)
test.Aa <- lista02_dat[test.Peso.SNP1,]
### Depois de separado o subconjunto, testo novamente
hist(test.Aa$Peso)
### Teste de variâncias (homoscedasticidade)
leveneTest(lista02_dat$Peso, as.factor(lista02_dat$rs4410381_A))
```
Em seguida, tentamos transformar os dados para que atendessem a normalidade, mesmo após diversas transformações ainda era impossível aceitar como normal os dados pelo teste de Shapiro:
```r=
# Testes de transformações pra variável PESO ----
## 1. Potência
shapiro.test(lista02_dat$Peso**2)
## 2. Log
shapiro.test(log10(lista02_dat$Peso))
shapiro.test(log(lista02_dat$Peso))
## 3. Raiz-quadrada e cúbica
shapiro.test(sqrt(lista02_dat$Peso))
shapiro.test(lista02_dat$Peso**1/3)
## 4. Box-cox
### Remover valores perdidos
lista02_semNAs <- lista02_dat[!is.na(lista02_dat$Peso), ]
lambda_Peso <- BoxCoxTrans(lista02_semNAs$Peso)
shapiro.test(predict(lambda_Peso, lista02_semNAs$Peso))
## Testes de outlier
### Tentar podar umas pontas pra "normalizar" os dados
plotNormalHistogram(lista02_semNAs$Peso)
vecPeso_Outliers <- !(lista02_semNAs$Peso > 130000 | lista02_semNAs$Peso < 40000)
lista02_outliers <- lista02_dat[vecPeso_Outliers, ]
### Transformação
shapiro.test(lista02_outliers$Peso)
lista02_outliers_semNAs <- lista02_outliers[!is.na(lista02_outliers$Peso), ]
lambda_Peso_Outliers <- BoxCoxTrans(lista02_outliers_semNAs$Peso)
shapiro.test(predict(lambda_Peso_Outliers, lista02_outliers_semNAs$Peso))
```
Considerando que não foi possível chegar a normalidade em nenhum caso, resolvemos assim tentar trocar de variável antes de lançar mão de estatística não-paramétrica. Testando para outras variáveis, quase todas deram distribuição não normal, exceto algumas poucas:
```r=
## Função transformVars
## Aplica todas as transformações possíveis a um dado vetor
## Retorna todos os p-valores de shapiro e levene
transformVars <- function(vector, factor=NA) {
# Bibliotecas utilizadas
require(emmeans)
require(readxl)
require(caret)
require(lmtest)
require(car)
require(rcompanion)
require(userfriendlyscience)
clean_vector <- vector[!is.na(vector)]
#print(clean_vector)
## Primeiro shapiro (sem modificações)
p00.clean <- shapiro.test(vector)$p.value
#print(p.clean.shapiro)
## Todas as transformações
### Ao quadrado
transform01.powr <- clean_vector**2
### Raiz quadrada
transform02.sqrt <- sqrt(clean_vector)
### Raiz cúbica
transform03.cbrt <- clean_vector**1/3
### Log 10
transform04.logt <- log10(clean_vector)
### Log NP
transform05.logn <- log1p(clean_vector)
### Box cox
transform06.boxc <- predict(BoxCoxTrans(clean_vector), clean_vector)
## Todos os shapiro
p01.power <- shapiro.test(transform01.powr)$p.value
p02.power <- shapiro.test(transform02.sqrt)$p.value
p03.power <- shapiro.test(transform03.cbrt)$p.value
p04.power <- shapiro.test(transform04.logt)$p.value
p05.power <- shapiro.test(transform05.logn)$p.value
p06.power <- shapiro.test(transform06.boxc)$p.value
## Salva os resultados
out.vector <- c(
"00.Raw.Values" = p00.clean,
"01.Powr.Transform" = p01.power,
"02.Sqrt.Transform" = p02.power,
"03.Cbrt.Transform" = p03.power,
"04.Logt.Transform" = p04.power,
"05.Logn.Transform" = p05.power,
"06.Boxc.Transform" = p06.power
)
## Dá a saída
print(out.vector)
}
## Matriz para os valores de shapiro
shapiro.all <- as.data.frame(cbind(
transformVars(lista02_dat[,02]),
transformVars(lista02_dat[,03]),
transformVars(lista02_dat[,04]),
transformVars(lista02_dat[,05]),
transformVars(lista02_dat[,06]),
transformVars(lista02_dat[,07]),
transformVars(lista02_dat[,08]),
transformVars(lista02_dat[,09]),
transformVars(lista02_dat[,10]),
transformVars(lista02_dat[,11]),
transformVars(lista02_dat[,12]),
transformVars(lista02_dat[,13])
))
## Atribui nomes de colunas
colnames(shapiro.all) <- colnames(lista02_dat[,c(2:13)])
```
E o resultado das transformações:
Entendemos que somente quanto o p-value é superior a 0.05 no teste de Shapiro é que podemos aceitar que a distribuição dos dados é normal; sendo assim prosseguiremos com a ANOVA usando a variável quantitativa CIRCABD.
Geramos, por fim, os coeficientes para todos os SNPs:
| Parameter | Estimate | Std. Error | t value | Pr(>\|t\|) |
|------------------------|----------|------------|----------|------------|
| Intercept_rs12628452_A | 87.28356 | 0.432821 | 201.6622 | 0 |
| Intercept_rs7289830_T | 87.21188 | 0.489369 | 178.2129 | 0 |
| Intercept_rs5746356_A | 87.15909 | 0.494296 | 176.3298 | 0 |
| Intercept_rs10154759_C | 87.23293 | 0.426257 | 204.6488 | 0 |
| Intercept_rs7354790_T | 87.46502 | 0.793841 | 110.1795 | 0 |
| Intercept_rs6423472_C | 87.58436 | 0.783981 | 111.7175 | 0 |
| Intercept_rs1041770_A | 86.66369 | 0.515957 | 167.967 | 0 |
| Intercept_rs9712893_G | 87.33205 | 0.436834 | 199.9207 | 0 |
| Intercept_rs6010318_A | 87.27397 | 0.431481 | 202.2662 | 0 |
| Intercept_rs7285246_T | 87.86282 | 0.518755 | 169.3725 | 0 |
| Intercept_rs8138488_G | 88.11679 | 0.600958 | 146.6273 | 0 |
| Intercept_rs9617160_T | 88.17009 | 0.529458 | 166.5289 | 0 |
| Intercept_rs2027649_G | 88.0368 | 0.486097 | 181.1094 | 0 |
| Intercept_rs6010418_A | 87.46216 | 0.450113 | 194.3115 | 0 |
| Intercept_rs1987558_T | 87.35 | 0.497113 | 175.7144 | 0 |
| Intercept_rs140378_G | 87.35624 | 0.457376 | 190.9944 | 0 |
| Intercept_rs131560_C | 87.37286 | 0.462826 | 188.7813 | 0 |
| Intercept_rs5748616_C | 87.35608 | 0.665194 | 131.3242 | 0 |
| Intercept_rs5994019_G | 87.33988 | 0.427905 | 204.1104 | 0 |
| Intercept_rs5748662_T | 87.0992 | 0.632583 | 137.6882 | 0 |
| Intercept_rs5994034_T | 87.35304 | 0.495371 | 176.3387 | 0 |
| Intercept_rs2212272_C | 87.06048 | 0.549497 | 158.4368 | 0 |
| Intercept_rs4010554_A | 87.42105 | 0.682545 | 128.081 | 0 |
| Intercept_rs4010560_C | 87.45596 | 0.438016 | 199.6639 | 0 |
| Intercept_rs4010550_G | 87.07796 | 0.631994 | 137.7829 | 0 |
| Intercept_rs11089179_G | 87.24806 | 0.619942 | 140.7358 | 0 |
| Intercept_rs9680776_A | 87.10484 | 0.777675 | 112.0067 | 0 |
| Intercept_rs2379965_C | 87.3035 | 0.761642 | 114.6253 | 0 |
| Intercept_rs2379981_G | 87.3427 | 0.446788 | 195.4904 | 0 |
| Intercept_rs4535153_C | 87.3651 | 0.447744 | 195.1228 | 0 |
| Intercept_rs5747940_A | 88.01754 | 0.808878 | 108.8144 | 0 |
| Intercept_rs5746647_G | 87.61393 | 0.464353 | 188.6794 | 0 |
| Intercept_rs16980739_A | 87.70562 | 0.46973 | 186.7152 | 0 |
| Intercept_rs9605923_T | 86.97955 | 0.526608 | 165.1696 | 0 |
| Intercept_rs5747999_G | 86.97397 | 0.569766 | 152.6486 | 0 |
| Intercept_rs11089263_C | 86.67931 | 0.716654 | 120.95 | 0 |
| Intercept_rs11089264_G | 86.77346 | 0.694891 | 124.8735 | 0 |
| Level01_rs2027649_G | -3.11218 | 0.989146 | -3.14633 | 0.001712 |
| Level01_rs9617160_T | -2.66791 | 0.967076 | -2.75874 | 0.005939 |
| Level01_rs8138488_G | -2.03297 | 0.888902 | -2.28706 | 0.022443 |
| Level02_rs1041770_A | 11.16965 | 5.006818 | 2.230888 | 0.025963 |
| Level01_rs7285246_T | -1.80448 | 0.943554 | -1.91243 | 0.056172 |
| Level02_rs4010560_C | 15.54404 | 8.616787 | 1.803925 | 0.07161 |
| Level01_rs16980739_A | -2.03186 | 1.130705 | -1.79699 | 0.072704 |
| Level02_rs11089263_C | 2.181617 | 1.224731 | 1.781303 | 0.075232 |
| Level01_rs1041770_A | 1.601621 | 0.93467 | 1.713569 | 0.086992 |
| Level02_rs9712893_G | 20.66795 | 12.20792 | 1.692996 | 0.090832 |
| Level01_rs5746647_G | -1.9205 | 1.140193 | -1.68437 | 0.092487 |
| Level01_rs4010560_C | -2.96616 | 1.79293 | -1.65437 | 0.098436 |
| Level02_rs11089264_G | 1.897966 | 1.244445 | 1.525151 | 0.127612 |
| Level01_rs10154759_C | 5.267073 | 3.549297 | 1.483976 | 0.138195 |
| Level02_rs5748662_T | 2.018451 | 1.468388 | 1.374604 | 0.169626 |
| Level01_rs5747940_A | -1.38454 | 1.010799 | -1.36975 | 0.171139 |
| Level02_rs2379981_G | 16.6573 | 12.21946 | 1.363178 | 0.173196 |
| Level01_rs6010318_A | 3.788527 | 3.087047 | 1.227233 | 0.220088 |
| Level02_rs7285246_T | -3.14853 | 2.714479 | -1.1599 | 0.246429 |
| Level02_rs5747999_G | 1.701706 | 1.531996 | 1.110777 | 0.266987 |
| Level02_rs5748616_C | 1.560584 | 1.412722 | 1.104664 | 0.269626 |
| Level02_rs9605923_T | 2.674292 | 2.452675 | 1.090358 | 0.275872 |
| Level01_rs6010418_A | -1.49157 | 1.551576 | -0.96133 | 0.336675 |
| Level02_rs2212272_C | 2.166789 | 2.666359 | 0.81264 | 0.416664 |
| Level02_rs131560_C | -5.70619 | 7.084915 | -0.8054 | 0.420822 |
| Level02_rs140378_G | -5.68957 | 7.065927 | -0.80521 | 0.420929 |
| Level02_rs7354790_T | 1.007029 | 1.257509 | 0.800813 | 0.423489 |
| Level01_rs9605923_T | 0.726766 | 0.912111 | 0.796795 | 0.425798 |
| Level01_rs1987558_T | -0.77347 | 1.001807 | -0.77207 | 0.440294 |
| Level01_rs6423472_C | -0.7415 | 0.985003 | -0.75279 | 0.451789 |
| Level02_rs16980739_A | -2.90562 | 3.89054 | -0.74684 | 0.455372 |
| Level02_rs1987558_T | 1.91087 | 2.587233 | 0.738577 | 0.460377 |
| Level02_rs9680776_A | 0.918599 | 1.332868 | 0.68919 | 0.490898 |
| Level01_rs5746356_A | 0.699605 | 1.030678 | 0.678781 | 0.497468 |
| Level02_rs5746647_G | 3.386067 | 4.997646 | 0.677532 | 0.498257 |
| Level01_rs5994019_G | -1.95099 | 2.910943 | -0.67023 | 0.5029 |
| Level01_rs2212272_C | 0.591393 | 0.901716 | 0.655853 | 0.512105 |
| Level02_rs4535153_C | 5.634899 | 8.65317 | 0.651195 | 0.515101 |
| Level02_rs4010550_G | 1.152812 | 1.804655 | 0.638799 | 0.523145 |
| Level01_rs4535153_C | -0.78371 | 1.391814 | -0.56308 | 0.573531 |
| Level01_rs7354790_T | -0.51185 | 1.025691 | -0.49903 | 0.6179 |
| Level02_rs9617160_T | -1.59866 | 3.315537 | -0.48217 | 0.629819 |
| Level01_rs5747999_G | 0.43959 | 0.912108 | 0.481949 | 0.62997 |
| Level01_rs11089263_C | 0.453849 | 0.949987 | 0.477742 | 0.632961 |
| Level02_rs4010554_A | 0.668056 | 1.398471 | 0.477705 | 0.632987 |
| Level02_rs6423472_C | 0.589551 | 1.24189 | 0.474721 | 0.635112 |
| Level01_rs9712893_G | -0.81205 | 1.779795 | -0.45626 | 0.648322 |
| Level01_rs2379981_G | -0.63682 | 1.397827 | -0.45558 | 0.648811 |
| Level02_rs2027649_G | 1.9632 | 4.323943 | 0.45403 | 0.649926 |
| Level01_rs5748616_C | -0.40646 | 0.904484 | -0.44938 | 0.653272 |
| Level01_rs4010554_A | -0.39611 | 0.917125 | -0.43191 | 0.66592 |
| Level02_rs6010318_A | -5.27397 | 12.23459 | -0.43107 | 0.666531 |
| Level01_rs5994034_T | -0.39673 | 0.985316 | -0.40264 | 0.687318 |
| Level01_rs11089264_G | 0.372142 | 0.944873 | 0.393854 | 0.693793 |
| Level01_rs131560_C | -0.34786 | 1.209857 | -0.28752 | 0.773787 |
| Level02_rs5747940_A | -0.30252 | 1.194663 | -0.25322 | 0.800158 |
| Level01_rs2379965_C | -0.21758 | 0.967425 | -0.22491 | 0.822108 |
| Level02_rs11089179_G | -0.29154 | 1.593709 | -0.18293 | 0.854896 |
| Level02_rs8138488_G | 0.277949 | 1.521253 | 0.18271 | 0.85507 |
| Level02_rs5994034_T | 0.646962 | 3.563572 | 0.181549 | 0.855981 |
| Level01_rs12628452_A | 0.327548 | 2.081993 | 0.157324 | 0.875028 |
| Level02_rs2379965_C | 0.210012 | 1.386799 | 0.151436 | 0.879671 |
| Level02_rs7289830_T | 0.399233 | 2.920312 | 0.136709 | 0.891294 |
| Level01_rs7289830_T | 0.121455 | 1.027019 | 0.11826 | 0.905891 |
| Level01_rs9680776_A | 0.108495 | 0.968545 | 0.112018 | 0.910836 |
| Level02_rs5746356_A | -0.20076 | 2.552533 | -0.07865 | 0.93733 |
| Level01_rs4010550_G | -0.06052 | 0.91178 | -0.06637 | 0.9471 |
| Level01_rs140378_G | -0.06421 | 1.236585 | -0.05192 | 0.958603 |
| Level01_rs11089179_G | -0.02644 | 0.886744 | -0.02982 | 0.97622 |
| Level01_rs5748662_T | -0.00294 | 0.894009 | -0.00329 | 0.997378 |
#### Resposta final
Foram significativos a 0,05 os 04 níveis que explicam em parte CIRCABD:
- Level01_rs2027649_G (Heterozigoto)
- Level01_rs9617160_T (Heterozigoto)
- Level01_rs8138488_G (Heterozigoto)
- Level02_rs1041770_A (Homozigoto dominante)
[TO DO] Obs: Lucas deu a dica para testar os resíduos; incluir ao final.
Código para gerar todas ANOVA:
```r=
## Somente os valores não NA com BoxCox
lista02_boxcox <- lista02_dat[!is.na(lista02_dat$CIRCABD),]
## Laço de fato
columnCounter = 14
## Matriz em branco pra começar
template.df <- data.frame()
for (column in lista02_boxcox[,c(14:50)]) {
## Usa cada coluna como fator
fitSNP <- lm(CIRCABD ~ as.factor(column),
data = lista02_boxcox)
## Guarda o sumário
fitSummary <- summary(fitSNP)
## Extrai os coeficientes do sumário
summaryData <- fitSummary$coefficients
## Extrai o nome do SNP e guarda na coluna do sumário
row.names(summaryData)[1] <- paste(
"Intercept", colnames(lista02_boxcox)[columnCounter],
sep = "_"
)
## Para o primeiro fator em relação à casela
row.names(summaryData)[2] <- paste(
"Level01", colnames(lista02_boxcox)[columnCounter],
sep = "_"
)
## Para o segundo fator em relação à casela, se existir
if (! is.na(row.names(summaryData)[3])) {
row.names(summaryData)[3] <- paste(
"Level02", colnames(lista02_boxcox)[columnCounter],
sep = "_"
)}
## Une os dados da nova ANOVA com a matriz geral
template.df <- rbind(template.df,
summaryData)
## Conta as colunas
columnCounter = columnCounter + 1
}
```
### Item "b)" - Manhattan plot
- Basta criar um gráfico de dispersão (scatter plot), em que o eixo X são cada SNP, e o eixo Y são os valores de -log10(p-value).
```r
### Dados a serem usados:
### Índice de p-valores:
itemb.df <- template.df[template.df$`Pr(>|t|)` > 0, ]
colnames(itemb.df)[4] <- "p_value"
## Salva o -log(p)
itemb.df$mlogp <- -log(itemb.df$p_value)
## Nossos cromossomos são os nomes dos SNP
clean_SNP_names <- gsub(pattern = "Level0._",
replacement = "",
x = row.names(itemb.df),
perl = F)
## Salva os nomes limpos dos SNPs
itemb.df$SNP <- clean_SNP_names
## Dispersão básico
require(ggplot2)
ggplot(itemb.df, aes(x=SNP, y=mlogp, color = mlogp)) +
geom_point() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
geom_hline(yintercept = 2.995, linetype="dashed", color="red") +
xlab("SNP Index") +
ylab("-log(p-value") +
labs(title = "Manhattan plot of SNPs")
```
### Item "c)" - Volcano plot do efeito aditivo
```r
map <- read.table("/home/leonardo/Seagate/Doutorado/Bioinformatica/Julia/Topico_2/map.txt", header = TRUE)
aditivo <- data.frame( map[-1:-13] )
for( i in 1:ncol(aditivo) ){
aditivo[ , i] = as.factor(aditivo [ , i])
}
#Recodificação para modelo que diretamente permite estimar
#efeitos Aditivo e de Dominância
library(car)
aditivo2 <- aditivo
dominante <- aditivo
for( i in 1:ncol(aditivo) ){
aditivo2[ , i] <- Recode(aditivo[ , i], "'0'=(-1);'1'=0;'2'=1",as.factor=FALSE)
dominante[ , i] <- Recode(aditivo[ , i], "'0'=0;'1'=1;'2'=0", as.factor=FALSE)
}
#LM
for(i in 1:ncol(aditivo2) ){
assign( paste0("lm",i),lm(map$GLICOSE ~ aditivo2[ ,i] + dominante[ ,i]) )
}
#P VALORES e XA*2
pvalores=c()
xa_2 = c()
for(i in 1:ncol(aditivo2)){
pvalor = summary( get( paste0("lm",i) ) ) $ coefficients[2,4]
pvalores = c(pvalores,pvalor)
xa_2i = 2*summary( get( paste0("lm",i) ) ) $ coefficients[2,1]
xa_2 = c(xa_2, xa_2i)
}
ex2_letrac = data.frame( cbind(xa_2,-log10(pvalores)) )
colnames(ex2_letrac) = c("xa2","logP")
library(ggplot2)
ggplot(ex2_letrac, aes(x=xa2, y=logP) ) +
geom_point() +
geom_hline(yintercept=-log10(0.05), col="red") +
theme_bw()
#Fim do codigo que eu usei (Leonardo)
```
### Item "d)" - Efeito de dominância em Volcano plot
Plot (Leonardo):
```
#EXERCICIO 2 LETRA D
pvalores_xd = c()
xd = c()
#É a mesma lógica da letra C, só foi necessário um loop
#pra parsear o p-valor de xd, e evitar um erro de OOB
for(i in 1:ncol(aditivo2)){
xd_i = get( paste0("lm",i) ) $ coefficients[3]
xd = c(xd, xd_i)
if (length(summary( get( paste0("lm",i) ) )$coefficients[,4]) == 3) {
pvalor_xd = summary( get( paste0("lm",i) ) )$coefficients[3,4]
pvalores_xd = c(pvalores_xd,pvalor_xd)
} else {
pvalor_xd = NA
pvalores_xd = c(pvalores_xd,pvalor_xd)
}
}
#União de colunas e verificação dos NAs nas mesmas linhas
dominante <- data.frame( cbind (-log10(pvalores_xd), xd))
colnames(dominante) <- c("logP", "xd")
which(is.na(dominante$logP))
which(is.na(dominante$xd))
#REMOCAO DE NA
dominante <- na.omit(dominante)
which(is.na(dominante$logP))
which(is.na(dominante$xd))
library(ggplot2)
ggplot(dominante, aes(x=xd, y=logP) ) +
geom_point() +
geom_hline(yintercept=-log10(0.05), col="red") +
theme_bw()
```
Gráfico de Dominância para os dados de glicose
### Item "e)" - Métodos meta-analíticos
### Resposta 2e)
- Valor p-global foi de 0.005050.
#### Plano de análise
1. Concatenar os válores da fórmula por estudo em uma única matriz
2. Operar e interpretar os Z-score e por fim o P-valor global
#### Variáveis
- $N_i$: número de pacientes
- $P_i$: p-valor de significância para efeito aditivo?
- $\Delta_i$: se o efeito de cada SNP deve ser interpretado como aditivo ou subtrativo; sinal do "fold-change"
- $\Phi$: número de SNPs com efeito não-zero, ou número de significativos.
#### Resolução no R
```r
#### Resolução 2e) ----
### Baseado em tamanho amostral
### Começamos dos dados brutos
map <- read.delim(
file = "./data/L02/map.txt",
header = TRUE,
sep = "\t"
)
## Matriz para salvar resultados
meta_out_df <- data.frame(matrix(ncol = 6, nrow = 1), stringsAsFactors = F)
## Nomes
colnames(meta_out_df) <- c(
"Treatment", "N_i", "P_i", "Delta_i", "Z_i", "W_i"
)
## Contador de coluna
column_ct <- 2
## Contador do Phi global (todos os significativos)
Phi_global = 0
## Loop para extrair os valores para cada estudo todos os SNPs
## Início do For1
for (treatment in map[,2:10]) {
#print(head(treatment))
## Limpeza da tabela
### Deixa só os casos completos
tmp.treatment <- cbind(
treatment,
map[,14:115]
)
## Ao final deve haver uma matriz:
# 10 linhas (1 por tratamento)
# 6 colunas:
# [ok] N_i
# [OK] P_i
# [Ok] Delta_i
# [Ok] Phi_i
# Z_i
# Wi
## Efeito sempre é positivo, isto é, o SNP irá adicionar numericamente
Delta_i = 1
## Qual é o nivel de significancia?
P_i = 0.05
## Assim, o loop irá obter sistematicamente,
## para cada coluna/tratamento os valores.
## Variável p_loop guarda os p_valores
p_loop_tmp <- c()
## Para cada SNP, conte quantos p-valores são significativos para
## Efeito de aditividade
for (snp in map[,14:115]) {
# print(snp)
ef_ad <- Recode(snp, "'0'=(-1);'1'=0;'2'=1", as.factor=FALSE)
ef_do <- Recode(snp, "'0'=(+0);'1'=1;'2'=0", as.factor=FALSE)
## Extrai os p-valores do aditivo e dominante
p_ad <- summary(lm(
treatment ~ ef_ad + ef_do, na.action = na.exclude))$coefficients[2,4]
## Guarda o p-valor
p_loop_tmp <- c(p_loop_tmp, p_ad)
# Fim do for2
}
## Quantos p-valores existem?
N_i = length(p_loop_tmp)
## Quantos foram significativos?
Phi_i = sum(p_loop_tmp < 0.05, na.rm = T)/N_i
## Salva os significativos no Phi_global
Phi_global = c(Phi_global, Phi_i)
## Calcular o Z_i
Z_i = (1/Phi_i)*(P_i/2)*(Delta_i)
## Calcular o W_i
W_i = sqrt(N_i)
## Salva os resultados
treat.name <- colnames(map)[column_ct]
out.vector <- c(treat.name, N_i, P_i, Delta_i, Z_i, W_i)
meta_out_df[column_ct-1, ] <- out.vector
colnames(meta_out_df) <- c(
"Treatment", "N_i", "P_i", "Delta_i", "Z_i", "W_i"
)
## Avança o nome de coluna
column_ct = column_ct + 1
# Fim do For1
}
# Matriz de saida
meta_out_df
# Calculo do Z global
Z_global_numerador <- sum(
as.double(meta_out_df$Z_i)*as.double(meta_out_df$W_i))
Z_global_denominador <- sqrt(sum(as.double(meta_out_df$W_i)**2))
Z_global = Z_global_numerador/Z_global_denominador
# Quantos SNPs tem efeito no total?
Phi_effect <- max(Phi_global)
# Calculo do P global
P_global = 2*Phi_effect*abs(-Z_global)
```
### Item "f)" - Interação de SNP
- 05 SNPs interagiram com o sexo em nosso modelo ANOVA: "rs2212272_C" "rs2845362_G" "rs4819884_T" "rs5746892_G" "rs5748636_A"
- Código em R:
```r
# Item f) ------
## Ver algumas interações
scatterplot(CIRCABD ~ rs12628452_A | SEX, data=map)
scatterplot(CIRCABD ~ rs10154759_C | SEX, data=map)
### Interação entre sexo e SNP para variável CIRCABD
fit.test <- lm(
CIRCABD ~ rs12628452_A + as.factor(SEX):rs12628452_A, dat=map)
summ <- summary(fit.test)
## Vars
keep <- data.frame(matrix(nrow=1, ncol=2), stringsAsFactors = F)
### Loop para tirar todas as interações
#### Para cada nome de coluna (snp)
for (coln in colnames(map[,14:115])) {
# Descubra qual é o número
colnumber <- grep(x = colnames(map), coln)
# Extraia os SNPs daquele numero
snp <- map[,colnumber]
# Faz o modelo estrutural (lm)
fit <- lm(CIRCABD ~ snp + as.factor(SEX):snp, data = map)
# Extrai os valores do modelo estrutural
summary <- summary(fit)
coef.interaction <- summary$coefficients[3,4]
# Salva os valores
out.vector <- c(coln, coef.interaction)
# Colar a coluna na matriz de saida
keep <- rbind(keep, out.vector)
}
# Remove os NAs
keep <- keep[complete.cases(keep),]
# Nomes de colunas gracinhas
colnames(keep) <- c("SNP", "interaction_p")
# Extrair os significativos
keep[as.double(keep[,"interaction_p" ]) < 0.05,][,1]
```
## Problema 03.
- É honesto comparar médias usando deltas/variações? No caso, da subtração $\Delta_f-\Delta_0$ do número de pulgões não são obrigatoriamente dados brutos.
### Item a) Premissas e ANOVA
- Não faço ideia como gerar o valor da casela. Em nosso entendimento, é necessário um valor fixado de referência para que possa ser extraído o verdadeiro valor da média.
- É impossível inferir assumindo que T1 é o valor de referência exceto algebricamente. Expando.
- Se T1 fosse o valor de referência, os valores de T2.4 e T3.4 seriam $(x+2)$ e $(x+16)$, o que ainda é uma forma algébrica e não um número discreto.
- Assim, imputei um dado com base nas médias de T1 e prossegui a análise conforme indicado: $T1.4 = 9$.
#### Suposições
- Foram testadas e atendidas as seguintes condições de normalidade:
- Distribuição normal dos resíduos;
- Homoscedasticidade.
- Entretanto, as amostras não são independentes, isto é, em cada planta a resposta de uma folha para atrair ou repelir pulgões pode perfeitamente depender uma folha da outra, sendo assim as suposições clássicas da ANOVA não podem ser atendidas.
#### [WIP] Gerando os dados brutos
Antes de realizar testes de médias, foi necessário gerar valores com as médias conforme as caselas.
#### Tabela de ANOVA.
Exibo a tabela de ANOVA sem destrinchamento:
#### Código R, potencialmente errado, simular valores?
```r
# Ex 3a ----
ex3 <- read.delim(file = "./data/L02/ex3.txt")
### Testes de normalidade
library(nortest)
ex3.aov <- aov(pulg~trat, data = ex3)
summary(lm(pulg~trat, data = ex3))
lillie.test(ex3.aov$res)
shapiro.test(ex3.aov$res)
ad.test(ex3.aov$res)
leveneTest(ex3.aov$res, group = ex3.aov$model$trat)
```
### Item b) Comentar sobre as fontes de variação (FVs).
- Pelo ajuste do modelo, entendemos que a maior fonte de variação foi a aplicação de óleo ($\beta$ = 10.6) foi a maior fonte de variação em relação aos outros tratamentos, especialmente esporos, que pouco variou.
### Item c) Teste de aleatorização.
Exemplo de solução (item 3.2.1)
https://compgenomr.github.io/book/how-to-test-for-differences-between-samples.html#randomization-based-testing-for-difference-of-the-means
Google Drive
Gist
Clipboard
Sign in with Wallet
