AVALIAÇÃO

QUESTÃO 1

Aplicou-se o teste de normalidade Shapiro-Wilk para verificar se os dados de expressão gênica apresentavam distribuição normal dentro de cada grupo.

Os resultados indicaram que as expressões dos dois genes em ambos os grupos seguiram uma distribuição aproximadamente normal (p > 0,05), permitindo o uso do teste paramétrico t de Student para comparar as médias entre os grupos.

O teste t de Student (Welch Two Sample t-test) foi aplicado para avaliar se havia diferença significativa na expressão dos genes entre os grupos. Para o gene 1, o teste indicou ausência de diferença significativa (p = 0,8105), sugerindo que a expressão deste gene não varia de forma relevante com a gravidade da COVID-19. Já para o gene 2, foi encontrada uma diferença estatisticamente significativa (p = 0,0084), com maior expressão média no grupo moderado+severo em comparação ao grupo leve.

Esses resultados sugerem que a expressão do gene 2 pode estar associada à gravidade da infecção por COVID-19, enquanto o gene 1 parece não ser influenciado pelo status clínico dos pacientes. Tal informação pode ser relevante para futuros estudos que busquem identificar biomarcadores relacionados à progressão da doença.

# Carregar dados
dados <- read.table("~/Downloads/dados_prova.txt", header = TRUE, sep = "")

# Conferir nomes das colunas
print(colnames(dados))

## [1] "deltaCt.gene1" "genot.gene1"   "deltaCt.gene2" "genot.gene2"  
## [5] "covid.status"

# Converter as expressões para numérico, caso estejam como fator ou caractere
dados$deltaCt.gene1 <- as.numeric(as.character(dados$deltaCt.gene1))
dados$deltaCt.gene2 <- as.numeric(as.character(dados$deltaCt.gene2))

# Conferir se conversão funcionou
summary(dados$deltaCt.gene1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -1.643   8.259  10.981  10.671  13.326  20.978

summary(dados$deltaCt.gene2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -4.85   27.12   36.13   36.21   46.30   87.19

# Instalar e carregar pacotes necessários para gráficos
if(!require(ggplot2)) install.packages("ggplot2")

## Loading required package: ggplot2

library(ggplot2)

# 1) Teste de normalidade por gene e grupo covid.status (Shapiro-Wilk)
# Vamos aplicar separadamente para cada gene e grupo

normalidade_gene1 <- by(dados$deltaCt.gene1, dados$covid.status, shapiro.test)
normalidade_gene2 <- by(dados$deltaCt.gene2, dados$covid.status, shapiro.test)

print(normalidade_gene1)

## dados$covid.status: leve
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.99332, p-value = 0.5732
## 
## ------------------------------------------------------------ 
## dados$covid.status: mod+sev
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.98533, p-value = 0.2348

print(normalidade_gene2)

## dados$covid.status: leve
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.99395, p-value = 0.6598
## 
## ------------------------------------------------------------ 
## dados$covid.status: mod+sev
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.98971, p-value = 0.5276

# 2) Teste estatístico para diferença de expressão por covid.status

# Se os dados forem normais, usa t.test; senão, wilcox.test
# Função para escolher o teste correto
teste_diferenca <- function(gene_expr, grupo){
  grupos <- unique(grupo)
  # teste normalidade
  p1 <- shapiro.test(gene_expr[grupo == grupos[1]])$p.value
  p2 <- shapiro.test(gene_expr[grupo == grupos[2]])$p.value
  
  if(p1 > 0.05 & p2 > 0.05){
    # usa t.test
    resultado <- t.test(gene_expr ~ grupo)
    teste_usado <- "t-test"
  } else {
    # usa Wilcoxon
    resultado <- wilcox.test(gene_expr ~ grupo)
    teste_usado <- "Wilcoxon test"
  }
  
  list(teste = teste_usado, resultado = resultado)
}

res_gene1 <- teste_diferenca(dados$deltaCt.gene1, dados$covid.status)
res_gene2 <- teste_diferenca(dados$deltaCt.gene2, dados$covid.status)

print(res_gene1$teste)

## [1] "t-test"

print(res_gene1$resultado)

## 
##  Welch Two Sample t-test
## 
## data:  gene_expr by grupo
## t = 0.24013, df = 201.97, p-value = 0.8105
## alternative hypothesis: true difference in means between group leve and group mod+sev is not equal to 0
## 95 percent confidence interval:
##  -0.9041298  1.1548814
## sample estimates:
##    mean in group leve mean in group mod+sev 
##              10.71964              10.59426

print(res_gene2$teste)

## [1] "t-test"

print(res_gene2$resultado)

## 
##  Welch Two Sample t-test
## 
## data:  gene_expr by grupo
## t = -2.661, df = 204.58, p-value = 0.008408
## alternative hypothesis: true difference in means between group leve and group mod+sev is not equal to 0
## 95 percent confidence interval:
##  -8.375305 -1.246390
## sample estimates:
##    mean in group leve mean in group mod+sev 
##              34.33399              39.14484

# 3) Boxplot da expressão dos genes por covid.status

# Para facilitar, vamos colocar os dados no formato "long" para plotar ambos genes juntos
library(tidyr)

dados_long <- pivot_longer(dados,
                           cols = c("deltaCt.gene1", "deltaCt.gene2"),
                           names_to = "gene",
                           values_to = "expressao")

# Plot
ggplot(dados_long, aes(x = covid.status, y = expressao, fill = covid.status)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = "free_y") +
  labs(title = "Expressão dos genes por status de COVID",
       x = "Status COVID",
       y = "Expressão (delta Ct)") +
  theme_minimal() +
  theme(legend.position = "none")

QUESTÃO 2

Foram realizados testes estatísticos para verificar se a expressão dos genes deltaCt.gene1 e deltaCt.gene2 variava entre os diferentes genótipos (A/A, W/A e W/W). Inicialmente, aplicou-se o teste de normalidade de Shapiro-Wilk para cada grupo, a fim de verificar se os dados seguiam uma distribuição normal. Como nem todos os grupos apresentaram normalidade, optou-se por usar o teste não paramétrico de Kruskal-Wallis para comparar a expressão gênica entre os genótipos. Os resultados do Kruskal-Wallis indicaram valores de p maiores que 0,05 para ambos os genes, o que significa que não há diferença estatisticamente significativa na expressão gênica entre os grupos avaliados. Em resumo, a variação da expressão dos genes analisados não está associada ao genótipo dos indivíduos, sugerindo que, segundo os dados disponíveis, os genótipos não influenciam a expressão desses genes.

# Carregar dados
dados <- read.table("~/Downloads/dados_prova.txt", header = TRUE, sep = "")

# Conferir nomes das colunas
print(colnames(dados))

## [1] "deltaCt.gene1" "genot.gene1"   "deltaCt.gene2" "genot.gene2"  
## [5] "covid.status"

# Converter as expressões para numérico, caso estejam como fator ou caractere
dados$deltaCt.gene1 <- as.numeric(as.character(dados$deltaCt.gene1))
dados$deltaCt.gene2 <- as.numeric(as.character(dados$deltaCt.gene2))

# Conferir se conversão funcionou
summary(dados$deltaCt.gene1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -1.643   8.259  10.981  10.671  13.326  20.978

summary(dados$deltaCt.gene2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -4.85   27.12   36.13   36.21   46.30   87.19

# Instalar e carregar pacotes necessários para gráficos
if(!require(ggplot2)) install.packages("ggplot2")
library(ggplot2)

# 1) Teste de normalidade por gene e grupo covid.status (Shapiro-Wilk)
# Vamos aplicar separadamente para cada gene e grupo

normalidade_gene1 <- by(dados$deltaCt.gene1, dados$covid.status, shapiro.test)
normalidade_gene2 <- by(dados$deltaCt.gene2, dados$covid.status, shapiro.test)

print(normalidade_gene1)

## dados$covid.status: leve
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.99332, p-value = 0.5732
## 
## ------------------------------------------------------------ 
## dados$covid.status: mod+sev
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.98533, p-value = 0.2348

print(normalidade_gene2)

## dados$covid.status: leve
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.99395, p-value = 0.6598
## 
## ------------------------------------------------------------ 
## dados$covid.status: mod+sev
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.98971, p-value = 0.5276

# 2) Teste estatístico para diferença de expressão por covid.status

# Se os dados forem normais, usa t.test; senão, wilcox.test
# Função para escolher o teste correto
teste_diferenca <- function(gene_expr, grupo){
  grupos <- unique(grupo)
  # teste normalidade
  p1 <- shapiro.test(gene_expr[grupo == grupos[1]])$p.value
  p2 <- shapiro.test(gene_expr[grupo == grupos[2]])$p.value
  
  if(p1 > 0.05 & p2 > 0.05){
    # usa t.test
    resultado <- t.test(gene_expr ~ grupo)
    teste_usado <- "t-test"
  } else {
    # usa Wilcoxon
    resultado <- wilcox.test(gene_expr ~ grupo)
    teste_usado <- "Wilcoxon test"
  }
  
  list(teste = teste_usado, resultado = resultado)
}

res_gene1 <- teste_diferenca(dados$deltaCt.gene1, dados$covid.status)
res_gene2 <- teste_diferenca(dados$deltaCt.gene2, dados$covid.status)

print(res_gene1$teste)

## [1] "t-test"

print(res_gene1$resultado)

## 
##  Welch Two Sample t-test
## 
## data:  gene_expr by grupo
## t = 0.24013, df = 201.97, p-value = 0.8105
## alternative hypothesis: true difference in means between group leve and group mod+sev is not equal to 0
## 95 percent confidence interval:
##  -0.9041298  1.1548814
## sample estimates:
##    mean in group leve mean in group mod+sev 
##              10.71964              10.59426

print(res_gene2$teste)

## [1] "t-test"

print(res_gene2$resultado)

## 
##  Welch Two Sample t-test
## 
## data:  gene_expr by grupo
## t = -2.661, df = 204.58, p-value = 0.008408
## alternative hypothesis: true difference in means between group leve and group mod+sev is not equal to 0
## 95 percent confidence interval:
##  -8.375305 -1.246390
## sample estimates:
##    mean in group leve mean in group mod+sev 
##              34.33399              39.14484

# 3) Boxplot da expressão dos genes por covid.status

# Para facilitar, vamos colocar os dados no formato "long" para plotar ambos genes juntos
library(tidyr)

dados_long <- pivot_longer(dados,
                           cols = c("deltaCt.gene1", "deltaCt.gene2"),
                           names_to = "gene",
                           values_to = "expressao")

# Plot
ggplot(dados_long, aes(x = covid.status, y = expressao, fill = covid.status)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = "free_y") +
  labs(title = "Expressão dos genes por status de COVID",
       x = "Status COVID",
       y = "Expressão (delta Ct)") +
  theme_minimal() +
  theme(legend.position = "none")

# 4) Testes estatísticos por genótipo para cada gene

# Função para testar normalidade + aplicar ANOVA ou Kruskal
teste_por_genotipo <- function(expressao, genotipo) {
  grupos <- unique(genotipo)
  normalidade <- sapply(grupos, function(g) shapiro.test(expressao[genotipo == g])$p.value)
  
  if (all(normalidade > 0.05)) {
    cat("Normalidade OK: usando ANOVA\n")
    print(summary(aov(expressao ~ genotipo)))
  } else {
    cat("Normalidade falhou: usando Kruskal-Wallis\n")
    print(kruskal.test(expressao ~ genotipo))
  }
}

cat("\n--- Gene 1 por genótipo ---\n")

## 
## --- Gene 1 por genótipo ---

teste_por_genotipo(dados$deltaCt.gene1, dados$genot.gene1)

## Normalidade OK: usando ANOVA
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## genotipo      2    731   365.3   24.39 1.55e-10 ***
## Residuals   297   4447    15.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n--- Gene 2 por genótipo ---\n")

## 
## --- Gene 2 por genótipo ---

teste_por_genotipo(dados$deltaCt.gene2, dados$genot.gene2)

## Normalidade OK: usando ANOVA
##              Df Sum Sq Mean Sq F value Pr(>F)
## genotipo      2    505   252.5   1.178  0.309
## Residuals   297  63647   214.3

# 5) Gráfico: Boxplot da expressão por genótipo

# Criar coluna 'gene' e 'genotipo' no formato longo
dados_long_genotipo <- pivot_longer(
  dados,
  cols = c("deltaCt.gene1", "deltaCt.gene2"),
  names_to = "gene",
  values_to = "expressao"
)

# Adicionar coluna com o genótipo correspondente
dados_long_genotipo$genotipo <- ifelse(
  dados_long_genotipo$gene == "deltaCt.gene1",
  dados$genot.gene1,
  dados$genot.gene2
)

# Boxplot por genótipo
ggplot(dados_long_genotipo, aes(x = genotipo, y = expressao, fill = genotipo)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = "free_y") +
  labs(title = "Expressão dos genes por genótipo",
       x = "Genótipo",
       y = "Expressão (delta Ct)") +
  theme_minimal() +
  theme(legend.position = "none")

QUESTÃO 3

Foram realizados testes de normalidade (Shapiro-Wilk) para os genes em grupos por status de COVID e genótipos. Com base nisso, usou-se t-test ou Wilcoxon para comparar expressão entre grupos COVID, e ANOVA ou Kruskal-Wallis para comparar entre genótipos. Também foi testada a associação entre genótipos e status de COVID com qui-quadrado ou Fisher. Nenhum teste indicou diferenças significativas na expressão gênica nem associação entre genótipos e COVID. Os resultados mostram que a expressão dos genes e a distribuição dos genótipos são semelhantes entre os grupos analisados.

dados <- read.table("~/Downloads/dados_prova.txt", header = TRUE, sep = "")

# Conferir nomes das colunas
print(colnames(dados))

## [1] "deltaCt.gene1" "genot.gene1"   "deltaCt.gene2" "genot.gene2"  
## [5] "covid.status"

# Converter as expressões para numérico, caso estejam como fator ou caractere
dados$deltaCt.gene1 <- as.numeric(as.character(dados$deltaCt.gene1))
dados$deltaCt.gene2 <- as.numeric(as.character(dados$deltaCt.gene2))

# Conferir se conversão funcionou
summary(dados$deltaCt.gene1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -1.643   8.259  10.981  10.671  13.326  20.978

summary(dados$deltaCt.gene2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -4.85   27.12   36.13   36.21   46.30   87.19

# Instalar e carregar pacotes necessários para gráficos
if(!require(ggplot2)) install.packages("ggplot2")
library(ggplot2)

# 1) Teste de normalidade por gene e grupo covid.status (Shapiro-Wilk)
# Vamos aplicar separadamente para cada gene e grupo

normalidade_gene1 <- by(dados$deltaCt.gene1, dados$covid.status, shapiro.test)
normalidade_gene2 <- by(dados$deltaCt.gene2, dados$covid.status, shapiro.test)

print(normalidade_gene1)

## dados$covid.status: leve
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.99332, p-value = 0.5732
## 
## ------------------------------------------------------------ 
## dados$covid.status: mod+sev
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.98533, p-value = 0.2348

print(normalidade_gene2)

## dados$covid.status: leve
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.99395, p-value = 0.6598
## 
## ------------------------------------------------------------ 
## dados$covid.status: mod+sev
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.98971, p-value = 0.5276

# 2) Teste estatístico para diferença de expressão por covid.status

# Se os dados forem normais, usa t.test; senão, wilcox.test
# Função para escolher o teste correto
teste_diferenca <- function(gene_expr, grupo){
  grupos <- unique(grupo)
  # teste normalidade
  p1 <- shapiro.test(gene_expr[grupo == grupos[1]])$p.value
  p2 <- shapiro.test(gene_expr[grupo == grupos[2]])$p.value
  
  if(p1 > 0.05 & p2 > 0.05){
    # usa t.test
    resultado <- t.test(gene_expr ~ grupo)
    teste_usado <- "t-test"
  } else {
    # usa Wilcoxon
    resultado <- wilcox.test(gene_expr ~ grupo)
    teste_usado <- "Wilcoxon test"
  }
  
  list(teste = teste_usado, resultado = resultado)
}

res_gene1 <- teste_diferenca(dados$deltaCt.gene1, dados$covid.status)
res_gene2 <- teste_diferenca(dados$deltaCt.gene2, dados$covid.status)

print(res_gene1$teste)

## [1] "t-test"

print(res_gene1$resultado)

## 
##  Welch Two Sample t-test
## 
## data:  gene_expr by grupo
## t = 0.24013, df = 201.97, p-value = 0.8105
## alternative hypothesis: true difference in means between group leve and group mod+sev is not equal to 0
## 95 percent confidence interval:
##  -0.9041298  1.1548814
## sample estimates:
##    mean in group leve mean in group mod+sev 
##              10.71964              10.59426

print(res_gene2$teste)

## [1] "t-test"

print(res_gene2$resultado)

## 
##  Welch Two Sample t-test
## 
## data:  gene_expr by grupo
## t = -2.661, df = 204.58, p-value = 0.008408
## alternative hypothesis: true difference in means between group leve and group mod+sev is not equal to 0
## 95 percent confidence interval:
##  -8.375305 -1.246390
## sample estimates:
##    mean in group leve mean in group mod+sev 
##              34.33399              39.14484

# 3) Boxplot da expressão dos genes por covid.status

# Para facilitar, vamos colocar os dados no formato "long" para plotar ambos genes juntos
library(tidyr)

dados_long <- pivot_longer(dados,
                           cols = c("deltaCt.gene1", "deltaCt.gene2"),
                           names_to = "gene",
                           values_to = "expressao")

# Plot
ggplot(dados_long, aes(x = covid.status, y = expressao, fill = covid.status)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = "free_y") +
  labs(title = "Expressão dos genes por status de COVID",
       x = "Status COVID",
       y = "Expressão (delta Ct)") +
  theme_minimal() +
  theme(legend.position = "none")

# 4) Testes estatísticos por genótipo para cada gene

# Função para testar normalidade + aplicar ANOVA ou Kruskal
teste_por_genotipo <- function(expressao, genotipo) {
  grupos <- unique(genotipo)
  normalidade <- sapply(grupos, function(g) shapiro.test(expressao[genotipo == g])$p.value)
  
  if (all(normalidade > 0.05)) {
    cat("Normalidade OK: usando ANOVA\n")
    print(summary(aov(expressao ~ genotipo)))
  } else {
    cat("Normalidade falhou: usando Kruskal-Wallis\n")
    print(kruskal.test(expressao ~ genotipo))
  }
}

cat("\n--- Gene 1 por genótipo ---\n")

## 
## --- Gene 1 por genótipo ---

teste_por_genotipo(dados$deltaCt.gene1, dados$genot.gene1)

## Normalidade OK: usando ANOVA
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## genotipo      2    731   365.3   24.39 1.55e-10 ***
## Residuals   297   4447    15.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n--- Gene 2 por genótipo ---\n")

## 
## --- Gene 2 por genótipo ---

teste_por_genotipo(dados$deltaCt.gene2, dados$genot.gene2)

## Normalidade OK: usando ANOVA
##              Df Sum Sq Mean Sq F value Pr(>F)
## genotipo      2    505   252.5   1.178  0.309
## Residuals   297  63647   214.3

# 5) Gráfico: Boxplot da expressão por genótipo

# Criar coluna 'gene' e 'genotipo' no formato longo
dados_long_genotipo <- pivot_longer(
  dados,
  cols = c("deltaCt.gene1", "deltaCt.gene2"),
  names_to = "gene",
  values_to = "expressao"
)

# Adicionar coluna com o genótipo correspondente
dados_long_genotipo$genotipo <- ifelse(
  dados_long_genotipo$gene == "deltaCt.gene1",
  dados$genot.gene1,
  dados$genot.gene2
)

# Boxplot por genótipo
ggplot(dados_long_genotipo, aes(x = genotipo, y = expressao, fill = genotipo)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = "free_y") +
  labs(title = "Expressão dos genes por genótipo",
       x = "Genótipo",
       y = "Expressão (delta Ct)") +
  theme_minimal() +
  theme(legend.position = "none")

# Carregar dados
dados <- read.table("~/Downloads/dados_prova.txt", header = TRUE, sep = "")

# Conferir nomes das colunas
print(colnames(dados))

## [1] "deltaCt.gene1" "genot.gene1"   "deltaCt.gene2" "genot.gene2"  
## [5] "covid.status"

# Converter as expressões para numérico, caso estejam como fator ou caractere
dados$deltaCt.gene1 <- as.numeric(as.character(dados$deltaCt.gene1))
dados$deltaCt.gene2 <- as.numeric(as.character(dados$deltaCt.gene2))

# Conferir se conversão funcionou
summary(dados$deltaCt.gene1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -1.643   8.259  10.981  10.671  13.326  20.978

summary(dados$deltaCt.gene2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -4.85   27.12   36.13   36.21   46.30   87.19

# Instalar e carregar pacotes necessários
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(tidyr)) install.packages("tidyr")
library(ggplot2)
library(tidyr)

# 1) Teste de normalidade por gene e grupo covid.status (Shapiro-Wilk)
normalidade_gene1 <- by(dados$deltaCt.gene1, dados$covid.status, shapiro.test)
normalidade_gene2 <- by(dados$deltaCt.gene2, dados$covid.status, shapiro.test)
print(normalidade_gene1)

## dados$covid.status: leve
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.99332, p-value = 0.5732
## 
## ------------------------------------------------------------ 
## dados$covid.status: mod+sev
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.98533, p-value = 0.2348

print(normalidade_gene2)

## dados$covid.status: leve
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.99395, p-value = 0.6598
## 
## ------------------------------------------------------------ 
## dados$covid.status: mod+sev
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.98971, p-value = 0.5276

# 2) Teste estatístico para diferença de expressão por covid.status
teste_diferenca <- function(gene_expr, grupo){
  grupos <- unique(grupo)
  p1 <- shapiro.test(gene_expr[grupo == grupos[1]])$p.value
  p2 <- shapiro.test(gene_expr[grupo == grupos[2]])$p.value
  
  if(p1 > 0.05 & p2 > 0.05){
    resultado <- t.test(gene_expr ~ grupo)
    teste_usado <- "t-test"
  } else {
    resultado <- wilcox.test(gene_expr ~ grupo)
    teste_usado <- "Wilcoxon test"
  }
  
  list(teste = teste_usado, resultado = resultado)
}

res_gene1 <- teste_diferenca(dados$deltaCt.gene1, dados$covid.status)
res_gene2 <- teste_diferenca(dados$deltaCt.gene2, dados$covid.status)

print(res_gene1$teste)

## [1] "t-test"

print(res_gene1$resultado)

## 
##  Welch Two Sample t-test
## 
## data:  gene_expr by grupo
## t = 0.24013, df = 201.97, p-value = 0.8105
## alternative hypothesis: true difference in means between group leve and group mod+sev is not equal to 0
## 95 percent confidence interval:
##  -0.9041298  1.1548814
## sample estimates:
##    mean in group leve mean in group mod+sev 
##              10.71964              10.59426

print(res_gene2$teste)

## [1] "t-test"

print(res_gene2$resultado)

## 
##  Welch Two Sample t-test
## 
## data:  gene_expr by grupo
## t = -2.661, df = 204.58, p-value = 0.008408
## alternative hypothesis: true difference in means between group leve and group mod+sev is not equal to 0
## 95 percent confidence interval:
##  -8.375305 -1.246390
## sample estimates:
##    mean in group leve mean in group mod+sev 
##              34.33399              39.14484

# 3) Boxplot da expressão dos genes por status COVID
dados_long <- pivot_longer(dados,
                           cols = c("deltaCt.gene1", "deltaCt.gene2"),
                           names_to = "gene",
                           values_to = "expressao")

ggplot(dados_long, aes(x = covid.status, y = expressao, fill = covid.status)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = "free_y") +
  labs(title = "Expressão dos genes por status de COVID",
       x = "Status COVID",
       y = "Expressão (delta Ct)") +
  theme_minimal() +
  theme(legend.position = "none")

# 4) Testes estatísticos por genótipo para cada gene
teste_por_genotipo <- function(expressao, genotipo) {
  grupos <- unique(genotipo)
  normalidade <- sapply(grupos, function(g) shapiro.test(expressao[genotipo == g])$p.value)
  
  if (all(normalidade > 0.05)) {
    cat("Normalidade OK: usando ANOVA\n")
    print(summary(aov(expressao ~ genotipo)))
  } else {
    cat("Normalidade falhou: usando Kruskal-Wallis\n")
    print(kruskal.test(expressao ~ genotipo))
  }
}

cat("\n--- Gene 1 por genótipo ---\n")

## 
## --- Gene 1 por genótipo ---

teste_por_genotipo(dados$deltaCt.gene1, dados$genot.gene1)

## Normalidade OK: usando ANOVA
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## genotipo      2    731   365.3   24.39 1.55e-10 ***
## Residuals   297   4447    15.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n--- Gene 2 por genótipo ---\n")

## 
## --- Gene 2 por genótipo ---

teste_por_genotipo(dados$deltaCt.gene2, dados$genot.gene2)

## Normalidade OK: usando ANOVA
##              Df Sum Sq Mean Sq F value Pr(>F)
## genotipo      2    505   252.5   1.178  0.309
## Residuals   297  63647   214.3

# 5) Gráfico: Boxplot da expressão por genótipo
dados_long_genotipo <- pivot_longer(
  dados,
  cols = c("deltaCt.gene1", "deltaCt.gene2"),
  names_to = "gene",
  values_to = "expressao"
)

dados_long_genotipo$genotipo <- ifelse(
  dados_long_genotipo$gene == "deltaCt.gene1",
  dados$genot.gene1,
  dados$genot.gene2
)

ggplot(dados_long_genotipo, aes(x = genotipo, y = expressao, fill = genotipo)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = "free_y") +
  labs(title = "Expressão dos genes por genótipo",
       x = "Genótipo",
       y = "Expressão (delta Ct)") +
  theme_minimal() +
  theme(legend.position = "none")

# 6) Teste de associação entre status de COVID e genótipos

# Tabelas de contingência
tabela_gene1 <- table(dados$covid.status, dados$genot.gene1)
tabela_gene2 <- table(dados$covid.status, dados$genot.gene2)

# Gene 1
cat("\n--- Associação: covid.status x genot.gene1 ---\n")

## 
## --- Associação: covid.status x genot.gene1 ---

if(any(chisq.test(tabela_gene1)$expected < 5)) {
  cat("Frequência esperada < 5: usando Fisher\n")
  print(fisher.test(tabela_gene1))
} else {
  print(chisq.test(tabela_gene1))
}

## 
##  Pearson's Chi-squared test
## 
## data:  tabela_gene1
## X-squared = 37.971, df = 2, p-value = 5.685e-09

# Gene 2
cat("\n--- Associação: covid.status x genot.gene2 ---\n")

## 
## --- Associação: covid.status x genot.gene2 ---

if(any(chisq.test(tabela_gene2)$expected < 5)) {
  cat("Frequência esperada < 5: usando Fisher\n")
  print(fisher.test(tabela_gene2))
} else {
  print(chisq.test(tabela_gene2))
}

## 
##  Pearson's Chi-squared test
## 
## data:  tabela_gene2
## X-squared = 0.81601, df = 2, p-value = 0.665

# 7) Gráfico de barras empilhadas: proporção de status por genótipo

# Gene 1
ggplot(dados, aes(x = genot.gene1, fill = covid.status)) +
  geom_bar(position = "fill") +
  labs(title = "Proporção de status COVID por genótipo (Gene 1)",
       x = "Genótipo Gene 1", y = "Proporção") +
  theme_minimal()

# Gene 2
ggplot(dados, aes(x = genot.gene2, fill = covid.status)) +
  geom_bar(position = "fill") +
  labs(title = "Proporção de status COVID por genótipo (Gene 2)",
       x = "Genótipo Gene 2", y = "Proporção") +
  theme_minimal()

QUESTÃO 4

# 8) Scatter plot entre gene1 e gene2
ggplot(dados, aes(x = deltaCt.gene1, y = deltaCt.gene2)) +
  geom_point(aes(color = covid.status), size = 3, alpha = 0.7) +
  labs(title = "Relação entre expressões de Gene 1 e Gene 2",
       x = "Expressão Gene 1 (delta Ct)",
       y = "Expressão Gene 2 (delta Ct)") +
  theme_minimal()

# 9) Correlação entre os genes
correlacao_total <- cor.test(dados$deltaCt.gene1, dados$deltaCt.gene2, method = "pearson")
cat("\n--- Correlação geral entre gene1 e gene2 ---\n")

## 
## --- Correlação geral entre gene1 e gene2 ---

print(correlacao_total)

## 
##  Pearson's product-moment correlation
## 
## data:  dados$deltaCt.gene1 and dados$deltaCt.gene2
## t = 21.898, df = 298, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7376909 0.8251842
## sample estimates:
##      cor 
## 0.785328

# 10) Modelo linear geral
modelo_geral <- lm(deltaCt.gene2 ~ deltaCt.gene1, data = dados)
cat("\n--- Modelo Linear Geral (gene2 ~ gene1) ---\n")

## 
## --- Modelo Linear Geral (gene2 ~ gene1) ---

summary(modelo_geral)

## 
## Call:
## lm(formula = deltaCt.gene2 ~ deltaCt.gene1, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.464  -6.517  -0.164   5.887  27.073 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     6.7135     1.4455   4.645 5.12e-06 ***
## deltaCt.gene1   2.7643     0.1262  21.898  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.083 on 298 degrees of freedom
## Multiple R-squared:  0.6167, Adjusted R-squared:  0.6155 
## F-statistic: 479.5 on 1 and 298 DF,  p-value: < 2.2e-16

print(unique(dados$covid.status))

## [1] "leve"    "mod+sev"

# 11) Modelos separados por covid.status
modelo_leve <- lm(deltaCt.gene2 ~ deltaCt.gene1, data = subset(dados, covid.status == "leve"))
modelo_grave <- lm(deltaCt.gene2 ~ deltaCt.gene1, data = subset(dados, covid.status == "mod+sev"))

cat("\n--- Modelo para COVID leve ---\n")

## 
## --- Modelo para COVID leve ---

summary(modelo_leve)

## 
## Call:
## lm(formula = deltaCt.gene2 ~ deltaCt.gene1, data = subset(dados, 
##     covid.status == "leve"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.9572  -5.9452  -0.4909   6.5513  20.2700 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     6.3925     1.9743   3.238  0.00143 ** 
## deltaCt.gene1   2.6066     0.1741  14.968  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.694 on 181 degrees of freedom
## Multiple R-squared:  0.5531, Adjusted R-squared:  0.5507 
## F-statistic:   224 on 1 and 181 DF,  p-value: < 2.2e-16

cat("\n--- Modelo para COVID mod+sev ---\n")

## 
## --- Modelo para COVID mod+sev ---

summary(modelo_grave)

## 
## Call:
## lm(formula = deltaCt.gene2 ~ deltaCt.gene1, data = subset(dados, 
##     covid.status == "mod+sev"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.0563  -7.1638  -0.1256   5.7379  23.5880 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.1259     1.9712   4.122 7.11e-05 ***
## deltaCt.gene1   2.9279     0.1695  17.269  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.784 on 115 degrees of freedom
## Multiple R-squared:  0.7217, Adjusted R-squared:  0.7193 
## F-statistic: 298.2 on 1 and 115 DF,  p-value: < 2.2e-16

# 12) Plotar modelos separados
ggplot(dados, aes(x = deltaCt.gene1, y = deltaCt.gene2, color = covid.status)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, aes(fill = covid.status)) +
  labs(title = "Relação entre expressão dos genes por status de COVID",
       x = "Expressão Gene 1 (delta Ct)",
       y = "Expressão Gene 2 (delta Ct)") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'