10
10
# --gene_corr_pval <gene-gene corr pvalues file>
11
11
# --sig_corr <genes correlation to signature file>
12
12
13
- options(show.error.messages = FALSE ,
14
- error = function () {
15
- cat(geterrmessage(), file = stderr())
16
- q(" no" , 1 , FALSE )
17
- }
13
+ options(
14
+ show.error.messages = FALSE ,
15
+ error = function () {
16
+ cat(geterrmessage(), file = stderr())
17
+ q(" no" , 1 , FALSE )
18
+ }
18
19
)
19
20
loc <- Sys.setlocale(" LC_MESSAGES" , " en_US.UTF-8" )
20
21
@@ -23,75 +24,76 @@ library(Hmisc)
23
24
24
25
# Arguments
25
26
option_list <- list (
26
- make_option(
27
- " --sep" ,
28
- default = " \t " ,
29
- type = " character" ,
30
- help = " File separator, must be the same for all input files [default : '%default' ]"
31
- ),
32
- make_option(
33
- " --colnames" ,
34
- default = TRUE ,
35
- type = " logical" ,
36
- help = " Consider first lines as header (must stand for all input files) [default : '%default' ]"
37
- ),
38
- make_option(
39
- " --expression_file" ,
40
- default = NA ,
41
- type = " character" ,
42
- help = " Input file that contains log2(CPM +1) expression values"
43
- ),
44
- make_option(
45
- " --signatures_file" ,
46
- default = NA ,
47
- type = " character" ,
48
- help = " Input file that contains cell signature"
49
- ),
50
- make_option(
51
- " --sig_corr" ,
52
- default = " sig_corr.tsv" ,
53
- type = " character" ,
54
- help = " signature correlations output [default : '%default' ]"
55
- ),
56
- make_option(
57
- " --gene_corr" ,
58
- default = " gene_corr.tsv" ,
59
- type = " character" ,
60
- help = " genes-genes correlations output [default : '%default' ]"
61
- ),
62
- make_option(
63
- " --gene_corr_pval" ,
64
- default = " gene_corr_pval.tsv" ,
65
- type = " character" ,
66
- help = " genes-genes correlations pvalues output [default : '%default' ]"
67
- )
27
+ make_option(
28
+ " --sep" ,
29
+ default = " \t " ,
30
+ type = " character" ,
31
+ help = " File separator, must be the same for all input files [default : '%default' ]"
32
+ ),
33
+ make_option(
34
+ " --colnames" ,
35
+ default = TRUE ,
36
+ type = " logical" ,
37
+ help = " Consider first lines as header (must stand for all input files) [default : '%default' ]"
38
+ ),
39
+ make_option(
40
+ " --expression_file" ,
41
+ default = NA ,
42
+ type = " character" ,
43
+ help = " Input file that contains log2(CPM +1) expression values"
44
+ ),
45
+ make_option(
46
+ " --signatures_file" ,
47
+ default = NA ,
48
+ type = " character" ,
49
+ help = " Input file that contains cell signature"
50
+ ),
51
+ make_option(
52
+ " --sig_corr" ,
53
+ default = " sig_corr.tsv" ,
54
+ type = " character" ,
55
+ help = " signature correlations output [default : '%default' ]"
56
+ ),
57
+ make_option(
58
+ " --gene_corr" ,
59
+ default = " gene_corr.tsv" ,
60
+ type = " character" ,
61
+ help = " genes-genes correlations output [default : '%default' ]"
62
+ ),
63
+ make_option(
64
+ " --gene_corr_pval" ,
65
+ default = " gene_corr_pval.tsv" ,
66
+ type = " character" ,
67
+ help = " genes-genes correlations pvalues output [default : '%default' ]"
68
+ )
68
69
)
69
70
70
71
opt <- parse_args(OptionParser(option_list = option_list ),
71
- args = commandArgs(trailingOnly = TRUE ))
72
+ args = commandArgs(trailingOnly = TRUE )
73
+ )
72
74
73
75
if (opt $ sep == " tab" ) {
74
- opt $ sep <- " \t "
76
+ opt $ sep <- " \t "
75
77
}
76
78
if (opt $ sep == " comma" ) {
77
- opt $ sep <- " ,"
79
+ opt $ sep <- " ,"
78
80
}
79
81
80
82
# Open files
81
83
data <- read.delim(
82
- opt $ expression_file ,
83
- header = opt $ colnames ,
84
- row.names = 1 ,
85
- sep = opt $ sep ,
86
- check.names = FALSE
84
+ opt $ expression_file ,
85
+ header = opt $ colnames ,
86
+ row.names = 1 ,
87
+ sep = opt $ sep ,
88
+ check.names = FALSE
87
89
)
88
90
signature <- read.delim(
89
- opt $ signatures_file ,
90
- header = TRUE ,
91
- stringsAsFactors = FALSE ,
92
- row.names = 1 ,
93
- sep = opt $ sep ,
94
- check.names = FALSE
91
+ opt $ signatures_file ,
92
+ header = TRUE ,
93
+ stringsAsFactors = FALSE ,
94
+ row.names = 1 ,
95
+ sep = opt $ sep ,
96
+ check.names = FALSE
95
97
)
96
98
97
99
@@ -105,39 +107,41 @@ data <- rbind(t(signature), data)
105
107
gene_corr <- rcorr(t(data ), type = " pearson" ) # transpose because we correlate genes, not cells
106
108
107
109
# Gene correlation with signature score
108
- gene_signature_corr <- cbind.data.frame(gene = colnames(gene_corr $ r ),
109
- Pearson_correlation = gene_corr $ r [, 1 ],
110
- p_value = gene_corr $ P [, 1 ])
110
+ gene_signature_corr <- cbind.data.frame(
111
+ gene = colnames(gene_corr $ r ),
112
+ Pearson_correlation = gene_corr $ r [, 1 ],
113
+ p_value = gene_corr $ P [, 1 ]
114
+ )
111
115
gene_signature_corr <- gene_signature_corr [order(gene_signature_corr [, 2 ], decreasing = TRUE ), ]
112
116
113
117
114
118
# ## Save files ###
115
119
116
120
write.table(
117
- format(gene_signature_corr , digits = 2 ),
118
- file = opt $ sig_corr ,
119
- sep = " \t " ,
120
- quote = FALSE ,
121
- col.names = TRUE ,
122
- row.names = FALSE
121
+ format(gene_signature_corr , digits = 2 ),
122
+ file = opt $ sig_corr ,
123
+ sep = " \t " ,
124
+ quote = FALSE ,
125
+ col.names = TRUE ,
126
+ row.names = FALSE
123
127
)
124
128
125
129
r_genes <- data.frame (gene = rownames(gene_corr $ r ), gene_corr $ r ) # add rownames as a variable for output
126
130
write.table(
127
- format(r_genes [- 1 , - 2 ], digits = 2 ),
128
- file = opt $ gene_corr ,
129
- sep = " \t " ,
130
- quote = FALSE ,
131
- col.names = TRUE ,
132
- row.names = FALSE
131
+ format(r_genes [- 1 , - 2 ], digits = 2 ),
132
+ file = opt $ gene_corr ,
133
+ sep = " \t " ,
134
+ quote = FALSE ,
135
+ col.names = TRUE ,
136
+ row.names = FALSE
133
137
)
134
138
135
139
p_genes <- data.frame (gene = rownames(gene_corr $ P ), gene_corr $ P ) # add rownames as a variable for output
136
140
write.table(
137
- format(p_genes [- 1 , - 2 ], digits = 2 ),
138
- file = opt $ gene_corr_pval ,
139
- sep = " \t " ,
140
- quote = FALSE ,
141
- col.names = TRUE ,
142
- row.names = FALSE
141
+ format(p_genes [- 1 , - 2 ], digits = 2 ),
142
+ file = opt $ gene_corr_pval ,
143
+ sep = " \t " ,
144
+ quote = FALSE ,
145
+ col.names = TRUE ,
146
+ row.names = FALSE
143
147
)
0 commit comments