From d41506172b503cce93f9bf6e3e5e40d298a82eb1 Mon Sep 17 00:00:00 2001
From: Evelia Lorena Coss-Navarrete
 <49405921+EveliaCoss@users.noreply.github.com>
Date: Thu, 29 Feb 2024 15:28:19 -0600
Subject: [PATCH] Add files via upload

---
 Practica_Dia3/scripts/DEG_analysis.R  | 22 +++++++-------
 Practica_Dia3/scripts/load_data_inR.R | 42 ++++++++++++++++++++++++---
 2 files changed, 49 insertions(+), 15 deletions(-)

diff --git a/Practica_Dia3/scripts/DEG_analysis.R b/Practica_Dia3/scripts/DEG_analysis.R
index 2ed0fdb..c97eed7 100644
--- a/Practica_Dia3/scripts/DEG_analysis.R
+++ b/Practica_Dia3/scripts/DEG_analysis.R
@@ -7,8 +7,8 @@
 # Primero correr el script "load_data_inR.R"
 # Usage: Correr las lineas en un nodo de prueba en el cluster.
 # Arguments:
-#   - Input: metadata.csv, cuentas de STAR (Terminacion ReadsPerGene.out.tab)
-#   - Output: Matriz de cuentas (CSV y RData)
+#   - Input: Cargar la variable raw_counts.RData que contiene la matriz de cuentas y la metadata
+#   - Output: DEG
 #######
 
 # qlogin 
@@ -20,17 +20,15 @@ library(DESeq2)
 
 # --- Load data -----
 # Cargar archivos
-indir <- "/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/STAR_output"
 outdir <- "/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/results/"
 figdir <- '/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/results/figures/'
 
 #Cargar variable "counts", proveniente del script "load_data_inR.R"
-load("/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/results/counts/STAR_counts.RData") 
+load("/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/results/counts/raw_counts.RData") 
 samples <- metadata$sample_id # Extraer los nombres de los Transcriptomas
 metadata$type <- as.factor(metadata$type) # convertir a factor
 
 # --- DEG ----
-counts <- counts[5:129239, ] # Filtramos los rows con informacion general sobre el mapeo
 counts <- counts[which(rowSums(counts) > 10),] #Seleccionamos genes con mas de 10 cuentas
 
 # Convertir al formato dds
@@ -57,7 +55,7 @@ dds <- DESeq(dds)
 # final dispersion estimates
 # fitting model and testing
 
-# Obtener la lista de coeficientes
+# Obtener la lista de coeficientes o contrastes
 resultsNames(dds)
 
 # [1] "Intercept"                 "type_PLS_15min_vs_CONTROL"
@@ -74,17 +72,19 @@ ntd <- normTransform(dds)
 # Normalizacion de las cuentas por logaritmo y podrias hacer el analisis usando este objeto en lugar del dds
 ddslog <- rlog(dds, blind = F) 
 
+# Opcion 3. vsd
+# Estima la tendencia de dispersion de los datos y calcula la varianza, hace una normalizacion de las 
+# cuentas con respecto al tamaño de la libreria
+vsdata <- vst(dds, blind = F) 
+
+## --- Deteccion de batch effect ----
+
 # Almacenar la grafica
 png(file = paste0(figdir, "PCA_rlog.png"))
 plt <- plotPCA(ddslog, intgroup = "type")
 print(plt)
 dev.off()
 
-# Opcion 3. vsd
-# Estima la tendencia de dispersion de los datos y calcula la varianza, hace una normalizacion de las 
-# cuentas con respecto al tamaño de la libreria
-vsdata <- vst(dds, blind = F) 
-
 # Almacenar la grafica
 png(file = paste0(figdir, "PCA_vsd.png"))
 plt <- plotPCA(vsdata, intgroup = "type")
diff --git a/Practica_Dia3/scripts/load_data_inR.R b/Practica_Dia3/scripts/load_data_inR.R
index 340f5a2..a86eaf1 100644
--- a/Practica_Dia3/scripts/load_data_inR.R
+++ b/Practica_Dia3/scripts/load_data_inR.R
@@ -16,10 +16,18 @@
 
 # --- Load data -----
 # Cargar archivos
-indir <- "/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/STAR_output"
+#indir <- "/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/STAR_output"
+indir <- "/mnt/Guanina/bioinfo24/data/STAR_output/"
 outdir <- "/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/results/"
+
+# Opcion A - moverme a la carpeta y buscar
 setwd(indir)
 files <- dir(pattern = "ReadsPerGene.out.tab")
+
+# Opcion B -  sin movernos de carpeta
+files <- dir(indir, pattern = "ReadsPerGene.out.tab")
+
+# crear matriz de cuentas
 counts <- c() # esta sera la matriz
 for(i in seq_along(files)){
   x <- read.table(file = files[i], sep = "\t", header = F, as.is = T)
@@ -28,17 +36,43 @@ for(i in seq_along(files)){
 }
 
 # Cargar Metadatos
-metadata <- read.csv("/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/metadata.csv", header = F)
-# Renombrar columnas con el ID de los transcriptomas
+metadata <- read.csv("/mnt/Guanina/bioinfo24/data/metadata.csv", header = F)
+# Renombrar columnas en la metadata
 colnames(metadata) <- c("sample_id", "type")
 # Convertir a formato dataframe
 counts <- as.data.frame(counts)
 rownames(counts) <- x[,1] # Renombrar las filas con el nombre de los genes
 colnames(counts) <- sub("_ReadsPerGene.out.tab", "", files)
 
+# Eliminar las 4 primeras filas
+# counts <- counts[5:129239, ] # Filtramos los rows con informacion general sobre el mapeo
+counts <- counts[-c(1:4),]
+
 # Almacenar metadata y matriz de cuentas
 save(metadata, counts, file = paste0(outdir, "counts/raw_counts.RData"))
 write.csv(counts, file = paste0(outdir,"counts/raw_counts.csv"))
 
 # Guardar informacion de ejecucion
-sessionInfo()
\ No newline at end of file
+sessionInfo()
+
+# R version 4.0.2 (2020-06-22)
+# Platform: x86_64-pc-linux-gnu (64-bit)
+# Running under: CentOS Linux 7 (Core)
+# 
+# Matrix products: default
+# BLAS:   /cm/shared/apps/r/4.0.2-studio/lib64/R/lib/libRblas.so
+# LAPACK: /cm/shared/apps/r/4.0.2-studio/lib64/R/lib/libRlapack.so
+# 
+# locale:
+#   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
+# [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
+# [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
+# [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
+# [9] LC_ADDRESS=C               LC_TELEPHONE=C
+# [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
+# 
+# attached base packages:
+#   [1] stats     graphics  grDevices utils     datasets  methods   base
+# 
+# loaded via a namespace (and not attached):
+#   [1] compiler_4.0.2 tools_4.0.2