annotate GSAR.R @ 1:8ff053661ae2 draft default tip

Uploaded
author mora-lab
date Thu, 20 May 2021 08:22:56 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
1 ###############################################################################
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
2 # title: Gene set analysis in R
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
3 # author: Xiaowei
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
4 # time: Mar.31 2021
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
5 ###############################################################################
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
6 #=================================================================
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
7 #how to pass parameters
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
8 #=================================================================
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
9 spec <- matrix(c("expr_file",'E', 1, 'character', 'Gene expression data which is an CSV file of expression values where rows correspond to genes and columns correspond to samples.',
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
10 "geneSet_file", 'G', 1, 'character', 'Gene set',
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
11 "design_file", 'D', 1, 'character', 'Design for the samples of expression data',
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
12 "min_size", 'I',1, 'numeric','a numeric value indicating the minimum allowed gene set size. Default value is 10.',
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
13 "max_size", 'A',1, 'numeric','a numeric value indicating the maximum allowed gene set size. Default value is 500.',
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
14 "test_method", 'T',1, 'character', "a character parameter indicating which statistical method to use for testing the gene sets. Must be one of 'GSNCAtest', 'WWtest', 'KStest', 'MDtest', 'RKStest', 'RMDtest'.",
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
15 "nperm_number", 'N', 1, 'numeric',"number of permutations used to estimate the null distribution of the test statistic. If not given, a default value 1000 is used.",
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
16 "cor_method", 'M', 1, 'character',"a character string indicating which correlation coefficient is to be computed. Possible values are 'pearson' (default), 'spearman' and 'kendall'. Default value is 'pearson'.",
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
17 "threshold_value", "V", 1, 'numeric', 'Threshold value for setting significant geneSet.',
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
18 "GSAR_output_p_value", 'R', 1, 'character',"P-value table",
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
19 "GSAR_output_plot", 'P', 1, 'character',"Plot genes relationships of significant pathway"
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
20
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
21 ),
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
22 byrow = TRUE, ncol = 5)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
23
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
24
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
25 if (!requireNamespace("getopt", quietly = TRUE))
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
26 install.packages("getopt")
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
27
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
28 opt <- getopt::getopt(spec)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
29
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
30 #----------------
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
31 #整理参数
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
32 #----------------
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
33 # expr_file
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
34 # geneSet_file
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
35 # design_file
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
36 if (is.null(opt$min_size)){min_size = 10}else{min_size = opt$min_size}
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
37 if (is.null(opt$max_size)){max_size = 500}else{max_size = opt$max_size}
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
38 if (is.null(opt$test_method)){test_method = "GSNCAtest"}else{test_method = opt$test_method}
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
39 if (is.null(opt$nperm_number)){nperm_number = 1000}else{nperm_number = opt$nperm_number}
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
40 if (is.null(opt$cor_method)){cor_method = "pearson"}else{cor_method = opt$cor_method}
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
41 if (is.null(opt$threshold_value)){threshold_value = 0.05}else{threshold_value = opt$threshold_value}
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
42
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
43
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
44 #================================================================
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
45 #run codes
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
46 #================================================================
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
47
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
48 #--- load package ------------------
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
49
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
50 suppressPackageStartupMessages(library(GSAR))
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
51 suppressPackageStartupMessages(library(GSEABase))
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
52 options(stringsAsFactors = FALSE)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
53 #---input --------------------------
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
54 # expr
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
55 data <- as.matrix(read.csv(opt$expr_file, row.names = 1))
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
56
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
57 # design
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
58 design <- read.csv(opt$design_file, row.names = 1)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
59 group <- design$group
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
60 label <- design$label
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
61
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
62 # geneSet
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
63 load(opt$geneSet_file)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
64 geneSetlist <- lapply(geneSet, geneIds)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
65 names(geneSetlist) <- names(geneSet)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
66
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
67 #-------GSAR -------------------------
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
68 # test.method <- c("GSNCAtest", "WWtest", "KStest", "MDtest", "RKStest", "RMDtest")
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
69 # “GSNCAtest”, “WWtest”, “KStest”, “MDtest”, “RKStest”, “RMDtest”
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
70 results <- TestGeneSets(object=data, group=group,
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
71 geneSets=geneSetlist, min.size=min_size,
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
72 max.size=max_size, test=test_method,
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
73 nperm = nperm_number)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
74
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
75
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
76
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
77 #================================================================
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
78 #output
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
79 #================================================================
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
80
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
81
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
82 # output p-value----------------------------------------------------
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
83 resutlts1 <- as.data.frame(t(as.data.frame(results)))
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
84 colnames(resutlts1) = "P_value"
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
85 resutlts1$geneSet <- rownames(resutlts1)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
86 resutlts1 <- resutlts1[,c("geneSet", "P_value")]
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
87 resutlts1 <- resutlts1[order(resutlts1$P_value),]
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
88 write.csv(resutlts1, file = opt$GSAR_output_p_value, row.names = FALSE)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
89 #-------------------------------------------------------------------
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
90
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
91
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
92
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
93 #plot -------------------------------------------------------------
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
94 sig.paths <- names(results[results <= threshold_value])
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
95
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
96
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
97 group1 <- unique(design[design$group == 1, "label"])
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
98 group2 <- unique(design[design$group == 2, "label"])
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
99 allgenes <- rownames(data)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
100
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
101 pdf(file = opt$GSAR_output_plot, width = 10.92)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
102 if (length(sig.paths) > 0){
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
103 for (sig.path in sig.paths){
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
104 path.index <- allgenes %in% geneSetlist[[sig.path]]
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
105 ## Plot MST2 for a pathway in two conditions
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
106 plotMST2.pathway(object=data[path.index,],
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
107 group=group, name=sig.path,
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
108 legend.size=1.2, #leg.x=-1.2, leg.y=2,
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
109 label.size=1, label.dist=0.8, cor.method=cor_method, group1.name = group1, group2.name = group2)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
110
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
111 rm(sig.path, path.index)
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
112 }
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
113
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
114 }
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
115 dev.off()
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
116
8ff053661ae2 Uploaded
mora-lab
parents:
diff changeset
117