R_scripts_data_science/DEGs_sepsis_FDR.R at main · pinheirofabiano/R_scripts_data_science · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
rm(list = ls())


source("https://raw.githubusercontent.com/nicolau/code-R/master/DEG_analysis.R")

setwd("/Users/fabiano/Desktop/projetos_em_andamento/INOVA_USP/projeto_AUTOPSIAS/data_analysis/results")


treatedGroupArray <- c("colon.sepsis", "cortex.sepsis", "heart.sepsis", "hippocampus.sepsis", "kidney.sepsis", "lung.sepsis")
controlGroupArray <- c("colon.control", "cortex.control", "heart.control", "hippocampus.control", "kidney.control", "lung.control")

FDRcut <- 0.05

dirOut <- "DEGs_FDR"
if(!dir.exists(dirOut)) {
  dir.create(dirOut, recursive = T)
}


for(i in 1:length(treatedGroupArray)) { ################
  # i <- 1

  counts <- read.table(file = "counts.tsv", header = T, sep = "\t")
  phenodata <- read.table(file = "phenodata.tsv", header = T, sep = "\t")

  treatedGroup <- treatedGroupArray[i] ################
  controlGroup <- controlGroupArray[i] ################
  message(paste0(treatedGroup, " vs ", controlGroup))

  phenodata$Class <- paste0(phenodata$Tissue, ".", phenodata$Class) ################

  phenodata <- phenodata[ c( which( phenodata$Class == controlGroup ),
                                 which( phenodata$Class == treatedGroup ) ), ]

  phenodata <- phenodata[ phenodata$Sample %in% colnames( counts ), ]
  #head(samplesinfo)
  counts <- counts[ ,c( "Symbol", as.character( phenodata$Sample ) ) ]
  phenodata <- phenodata[ phenodata$Sample %in% colnames( counts ), ]

  rownames( counts ) <- as.character( counts$Symbol )
  counts$Symbol <- NULL

  message(paste0("Control: ", controlGroup))
  head(phenodata)

  results_DEG <- DEG_analysis(data = counts, samplesinfo = phenodata, nontreated = controlGroup,
                              treated = treatedGroup, method = "edgeR", disp = 0, verbose = T)

  write.table(x = data.frame(Symbol = rownames(results_DEG), results_DEG), file = paste0(dirOut, "/", treatedGroup, "_vs_", controlGroup, ".tsv"), quote = F, sep = "\t", row.names = F)

  DEGs <- results_DEG[results_DEG$FDR < FDRcut,]
  write.table(x = data.frame(Symbol = rownames(DEGs), DEGs),
              file = paste0(dirOut, "/", treatedGroup, "_vs_", controlGroup, "_FDR_", FDRcut, "_all_degs.tsv"),
              quote = F, sep = "\t", row.names = F)

  DEGs <- results_DEG[results_DEG$FDR < FDRcut & results_DEG$logFC > 0,]
  write.table(x = data.frame(Symbol = rownames(DEGs), DEGs),
              file = paste0(dirOut, "/", treatedGroup, "_vs_", controlGroup, "_FDR_", FDRcut, "_up-regulated.tsv"),
              quote = F, sep = "\t", row.names = F)

  DEGs <- results_DEG[results_DEG$FDR < FDRcut & results_DEG$logFC < 0,]
  write.table(x = data.frame(Symbol = rownames(DEGs), DEGs),
              file = paste0(dirOut, "/", treatedGroup, "_vs_", controlGroup, "_FDR_", FDRcut, "_down-regulated.tsv"),
              quote = F, sep = "\t", row.names = F)

  message("")
}