Mercurial > repos > mora-lab > mogsa

comparison mogsa.R @ 1:1f48d23544a5 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mora-lab
date Thu, 20 May 2021 08:48:44 +0000
parents
children
comparison
equal deleted inserted replaced
0:d3eba0ce0908 1:1f48d23544a5
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