@@ -5,43 +5,43 @@ option.list <- list(
5
5
make_option(" --pair_id" , dest = " pair_id" ),
6
6
make_option(" --abs_lib_dir" , dest = " abs_lib_dir" ),
7
7
make_option(" --pcawg" , dest = " pcawg" , default = " FALSE" ))
8
-
8
+
9
9
opt <- parse_args(OptionParser(option_list = option.list ))
10
10
print(opt )
11
-
11
+
12
12
Called_Absolute = opt [[" Called_Absoute" ]]
13
13
CGA_DIR_ABS = opt [[" abs_lib_dir" ]]
14
14
pair_id = opt [[" pair_id" ]]
15
15
PCAWG = opt [[" pcawg" ]]
16
-
16
+
17
17
Out_Name = file.path(paste(pair_id ," tsv" ,sep = " ." ))
18
-
18
+
19
19
print( paste(" sourcing files in " , CGA_DIR_ABS , sep = " " ))
20
-
20
+
21
21
rr = dir(CGA_DIR_ABS ,full.names = TRUE )
22
22
for ( i in 1 : length(rr ) ) {
23
23
source(rr [i ])}
24
24
load(Called_Absolute )
25
25
input_purity = seg.obj $ mode.res $ mode.tab [[1 ]]
26
-
26
+
27
27
allele.segs = get_hom_pairs_segtab(seg.obj )
28
-
28
+
29
29
ExtractSampleObs <<- AllelicExtractSampleObs
30
30
GetCopyRatioComb <<- AllelicGetCopyRatioComb
31
31
SCNA_model = seg.obj [[" mode.res" ]][[" mode_SCNA_models" ]][[1 ]]
32
32
tree_clust = resort_tree_clusters( SCNA_model , SCNA_model [[" seg_CCF_DP" ]][[" tree_clust" ]] )
33
-
33
+
34
34
SCNA_model [[" seg_CCF_DP" ]][[" tree_clust" ]] = tree_clust
35
35
SCNA_model [[" seg_CCF_DP" ]][[" seg_clust_tab" ]] = get_seg_clust_tab( SCNA_model )
36
36
mut_cols = SCNA_model [[" seg_CCF_DP" ]][[" tree_clust" ]][[" assign" ]]
37
-
37
+
38
38
mode.tab <- seg.obj [[" mode.res" ]][[" mode.tab" ]]
39
-
39
+
40
40
res = get_b_and_delta( mode.tab [1 , " alpha" ], mode.tab [1 , " tau" ] )
41
41
delta = res $ delta
42
42
b = res $ b
43
43
segs_d0 = deconstruct_SCNAs( SCNA_model , seg.obj , allele.segs , b , delta )
44
-
44
+
45
45
AS.seg.ix = allele.segs [, c(" seg.ix.1" , " seg.ix.2" )]
46
46
d0.allele.segs = allele.segs
47
47
d0.allele.segs [," A1.Seg.CN" ] = segs_d0 [ AS.seg.ix [,1 ] ]
@@ -52,34 +52,34 @@ d0.allele.segs["seg.ix.2"] = NULL
52
52
d0.allele.segs [" W" ] = NULL
53
53
d0.allele.segs = rename(d0.allele.segs , c(" Start.bp" = " Start" ,
54
54
" End.bp" = " End" ," n_probes" = " N_probes" ," length" = " Length" ," A1.Seg.sem" = " A1.Sigma" ," A2.Seg.sem" = " A2.Sigma" ))
55
-
55
+
56
56
write.table(d0.allele.segs , file = Out_Name , append = FALSE , quote = FALSE , sep = " \t " ,
57
57
eol = " \n " , na = " NA" , dec = " ." , row.names = FALSE ,
58
58
col.names = TRUE )
59
59
if (PCAWG == " TRUE" ){
60
60
Out_Name = file.path(paste(pair_id ," segments.txt" ,sep = " _" ))
61
61
load(Called_Absolute )
62
62
input_purity = seg.obj $ mode.res $ mode.tab [[1 ]]
63
-
63
+
64
64
ExtractSampleObs <<- AllelicExtractSampleObs
65
65
GetCopyRatioComb <<- AllelicGetCopyRatioComb
66
66
SCNA_model = seg.obj [[" mode.res" ]][[" mode_SCNA_models" ]][[1 ]]
67
67
tree_clust = resort_tree_clusters( SCNA_model , SCNA_model [[" seg_CCF_DP" ]][[" tree_clust" ]] )
68
-
68
+
69
69
SCNA_model [[" seg_CCF_DP" ]][[" tree_clust" ]] = tree_clust
70
70
SCNA_model [[" seg_CCF_DP" ]][[" seg_clust_tab" ]] = get_seg_clust_tab( SCNA_model )
71
71
mut_cols = SCNA_model [[" seg_CCF_DP" ]][[" tree_clust" ]][[" assign" ]]
72
-
72
+
73
73
mode.tab <- seg.obj [[" mode.res" ]][[" mode.tab" ]]
74
-
74
+
75
75
res = get_b_and_delta( mode.tab [1 , " alpha" ], mode.tab [1 , " tau" ] )
76
76
delta = res $ delta
77
77
b = res $ b
78
78
segs_d0 = deconstruct_SCNAs( SCNA_model , seg.obj , allele.segs , b , delta )
79
-
80
-
79
+
80
+
81
81
sc_tab = seg.obj [[" mode.res" ]][[" subclonal_SCNA_res" ]][[" subclonal_SCNA_tab" ]][1 ,,]
82
-
82
+
83
83
nc = ncol(SCNA_model [[" collapsed_CCF_dens" ]])
84
84
GRID = cumsum( c(0 , rep(1 / (nc - 1 ),(nc - 1 ))))
85
85
clust_dens = SCNA_model [[" seg_CCF_DP" ]][[" tree_clust" ]][[" CCF_dens" ]]
@@ -88,53 +88,49 @@ input_purity=seg.obj$mode.res$mode.tab[[1]]
88
88
# clust_dens_str = paste(as.character(one), collapse=", ")
89
89
ccf_hats = max.col(clust_dens )/ 100
90
90
mut_cols = SCNA_model [[" seg_CCF_DP" ]][[" tree_clust" ]][[" assign" ]]
91
-
91
+
92
92
sc_ix = mut_cols != which(clust_dens [,1 ]== 1 )
93
93
sc_ix [is.na(sc_ix ) & sc_tab [," subclonal_ix" ]== 0 ]= 0
94
94
sc_ix [sc_ix == FALSE ]= 0
95
95
sc_ix [sc_ix == TRUE ]= 1
96
-
96
+
97
97
# sc_ix_noNA=sc_ix
98
98
# sc_ix_noNA[is.na(sc_ix)]=-1
99
99
# segs_d0[sc_ix_noNA==0]=round(segs_d0[sc_ix_noNA==0])
100
-
100
+
101
101
AS.seg.ix = allele.segs [, c(" seg.ix.1" , " seg.ix.2" )]
102
102
d0.allele.segs = allele.segs
103
103
d0.allele.segs [," A1.Seg.CN" ] = round(segs_d0 [ AS.seg.ix [,1 ] ])
104
104
d0.allele.segs [," A2.Seg.CN" ] = round(segs_d0 [ AS.seg.ix [,2 ] ])
105
-
105
+
106
106
d0.allele.segs [" W" ] = NULL
107
107
d0.allele.segs = rename(d0.allele.segs , c(" Start.bp" = " Start" ,
108
108
" End.bp" = " End" ," n_probes" = " N_probes" ," length" = " Length" ," A1.Seg.sem" = " A1.Sigma" ," A2.Seg.sem" = " A2.Sigma" ))
109
-
109
+
110
110
d0.allele.segs [" seg.ix.1" ] = NULL
111
111
d0.allele.segs [" seg.ix.2" ] = NULL
112
112
d0.allele.segs [" A1.Sigma" ] = NULL
113
113
d0.allele.segs [" A2.Sigma" ] = NULL
114
-
115
-
116
-
114
+
117
115
A1.Subclonal.ix = sc_ix [AS.seg.ix [,1 ]]
118
116
A2.Subclonal.ix = sc_ix [AS.seg.ix [,2 ]]
119
117
A1.Subclonal.CN = sc_tab [," qs" ][AS.seg.ix [,1 ]]
120
118
A2.Subclonal.CN = sc_tab [," qs" ][AS.seg.ix [,2 ]]
121
119
A1.Clonal.CN = sc_tab [," qc" ][AS.seg.ix [,1 ]]
122
120
A2.Clonal.CN = sc_tab [," qc" ][AS.seg.ix [,2 ]]
123
-
124
-
125
-
121
+
126
122
A1.CCF_hat = ccf_hats [mut_cols [AS.seg.ix [,1 ]]]
127
123
A2.CCF_hat = ccf_hats [mut_cols [AS.seg.ix [,2 ]]]
128
124
A1.CCF_dens = clust_dens_str [mut_cols [AS.seg.ix [,1 ]]]
129
125
A2.CCF_dens = clust_dens_str [mut_cols [AS.seg.ix [,2 ]]]
130
-
131
-
126
+
127
+
132
128
A1.CCF_hat [is.na(A1.CCF_hat )]= sc_tab [," CCF_hat" ][AS.seg.ix [,1 ]][is.na(A1.CCF_hat )]
133
129
A2.CCF_hat [is.na(A2.CCF_hat )]= sc_tab [," CCF_hat" ][AS.seg.ix [,2 ]][is.na(A2.CCF_hat )]
134
-
130
+
135
131
d0.allele.segs [," CCF" ]= 1
136
132
d0.allele.segs [," CCF_dens" ]= clonal_clust_dens_str
137
-
133
+
138
134
clonal.ix = intersect(which(A1.Subclonal.ix == 0 ),which(A2.Subclonal.ix == 0 ))
139
135
a1.sub.ix = intersect(which(A1.Subclonal.ix != 0 ),which(A2.Subclonal.ix == 0 ))
140
136
a2.sub.ix = intersect(which(A1.Subclonal.ix == 0 ),which(A2.Subclonal.ix != 0 ))
@@ -143,49 +139,49 @@ input_purity=seg.obj$mode.res$mode.tab[[1]]
143
139
A1_Sub_Clonal_Seg = d0.allele.segs [intersect(which(A1.Subclonal.ix != 0 ),which(A2.Subclonal.ix == 0 )),]
144
140
A2_Sub_Clonal_Seg = d0.allele.segs [intersect(which(A1.Subclonal.ix == 0 ),which(A2.Subclonal.ix != 0 )),]
145
141
Both_Sub_Clonal_Seg = d0.allele.segs [intersect(which(A1.Subclonal.ix != 0 ),which(A2.Subclonal.ix != 0 )),]
146
-
147
-
142
+
143
+
148
144
d0.allele.segs [," Historically.Clonal" ]= 0
149
-
145
+
150
146
if (length(a1.sub.ix ) > 0 ){
151
147
d0.allele.segs [a1.sub.ix ," A1.Seg.CN" ]= A1.Clonal.CN [a1.sub.ix ]
152
148
d0.allele.segs [a1.sub.ix ," CCF" ]= 1 # -A1.CCF_hat[a1.sub.ix]
153
149
d0.allele.segs [a1.sub.ix ," CCF_dens" ]= clonal_clust_dens_str
154
150
d0.allele.segs [a1.sub.ix ," Historically.Clonal" ]= 1
155
-
156
-
151
+
152
+
157
153
a1_segs = d0.allele.segs [a1.sub.ix ,]
158
154
a1_segs [" A1.Seg.CN" ]= A1.Subclonal.CN [a1.sub.ix ]
159
155
a1_segs [" CCF" ]= A1.CCF_hat [a1.sub.ix ]
160
156
a1_segs [" CCF_dens" ]= A1.CCF_dens [a1.sub.ix ]
161
157
a1_segs [" Historically.Clonal" ]= 0
162
-
158
+
163
159
}else {
164
160
a1_segs = NULL }
165
161
if (length(a2.sub.ix ) > 0 ){
166
162
d0.allele.segs [a2.sub.ix ," A2.Seg.CN" ]= A2.Clonal.CN [a2.sub.ix ]
167
163
d0.allele.segs [a2.sub.ix ," CCF" ]= 1 # -A2.CCF_hat[a2.sub.ix]
168
164
d0.allele.segs [a2.sub.ix ," CCF_dens" ]= clonal_clust_dens_str
169
165
d0.allele.segs [a2.sub.ix ," Historically.Clonal" ]= 1
170
-
166
+
171
167
a2_segs = d0.allele.segs [a2.sub.ix ,]
172
168
a2_segs [" A2.Seg.CN" ]= A2.Subclonal.CN [a2.sub.ix ]
173
169
a2_segs [" CCF" ]= A2.CCF_hat [a2.sub.ix ]
174
170
a2_segs [" CCF_dens" ]= A2.CCF_dens [a2.sub.ix ]
175
171
a2_segs [" Historically.Clonal" ]= 0
176
172
}else {
177
-
173
+
178
174
a2_segs = NULL }
179
175
if (length(all.sub.ix ) > 0 ){
180
-
176
+
181
177
d0.allele.segs [all.sub.ix ," A1.Seg.CN" ]= A1.Clonal.CN [all.sub.ix ]
182
178
d0.allele.segs [all.sub.ix ," A2.Seg.CN" ]= A2.Clonal.CN [all.sub.ix ]
183
179
matched = A1.CCF_hat [all.sub.ix ] == A2.CCF_hat [all.sub.ix ]
184
-
180
+
185
181
d0.allele.segs [all.sub.ix [! matched ]," CCF" ]= 1 # -A1.CCF_hat[all.sub.ix[!matched]]-A2.CCF_hat[all.sub.ix[!matched]]
186
182
d0.allele.segs [all.sub.ix [! matched ]," CCF_dens" ]= clonal_clust_dens_str
187
183
d0.allele.segs [all.sub.ix [! matched ]," Historically.Clonal" ]= 1
188
-
184
+
189
185
d0.allele.segs [all.sub.ix [matched ]," CCF" ]= 1 # -A1.CCF_hat[all.sub.ix[matched]]
190
186
d0.allele.segs [all.sub.ix [matched ]," CCF_dens" ]= clonal_clust_dens_str
191
187
d0.allele.segs [all.sub.ix [matched ]," Historically.Clonal" ]= 1
@@ -195,17 +191,17 @@ input_purity=seg.obj$mode.res$mode.tab[[1]]
195
191
all_sub_a1_segs [" CCF" ]= A1.CCF_hat [all.sub.ix [! matched ]]
196
192
all_sub_a1_segs [" CCF_dens" ]= A1.CCF_dens [all.sub.ix [! matched ]]
197
193
all_sub_a1_segs [" Historically.Clonal" ]= 0
198
-
194
+
199
195
all_sub_a2_segs = d0.allele.segs [all.sub.ix [! matched ],]
200
196
all_sub_a2_segs [" A2.Seg.CN" ]= A2.Subclonal.CN [all.sub.ix [! matched ]]
201
197
all_sub_a2_segs [" CCF" ]= A2.CCF_hat [all.sub.ix [! matched ]]
202
198
all_sub_a2_segs [" CCF_dens" ]= A2.CCF_dens [all.sub.ix [! matched ]]
203
199
all_sub_a2_segs [" Historically.Clonal" ]= 0
204
200
}else {
205
-
201
+
206
202
all_sub_a1_segs = NULL
207
203
all_sub_a2_segs = NULL
208
-
204
+
209
205
}
210
206
if (length(all.sub.ix [matched ])> 0 ){
211
207
euqal_sub_seg = d0.allele.segs [all.sub.ix [matched ],]
@@ -214,30 +210,30 @@ input_purity=seg.obj$mode.res$mode.tab[[1]]
214
210
euqal_sub_seg [" CCF" ]= A1.CCF_hat [all.sub.ix [matched ]]
215
211
euqal_sub_seg [" CCF_dens" ]= A1.CCF_dens [all.sub.ix [matched ]]
216
212
euqal_sub_seg [" Historically.Clonal" ]= 0
217
-
213
+
218
214
}else {
219
-
215
+
220
216
euqal_sub_seg = NULL }}
221
-
217
+
222
218
else {
223
-
219
+
224
220
all_sub_a1_segs = NULL
225
221
all_sub_a2_segs = NULL
226
222
euqal_sub_seg = NULL
227
223
}
228
-
224
+
229
225
comb_segtab = rbind(d0.allele.segs ,a1_segs ,a2_segs ,all_sub_a1_segs ,all_sub_a2_segs ,euqal_sub_seg )
230
-
226
+
231
227
comb_segtab = comb_segtab [with(comb_segtab , order(as.integer(comb_segtab [[" Chromosome" ]]), as.integer(comb_segtab [[" Start" ]]))), ]
232
-
228
+
233
229
comb_segtab [," copy_number" ]= comb_segtab [[" A1.Seg.CN" ]]+ comb_segtab [[" A2.Seg.CN" ]]
234
230
comb_segtab [" N_probes" ]= NULL
235
231
comb_segtab [" Length" ]= NULL
236
-
232
+
237
233
comb_segtab = rename(comb_segtab , c(" Chromosome" = " chromosome" ," Start" = " start" ,
238
234
" End" = " end" ," A1.Seg.CN" = " minor_cn" ," A2.Seg.CN" = " major_cn" ," CCF" = " ccf" ," Historically.Clonal" = " historically_clonal" ))
239
235
comb_segtab = comb_segtab [,c(1 ,2 ,3 ,9 ,4 ,5 ,6 ,8 )]
240
-
236
+
241
237
write.table(comb_segtab , file = Out_Name , append = FALSE , quote = FALSE , sep = " \t " ,
242
238
eol = " \n " , na = " NA" , dec = " ." , row.names = FALSE ,
243
239
col.names = TRUE )
0 commit comments