|
0
|
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
|
