Mercurial > repos > mora-lab > gsva

comparison GSVA.R @ 0:f94ef9b31552 draft

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:30:22 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f94ef9b31552
1 ###############################################################################
2 # title: Gene set variation analysis
3 # author: Xiaowei
4 # time: Mar.31 2021
5 ###############################################################################
6 #=================================================================
7 #how to pass parameters
8 #=================================================================
9 spec <- matrix(c("expr", 'E', 1, 'character', 'Gene expression data which is an CSV file of expression values where rows correspond to genes and columns correspond to samples.',
10 "geneSet", 'G', 1, 'character', 'Gene set',
11 'gene_Identifier_class','C', 1, 'character', 'Gene Identifier class of GeneSet',
12 'method', "M", 1,'character', 'One of gsva, ssgsea, zscore, plage',
13 'img_file', 'I', 1,'character', 'img_file',
14 'img_type', 'T', 1,'character', 'PDF, PNG, JPEG',
15 'img_width', 'W', 1, 'integer', 'the img file width',
16 'img_height', 'H', 1, 'integer', 'the img file height',
17 'GSVA_result', 'R', 1, 'character', 'Result of GSVA, an CSV file.'
18
19 ),
20 byrow = TRUE, ncol = 5)
21
22
23 if (!requireNamespace("getopt", quietly = TRUE))
24 install.packages("getopt")
25
26 opt <- getopt::getopt(spec)
27
28 #----------------
29 #整理参数
30 #----------------
31
32 if(is.null(opt$gene_Identifier_class)){gene_Identifier_class = 'Symbol'}else{gene_Identifier_class = opt$gene_Identifier_class }
33 if(is.null(opt$method)){opt$method = 'gsva'}
34 if(is.null(opt$img_type)){opt$img_type = 'PNG'}
35 if(is.null(opt$img_width)){img_width = 900}else{img_width = opt$img_width}
36 if(is.null(opt$img_height)){img_height = 900}else{img_height = opt$img_height}
37
38 #================================================================
39 #run codes
40 #================================================================
41 gsva_input_data <- read.csv(opt$expr, row.names = 1)
42
43 # if (gene_Identifier_class == 'Symbol'){
44 # geneset <- GSEABase::getGmt(opt$geneSet,
45 # geneIdType = GSEABase::SymbolIdentifier())
46 # }else{
47 # geneset <- GSEABase::getGmt(opt$geneSet,
48 # geneIdType = GSEABase::EntrezIdentifier())
49 # }
50
51 load(opt$geneSet)
52
53 result <- GSVA::gsva(as.matrix(gsva_input_data), geneSet, mx.diff=FALSE,
54 verbose=FALSE, parallel.sz=2, method = opt$method)
55 #================================================================
56 #output
57 #================================================================
58
59 write.csv(result, file = opt$GSVA_result)
60
61 if (opt$img_type == 'PNG'){
62 png(filename = opt$img_file, width = img_width, height = img_height)
63 }else if(opt$img_type == 'JPG'){
64 jpeg(filename = opt$img_file, width = img_width, height = img_height)
65 }else{
66 pdf(file = opt$img_file, width = img_width, height = img_height)
67 }
68
69 pheatmap::pheatmap(result, scale = "row", main = "heatmap", show_colnames=T)
70 dev.off()
71
72
73
74
75