Mercurial > repos > xiaowei > gsva

comparison GSVA.R @ 1:acd8a43b0973 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
GSVA.R --2021.4.01
author xiaowei
date Thu, 01 Apr 2021 10:19:31 +0000
parents
children
comparison
equal deleted inserted replaced
0:5b11b0530ee6 1:acd8a43b0973
1 #=================================================================
2 #how to pass parameters
3 #=================================================================
4 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.',
5 "geneSet", 'G', 1, 'character', 'Gene set',
6 'gene_Identifier_class','C', 1, 'character', 'Gene Identifier class of GeneSet',
7 'method', "M", 1,'character', 'One of gsva, ssgsea, zscore, plage',
8 'img_file', 'I', 1,'character', 'img_file',
9 'img_type', 'T', 1,'character', 'PDF, PNG, JPEG',
10 'img_width', 'W', 1, 'integer', 'the img file width',
11 'img_height', 'H', 1, 'integer', 'the img file height',
12 'GSVA_result', 'R', 1, 'character', 'Result of GSVA, an CSV file.'
13
14 ),
15 byrow = TRUE, ncol = 5)
16
17
18 if (!requireNamespace("getopt", quietly = TRUE))
19 install.packages("getopt")
20
21 opt <- getopt::getopt(spec)
22
23 #----------------
24 #整理参数
25 #----------------
26
27 if(is.null(opt$gene_Identifier_class)){gene_Identifier_class = 'Symbol'}else{gene_Identifier_class = opt$gene_Identifier_class }
28 if(is.null(opt$method)){opt$method = 'gsva'}
29 if(is.null(opt$img_type)){opt$img_type = 'PNG'}
30 if(is.null(opt$img_width)){img_width = 900}else{img_width = opt$img_width}
31 if(is.null(opt$img_height)){img_height = 900}else{img_height = opt$img_height}
32
33 #================================================================
34 #run codes
35 #================================================================
36 gsva_input_data <- read.csv(opt$expr, row.names = 1)
37
38 # if (gene_Identifier_class == 'Symbol'){
39 # geneset <- GSEABase::getGmt(opt$geneSet,
40 # geneIdType = GSEABase::SymbolIdentifier())
41 # }else{
42 # geneset <- GSEABase::getGmt(opt$geneSet,
43 # geneIdType = GSEABase::EntrezIdentifier())
44 # }
45
46 load(opt$geneSet)
47
48 result <- GSVA::gsva(as.matrix(gsva_input_data), geneSet, mx.diff=FALSE,
49 verbose=FALSE, parallel.sz=2, method = opt$method)
50 #================================================================
51 #output
52 #================================================================
53
54 write.csv(result, file = opt$GSVA_result)
55
56 if (opt$img_type == 'PNG'){
57 png(filename = opt$img_file, width = img_width, height = img_height)
58 }else if(opt$img_type == 'JPG'){
59 jpeg(filename = opt$img_file, width = img_width, height = img_height)
60 }else{
61 pdf(file = opt$img_file, width = img_width, height = img_height)
62 }
63
64 pheatmap::pheatmap(result, scale = "row", main = "heatmap", show_colnames=T)
65 dev.off()
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88 #参考：
89 #https://nbviewer.jupyter.org/github/mora-lab/benchmarks/blob/master/single-sample/workflows/Tutorial%20of%20GSVA%20using%20data%20GSE10245.ipynb
90
91 #<requirement type="package" version="1.0.12">r-pheatmap</requirement>
92
93 # library(GSVA)
94 #
95 # p <- 20000 ## number of genes
96 # n <- 30 ## number of samples
97 # nGS <- 100 ## number of gene sets
98 # min.sz <- 10 ## minimum gene set size
99 # max.sz <- 100 ## maximum gene set size
100 # X <- matrix(rnorm(p*n), nrow=p, dimnames=list(1:p, 1:n))
101 # dim(X)
102 # gs <- as.list(sample(min.sz:max.sz, size=nGS, replace=TRUE)) ## sample gene set sizes
103 # gs <- lapply(gs, function(n, p) sample(1:p, size=n, replace=FALSE), p) ## sample gene sets
104 # es.max <- gsva(X, gs, mx.diff=FALSE, verbose=FALSE, parallel.sz=1)
105 # es.dif <- gsva(X, gs, mx.diff=TRUE, verbose=FALSE, parallel.sz=1)
106
107 #boss说只需要输出行为pathway列为sample的矩阵和其heatmap就可以了
108
109 #
110 # gsva(expr, #基因表达数据，可以是SummarizedExperiment或ExpressionSet对象，也可以是行为基因列为样本的矩阵
111 # gset.idx.list, #基因集，list或GeneSetCollection 对象
112 # annotation, #基因表达数据为SummarizedExperiment对象时，用来挑选矩阵中的分子数据； 当基因表达数据为matrix，基因集为GeneSetCollection对象时，annotation参数用来注释矩阵的注释信息。 当表达数据是ExpressionSet时，该参数忽略。
113 # method=c("gsva", "ssgsea", "zscore", "plage"), #用于估计每个样本的基因集富集分数的方法
114 # kcdf=c("Gaussian", "Poisson", "none"), #表示样本间表达式水平累积分布函数非参数估计时使用的内核字符串。当基因表达值是连续时(对数)，使用"Gaussian", 当基因表达值是整数计数时(非对数)，使用"Poisson"
115 # abs.ranking=FALSE, #当mx.diff=TRUE时才会用到该参数。 abs.ranking=FALSE(默认值)，使用修正的Kuiper统计量来计算富集分数，取最大的正随机游走偏差和负随机游走偏差之间的量级差。 当abs.ranking=TRUE时，使用最大的正的和负的随机游走偏差之和的原始Kuiper统计量。
116 # min.sz=1, #gene set的最小基因数
117 # max.sz=Inf, #gene set的最大基因数
118 # parallel.sz=1L, #并行计算时要使用的执行线程数。 参数BPPARAM允许设置并行后端并优化其配置。
119 # mx.diff=TRUE, #提出从KS随机游走统计量计算富集统计量的两种方法。 mx.diff=FALSE, 计算为随机游走到0的最大距离。 TRUE（默认）时，计算最大的正的和负的随机游走偏差之间的幅度差
120 # tau=switch(method, gsva=1, ssgsea=0.25, NA), #定义由gsva和ssgsea方法执行的随机游走中尾部权重的指数。
121 # ssgsea.norm=TRUE, #当为TRUE时，运行ssgsea时，通过最小值和最大值之间的绝对差对分数进行标准化。否则忽略这步标准化
122 # verbose=TRUE, 是否打印输出每一步计算，默认是false
123 # BPPARAM=SerialParam(progressbar=verbose)) #BiocParallelParam类的一个对象，指定与此函数内的一些任务和计算的并行执行相关的参数。
124
125
126 ###########################
127 #Input
128 ###########################
129 #表达数据
130 #gene set， kegg的那些
131 #method "gsva", "ssgsea", "zscore", "plage" 之一
132
133 ######################################################
134 #boss要求输出-- 矩阵数据，热图， limma做的通路差异表达
135 ######################################################
136
137
138 #要求第一列是基因名称，第二列开始是样本
139 # gsva_input_data <- read.csv("GSVA_input.csv", row.names = 1)
140 #
141 # geneset <- GSEABase::getGmt("example_pathway.gmt",
142 # geneIdType = SymbolIdentifier())
143 #
144 # result <- gsva(as.matrix(gsva_input_data), geneset, mx.diff=FALSE,
145 # verbose=FALSE, parallel.sz=2, method = "plage")
146 #
147 # write.csv(result, file = "GSVA_result.csv")
148
149
150