#creador: Sofía Benítez #Grafica de Volcano
intalacion y carga de la paqueteria
#install.packages(“pacman”)
library(pacman)
p_load(“readr”, #para llamar las bases de datos “ggplot2”, #para graficar “ggrepel”, #para etiquetar datos en una gráfica “matrixTests”, # para realizar la prueba estadística “dplyr”) # facilita el manejo de datos
datos <- read_csv(file=“https://raw.githubusercontent.com/ManuelLaraMVZ/Transcript-mica/main/datos_miRNAs.csv”)
head(datos) # extacción de genes controles (referencia)
Controles <- datos %>% filter(Condicion==“Control”)
head(Controles)
promedio_controles <- Controles %>% summarise(Mean_C1 = mean(Cx1), Mean_C2 = mean(Cx2), Mean_C3 = mean(Cx3), Mean_T1 = mean(T1), Mean_T2 = mean(T2), Mean_T3 = mean(T3)) %>% mutate(Gen=“Promedio_controles”) %>% select(7,1,2,3,4,5,6) promedio_controles
#extraer los genes de la tabla “datos”
genes <- datos %>% filter(Condicion==“Target”) %>% select(-2) head(genes) ########
#Sacar el 2^-DCT
DCT <- genes %>% mutate(DCT_C1=2^-(Cx1-promedio_controles\(Mean_C1), DCT_C2=2^-(Cx2-promedio_controles\)Mean_C2), DCT_C3=2^-(Cx3-promedio_controles\(Mean_C3), DCT_T1=2^-(T1-promedio_controles\)Mean_T1), DCT_T2=2^-(T2-promedio_controles\(Mean_T2), DCT_T3=2^-(T3-promedio_controles\)Mean_T3)) %>% select(-2,-3,-4,-5,-6,-7)
DCT
promedio_genes <- DCT %>% mutate(Mean_DCT_Cx=(DCT_C1+DCT_C2+DCT_C3)/3, Mean_DCT_Tx=(DCT_T1+DCT_T2+DCT_T3)/3)
promedio_genes
#Definir prueba estadistica: t-Student
tvalue_gen <- row_t_welch(promedio_genes[,c(“DCT_C1”, “DCT_C2”, “DCT_C3”)], promedio_genes[,c(“DCT_T1”, “DCT_T2”, “DCT_T3”)])
View(tvalue_gen)
FCyPV <- promedio_genes %>% select(1, 8, 9) %>% mutate(p_value = tvalue_gen$pvalue, Fold_change = Mean_DCT_Tx/Mean_DCT_Cx) %>% select(1, 4, 5)
FCyPV
Logs <- FCyPV %>% mutate(LPV = -log10(p_value), LFC = log2(Fold_change)) %>% select(1, 4, 5)
Logs
vulcano <- ggplot(Logs, mapping = aes(x = LFC, y = LPV))+ geom_point(size = 2, color = “gray”)+ theme_classic()+ labs(title = “Analisis comparativo de miRNAs”, caption = “Creador: Sofia Benitez”, x = “Log2 (2-DDCT)”, y = “-Log10 (valor de p)”)
vulcano ######################
#limites
limite_p <- 0.05
limite_FC <- 1.5
down_regulated <- Logs %>% filter(LFC < -log2(limite_FC), LPV > -log10(limite_p))
down_regulated
up_regulated <- Logs %>% filter(LFC > log2(limite_FC), LPV > -log10(limite_p))
up_regulated
top_down_regulated <- down_regulated %>% arrange(desc(LPV)) %>% head(5)
top_down_regulated
top_up_regulated <- up_regulated %>% arrange(desc(LPV)) %>% head(5)
top_up_regulated
#mejorar la grafica
vulcano2 <- vulcano + geom_hline( yintercept = -log10(limite_p), linetype = “dashed” ) + geom_vline( xintercept = c(-log2(limite_FC), log2(limite_FC)), linetype = “dashed” )
vulcano2
vulcano3 <- vulcano2+ geom_point(data = up_regulated, x = up_regulated\(LFC, y = up_regulated\)LPV, alpha = 1, size = 3, color = “#E64B3B”)+ geom_point(data = down_regulated, x = down_regulated\(LFC, y = down_regulated\)LPV, alpha = 1, size = 3, color = “#3C5488B2”)
vulcano3
vulcano4 <- vulcano3+ geom_label_repel(data = top_up_regulated, mapping = aes(x = top_up_regulated\(LFC, y = top_up_regulated\)LPV), label = top_up_regulated\(Gen, max.overlaps = 50)+ geom_label_repel(data = top_down_regulated, mapping = aes(x = top_down_regulated\)LFC, y = top_down_regulated\(LPV), label = top_down_regulated\)Gen, max.overlaps = 50)
vulcano4
ggsave(“vulcano4.jpeg”, plot = vulcano4, height = 5, width = 6, dpi = 300)