|
1
|
1 ###############################################################################
|
|
|
2 # title: mogsa
|
|
|
3 # author: Xiaowei
|
|
|
4 # time: Mar.31 2021
|
|
|
5 ###############################################################################
|
|
|
6 #=================================================================
|
|
|
7 #how to pass parameters
|
|
|
8 #=================================================================
|
|
|
9 spec <- matrix(c("data_file", "D",1, "character", "mogsa data file, xlsx file",
|
|
|
10 "geneSet_file", "G",1, "character", "gene set /pathway, rdata file, geneSet object",
|
|
|
11 "design_file", "S",1, "character", "CSV file, colcode, label, sample",
|
|
|
12
|
|
|
13 "PC_number", "P",1, "integer", "numbers for PCs",
|
|
|
14 "w_data","W",1, "character", "uniform, lambdal, inertia",
|
|
|
15 "proc_row", "O", 1, "character", "cnetre, center_ssq1, center_ssqN, center_ssqNm1",
|
|
|
16 "ks_B","B",1, "character", "An integer to indicate the number of bootstrapping samples to calculated the p-value of KS statistic.",
|
|
|
17 "p_adjust_method", "M", 1, "character", "p values adjust method",
|
|
|
18 "output_file1", "R", 1, "character", "genesetScoreMatrix CSV file",
|
|
|
19 "output_file2", "R2", 1, "character", "P value Matrix CSV file"
|
|
|
20
|
|
|
21 ),
|
|
|
22 byrow = TRUE, ncol = 5)
|
|
|
23
|
|
|
24
|
|
|
25 if (!requireNamespace("getopt", quietly = TRUE))
|
|
|
26 install.packages("getopt")
|
|
|
27
|
|
|
28 opt <- getopt::getopt(spec)
|
|
|
29
|
|
|
30 if (is.null(opt$PC_number)){PC_number = 3}else{PC_number = opt$PC_number}
|
|
|
31 if (is.null(opt$w_data)){w_data = "inertia"}else{w_data = opt$w_data}
|
|
|
32 if (is.null(opt$proc_row)){proc_row = "center_ssq1"}else{proc_row = opt$proc_row}
|
|
|
33 if (is.null(opt$ks_B)){ks_B = 1000}else{ks_B = opt$ks_B}
|
|
|
34 if (is.null(opt$p_adjust_method)){p_adjust_method = "fdr"}else{p_adjust_method = opt$p_adjust_method}
|
|
|
35
|
|
|
36 ##########################################################
|
|
|
37 # load packages
|
|
|
38 ##########################################################
|
|
|
39 library(mogsa)
|
|
|
40 library(openxlsx) #for load the xlsx file
|
|
|
41 ##########################################################
|
|
|
42 # Parameters
|
|
|
43 ##########################################################
|
|
|
44
|
|
|
45 sheets <- getSheetNames(opt$data_file)
|
|
|
46 mogsa_data <- vector(mode = "list", length = length(sheets))
|
|
|
47 names(mogsa_data) = sheets
|
|
|
48 for (sht in sheets){
|
|
|
49 mogsa_data[[sht]] <- read.xlsx(opt$data_file, sheet = sht, colNames = T, rowNames = T)
|
|
|
50 }
|
|
|
51
|
|
|
52
|
|
|
53
|
|
|
54 # design
|
|
|
55 design <- read.csv(opt$design_file, header = T, stringsAsFactors = F) #分类别、分组别
|
|
|
56 groups <- as.factor(design$label)
|
|
|
57 #if ("colcode" %in% colnames(design)){colcode <- design$colcode}else{colcode = NULL}
|
|
|
58
|
|
|
59
|
|
|
60 # geneSet
|
|
|
61 load(opt$geneSet_file)
|
|
|
62 sup_data <- prepSupMoa(mogsa_data, geneSets=geneSet, minMatch = 1)
|
|
|
63
|
|
|
64
|
|
|
65 ##########################################################
|
|
|
66 # mogsaBasicRun
|
|
|
67 ##########################################################
|
|
|
68 mgsa1 <- mogsa(x = mogsa_data, sup=sup_data, nf=PC_number,
|
|
|
69 proc.row = proc_row, w.data = w_data, statis = TRUE,
|
|
|
70 ks.B = ks_B, p.adjust.method = p_adjust_method)
|
|
|
71
|
|
|
72
|
|
|
73
|
|
|
74 ##########################################################
|
|
|
75 # gene Set Score (GSS) Matrix 基因集打分矩阵 每个值表示该基因集在该样本的总体活性水平
|
|
|
76 ##########################################################
|
|
|
77 # get the score matrix
|
|
|
78 scores <- getmgsa(mgsa1, "score")
|
|
|
79 write.csv(scores, file = opt$output_file1)
|
|
|
80
|
|
|
81
|
|
|
82 ##########################################################
|
|
|
83 # P value Matrix
|
|
|
84 ##########################################################
|
|
|
85 # 获取p值矩阵，行是基因集/pathway名称,列为样本名称
|
|
|
86 p.mat <- getmgsa(mgsa1, "p.val") # get p value matrix
|
|
|
87 write.csv(p.mat, file = opt$output_file2)
|
|
|
88
|
|
|
89
|
