cyttools/compDiffExpr.R at master · bc2zb/cyttools · 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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#!/usr/bin/env Rscript

require(docopt)
require(methods)

"
Usage:
compDiffExpr.R (-h | --help | --version)
compDiffExpr.R DIR

Description:   This script calculates differential abundance using a feature table
Options:
--version       Show the current version.

Arguments:

DIR             Provide directory where cyttools.args.Rdata file is located
" -> doc


args <- docopt(doc)

ARGS_DIR <- args$DIR

cat("\nLoading arguments from", ARGS_DIR, "\n")

load(paste(ARGS_DIR, "cyttools.args.Rdata", sep = ""))

RESULTS_DIR <- args$OUT

source("cyttoolsFunctions.R")

panelDesign <- read.delim(args$PANEL)
targets <- read.delim(args$METADATA)

colsToCheck <- c("TimePoint", "Condition", "SampleID", "FileName", "Group")
if(checkDesignCols(targets, colsToCheck)){
  missingCols <- colsToCheck[which(colsToCheck %in% colnames(targets) == F)]
  cat("\n\nERROR: PANEL file does not include required columns.
      \n\nMissing Columns:", missingCols,
      "\n\nPlease run cyttools.R --makeMetaDataBlank to generate compatible meta data file.\n\nStopping cyttools.R\n\n")
  q()
}

targets$TimePoint <- factor(targets$TimePoint, levels = unique(targets$TimePoint))
targets$Condition <- factor(targets$Condition, levels = unique(targets$Condition))
targets$SampleID <- factor(targets$SampleID, levels = unique(targets$SampleID))

exprDesign <- targets$Condition

orderList <- gsub("\\s", ".", targets$FileName)

nodeExprTable <- read.delim(args$FEATURETABLE, row.names = 1)
nodeExprTable <- nodeExprTable[,orderVectorByListOfTerms(colnames(nodeExprTable), c(colnames(nodeExprTable)[1:2],
                                                                                    orderList))]

nodeExprTable <- nodeExprTable %>%
  mutate(RowID = paste(ConsensusCluster, Metal, sep = "_")) %>%
  select(-ConsensusCluster,
         -Metal) %>%
  column_to_rownames("RowID")

diffExprStatsTable <- tibble()

for (baseline in levels(exprDesign)){

  tmpExprDesign <- relevel(exprDesign, baseline)
  design <- model.matrix(~targets$Group + tmpExprDesign)
  colnames(design) <- gsub("tmpExprDesign", "Cnd.", colnames(design))
  colnames(design) <- gsub("targets\\$SampleID|targets\\$Group", "BatchEffect", colnames(design))

  fit <- lmFit(nodeExprTable, design = design)
  res <- eBayes(fit)

  for(experimental in (levels(tmpExprDesign)[-1])){
    topTable <- topTable(res, coef = paste0("Cnd.", experimental), number = Inf) %>%
      rownames_to_column("ClusterID") %>%
      separate(ClusterID, c("ClusterID", "Observation"), "_") %>%
      mutate(Condition = rep(experimental, nrow(.)),
             Baseline = rep(baseline, nrow(.)))
    diffExprStatsTable <- diffExprStatsTable %>% bind_rows(topTable)

  }

}

diffExprStatsTable <- diffExprStatsTable %>%
  select(ClusterID, Observation, P.Value, adj.P.Val, Condition, Baseline, logFC) %>%
  mutate(PValue = P.Value,
         FDR = adj.P.Val,
         adj.P.Val = NULL,
         P.Value = NULL)

nodeExprStatsFile <- paste(RESULTS_DIR, "nodeDifferentialExpressionTable.txt", sep = "")

write.table(diffExprStatsTable, nodeExprStatsFile, sep = "\t", quote = F, row.names = F)

# workspaceFile <- paste(RESULTS_DIR, "compDiffExprWorkspace.Rdata", sep = "")
#
# save.image(file = workspaceFile)