Skip to content

Commit 1f0805c

Browse files
xihaolixihaoli
committed
STAARpipelineSummary v0.9.7
1 parent 748cd72 commit 1f0805c

File tree

4 files changed

+309
-4
lines changed

4 files changed

+309
-4
lines changed

‎R/Gene_Centric_Coding_Info.R

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -87,14 +87,14 @@ Gene_Centric_Coding_Info <- function(category=c("plof","plof_ds","missense","dis
8787

8888
if(category=="ptv")
8989
{
90-
results <- info_synonymous(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,gene_name=gene_name,known_loci=known_loci,rare_maf_cutoff=rare_maf_cutoff,
90+
results <- info_ptv(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,gene_name=gene_name,known_loci=known_loci,rare_maf_cutoff=rare_maf_cutoff,
9191
method_cond=method_cond,QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
9292
Annotation_dir=Annotation_dir,Annotation_name_catalog=Annotation_name_catalog,Annotation_name=Annotation_name)
9393
}
9494

9595
if(category=="ptv_ds")
9696
{
97-
results <- info_synonymous(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,gene_name=gene_name,known_loci=known_loci,rare_maf_cutoff=rare_maf_cutoff,
97+
results <- info_ptv_ds(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,gene_name=gene_name,known_loci=known_loci,rare_maf_cutoff=rare_maf_cutoff,
9898
method_cond=method_cond,QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
9999
Annotation_dir=Annotation_dir,Annotation_name_catalog=Annotation_name_catalog,Annotation_name=Annotation_name)
100100
}

‎R/Gene_Centric_Coding_Results_Summary.R

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -166,6 +166,158 @@ Gene_Centric_Coding_Results_Summary <- function(agds_dir,gene_centric_coding_job
166166
# synonymous
167167
results_synonymous_genome <- results_synonymous_genome[results_synonymous_genome[,"cMAC"]>cMAC_cutoff,]
168168

169+
### recalculate missense pvalue
170+
if(cMAC_cutoff > 0)
171+
{
172+
genes_name_disruptive_missense <- as.vector(unlist(results_disruptive_missense_genome[,1]))
173+
genes_name_missense <- as.vector(unlist(results_missense_genome[,1]))
174+
175+
recal_id <- (1:dim(results_missense_genome)[1])[(!genes_name_missense%in%genes_name_disruptive_missense)]
176+
177+
for(kk in recal_id)
178+
{
179+
print(kk)
180+
results_m <- results_missense_genome[kk,]
181+
182+
if(use_SPA)
183+
{
184+
## disruptive missense cMAC < cut_off, set p-value to 1
185+
results_missense_genome[kk,"Burden(1,25)-Disruptive"] <- 1
186+
results_missense_genome[kk,"Burden(1,1)-Disruptive"] <- 1
187+
188+
apc_num <- (length(results_m)-10)/2
189+
p_seq <- c(1:apc_num,1:apc_num+(apc_num+1))
190+
191+
## calculate STAAR-B
192+
pvalues_sub <- as.numeric(results_m[6:length(results_m)][p_seq])
193+
194+
if(sum(is.na(pvalues_sub))>0)
195+
{
196+
if(sum(is.na(pvalues_sub))==length(pvalues_sub))
197+
{
198+
results_m["STAAR-B"] <- 1
199+
}else
200+
{
201+
## not all NAs
202+
pvalues_sub <- pvalues_sub[!is.na(pvalues_sub)]
203+
if(sum(pvalues_sub[pvalues_sub<1])>0)
204+
{
205+
## not all ones
206+
results_m["STAAR-B"] <- CCT(pvalues_sub[pvalues_sub<1])
207+
208+
}else
209+
{
210+
results_m["STAAR-B"] <- 1
211+
212+
}
213+
}
214+
}else
215+
{
216+
if(sum(pvalues_sub[pvalues_sub<1])>0)
217+
{
218+
results_m["STAAR-B"] <- CCT(pvalues_sub[pvalues_sub<1])
219+
}else
220+
{
221+
results_m["STAAR-B"] <- 1
222+
}
223+
}
224+
225+
results_missense_genome[kk,"STAAR-B"] <- results_m["STAAR-B"]
226+
227+
## calculate STAAR-B(1,25)
228+
pvalues_sub <- as.numeric(results_m[6:length(results_m)][c(1:apc_num)])
229+
if(sum(is.na(pvalues_sub))>0)
230+
{
231+
if(sum(is.na(pvalues_sub))==length(pvalues_sub))
232+
{
233+
results_m["STAAR-B(1,25)"] <- 1
234+
}else
235+
{
236+
## not all NAs
237+
pvalues_sub <- pvalues_sub[!is.na(pvalues_sub)]
238+
if(sum(pvalues_sub[pvalues_sub<1])>0)
239+
{
240+
## not all ones
241+
results_m["STAAR-B(1,25)"] <- CCT(pvalues_sub[pvalues_sub<1])
242+
243+
}else
244+
{
245+
results_m["STAAR-B(1,25)"] <- 1
246+
247+
}
248+
}
249+
}else
250+
{
251+
if(sum(pvalues_sub[pvalues_sub<1])>0)
252+
{
253+
results_m["STAAR-B(1,25)"] <- CCT(pvalues_sub[pvalues_sub<1])
254+
}else
255+
{
256+
results_m["STAAR-B(1,25)"] <- 1
257+
}
258+
}
259+
260+
results_missense_genome[kk,"STAAR-B(1,25)"] <- results_m["STAAR-B(1,25)"]
261+
262+
## calculate STAAR-B(1,1)
263+
pvalues_sub <- as.numeric(results_m[6:length(results_m)][c(1:apc_num+(apc_num+1))])
264+
if(sum(is.na(pvalues_sub))>0)
265+
{
266+
if(sum(is.na(pvalues_sub))==length(pvalues_sub))
267+
{
268+
results_m["STAAR-B(1,1)"] <- 1
269+
}else
270+
{
271+
## not all NAs
272+
pvalues_sub <- pvalues_sub[!is.na(pvalues_sub)]
273+
if(sum(pvalues_sub[pvalues_sub<1])>0)
274+
{
275+
## not all ones
276+
results_m["STAAR-B(1,1)"] <- CCT(pvalues_sub[pvalues_sub<1])
277+
278+
}else
279+
{
280+
results_m["STAAR-B(1,1)"] <- 1
281+
282+
}
283+
}
284+
}else
285+
{
286+
if(sum(pvalues_sub[pvalues_sub<1])>0)
287+
{
288+
results_m["STAAR-B(1,1)"] <- CCT(pvalues_sub[pvalues_sub<1])
289+
}else
290+
{
291+
results_m["STAAR-B(1,1)"] <- 1
292+
}
293+
}
294+
295+
results_missense_genome[kk,"STAAR-B(1,1)"] <- results_m["STAAR-B(1,1)"]
296+
}else
297+
{
298+
## disruptive missense cMAC < cut_off, set p-value to 1
299+
results_missense_genome[kk,"Burden(1,25)-Disruptive"] <- 1
300+
results_missense_genome[kk,"Burden(1,1)-Disruptive"] <- 1
301+
results_missense_genome[kk,"SKAT(1,25)-Disruptive"] <- 1
302+
results_missense_genome[kk,"SKAT(1,1)-Disruptive"] <- 1
303+
results_missense_genome[kk,"ACAT-V(1,25)-Disruptive"] <- 1
304+
results_missense_genome[kk,"ACAT-V(1,1)-Disruptive"] <- 1
305+
306+
apc_num <- (length(results_m)-19)/6
307+
308+
p_seq <- c(1:apc_num,1:apc_num+(apc_num+1),1:apc_num+2*(apc_num+1),1:apc_num+3*(apc_num+1),1:apc_num+4*(apc_num+1),1:apc_num+5*(apc_num+1))
309+
310+
results_m["STAAR-O"] <- CCT(as.numeric(results_m[6:length(results_m)][p_seq]))
311+
results_m["STAAR-S(1,25)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num)]))
312+
results_m["STAAR-S(1,1)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+(apc_num+1))]))
313+
results_m["STAAR-B(1,25)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+2*(apc_num+1))]))
314+
results_m["STAAR-B(1,1)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+3*(apc_num+1))]))
315+
results_m["STAAR-A(1,25)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+4*(apc_num+1))]))
316+
results_m["STAAR-A(1,1)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+5*(apc_num+1))]))
317+
}
318+
}
319+
}
320+
169321
###### whole-genome results
170322
# plof
171323
save(results_plof_genome,file=paste0(output_path,"plof.Rdata"))

‎R/Gene_Centric_Coding_Results_Summary_incl_ptv.R

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -189,6 +189,160 @@ Gene_Centric_Coding_Results_Summary_incl_ptv <- function(agds_dir,gene_centric_c
189189
# ptv + disruptive missense
190190
results_ptv_ds_genome <- results_ptv_ds_genome[results_ptv_ds_genome[,"cMAC"]>cMAC_cutoff,]
191191

192+
193+
### recalculate missense pvalue
194+
if(cMAC_cutoff > 0)
195+
{
196+
genes_name_disruptive_missense <- as.vector(unlist(results_disruptive_missense_genome[,1]))
197+
genes_name_missense <- as.vector(unlist(results_missense_genome[,1]))
198+
199+
recal_id <- (1:dim(results_missense_genome)[1])[(!genes_name_missense%in%genes_name_disruptive_missense)]
200+
201+
for(kk in recal_id)
202+
{
203+
print(kk)
204+
results_m <- results_missense_genome[kk,]
205+
206+
if(use_SPA)
207+
{
208+
## disruptive missense cMAC < cut_off, set p-value to 1
209+
results_missense_genome[kk,"Burden(1,25)-Disruptive"] <- 1
210+
results_missense_genome[kk,"Burden(1,1)-Disruptive"] <- 1
211+
212+
apc_num <- (length(results_m)-10)/2
213+
p_seq <- c(1:apc_num,1:apc_num+(apc_num+1))
214+
215+
## calculate STAAR-B
216+
pvalues_sub <- as.numeric(results_m[6:length(results_m)][p_seq])
217+
218+
if(sum(is.na(pvalues_sub))>0)
219+
{
220+
if(sum(is.na(pvalues_sub))==length(pvalues_sub))
221+
{
222+
results_m["STAAR-B"] <- 1
223+
}else
224+
{
225+
## not all NAs
226+
pvalues_sub <- pvalues_sub[!is.na(pvalues_sub)]
227+
if(sum(pvalues_sub[pvalues_sub<1])>0)
228+
{
229+
## not all ones
230+
results_m["STAAR-B"] <- CCT(pvalues_sub[pvalues_sub<1])
231+
232+
}else
233+
{
234+
results_m["STAAR-B"] <- 1
235+
236+
}
237+
}
238+
}else
239+
{
240+
if(sum(pvalues_sub[pvalues_sub<1])>0)
241+
{
242+
results_m["STAAR-B"] <- CCT(pvalues_sub[pvalues_sub<1])
243+
}else
244+
{
245+
results_m["STAAR-B"] <- 1
246+
}
247+
}
248+
249+
results_missense_genome[kk,"STAAR-B"] <- results_m["STAAR-B"]
250+
251+
## calculate STAAR-B(1,25)
252+
pvalues_sub <- as.numeric(results_m[6:length(results_m)][c(1:apc_num)])
253+
if(sum(is.na(pvalues_sub))>0)
254+
{
255+
if(sum(is.na(pvalues_sub))==length(pvalues_sub))
256+
{
257+
results_m["STAAR-B(1,25)"] <- 1
258+
}else
259+
{
260+
## not all NAs
261+
pvalues_sub <- pvalues_sub[!is.na(pvalues_sub)]
262+
if(sum(pvalues_sub[pvalues_sub<1])>0)
263+
{
264+
## not all ones
265+
results_m["STAAR-B(1,25)"] <- CCT(pvalues_sub[pvalues_sub<1])
266+
267+
}else
268+
{
269+
results_m["STAAR-B(1,25)"] <- 1
270+
271+
}
272+
}
273+
}else
274+
{
275+
if(sum(pvalues_sub[pvalues_sub<1])>0)
276+
{
277+
results_m["STAAR-B(1,25)"] <- CCT(pvalues_sub[pvalues_sub<1])
278+
}else
279+
{
280+
results_m["STAAR-B(1,25)"] <- 1
281+
}
282+
}
283+
284+
results_missense_genome[kk,"STAAR-B(1,25)"] <- results_m["STAAR-B(1,25)"]
285+
286+
## calculate STAAR-B(1,1)
287+
pvalues_sub <- as.numeric(results_m[6:length(results_m)][c(1:apc_num+(apc_num+1))])
288+
if(sum(is.na(pvalues_sub))>0)
289+
{
290+
if(sum(is.na(pvalues_sub))==length(pvalues_sub))
291+
{
292+
results_m["STAAR-B(1,1)"] <- 1
293+
}else
294+
{
295+
## not all NAs
296+
pvalues_sub <- pvalues_sub[!is.na(pvalues_sub)]
297+
if(sum(pvalues_sub[pvalues_sub<1])>0)
298+
{
299+
## not all ones
300+
results_m["STAAR-B(1,1)"] <- CCT(pvalues_sub[pvalues_sub<1])
301+
302+
}else
303+
{
304+
results_m["STAAR-B(1,1)"] <- 1
305+
306+
}
307+
}
308+
}else
309+
{
310+
if(sum(pvalues_sub[pvalues_sub<1])>0)
311+
{
312+
results_m["STAAR-B(1,1)"] <- CCT(pvalues_sub[pvalues_sub<1])
313+
}else
314+
{
315+
results_m["STAAR-B(1,1)"] <- 1
316+
}
317+
}
318+
319+
results_missense_genome[kk,"STAAR-B(1,1)"] <- results_m["STAAR-B(1,1)"]
320+
}else
321+
{
322+
## disruptive missense cMAC < cut_off, set p-value to 1
323+
results_missense_genome[kk,"Burden(1,25)-Disruptive"] <- 1
324+
results_missense_genome[kk,"Burden(1,1)-Disruptive"] <- 1
325+
results_missense_genome[kk,"SKAT(1,25)-Disruptive"] <- 1
326+
results_missense_genome[kk,"SKAT(1,1)-Disruptive"] <- 1
327+
results_missense_genome[kk,"ACAT-V(1,25)-Disruptive"] <- 1
328+
results_missense_genome[kk,"ACAT-V(1,1)-Disruptive"] <- 1
329+
330+
apc_num <- (length(results_m)-19)/6
331+
332+
p_seq <- c(1:apc_num,1:apc_num+(apc_num+1),1:apc_num+2*(apc_num+1),1:apc_num+3*(apc_num+1),1:apc_num+4*(apc_num+1),1:apc_num+5*(apc_num+1))
333+
334+
results_m["STAAR-O"] <- CCT(as.numeric(results_m[6:length(results_m)][p_seq]))
335+
results_m["STAAR-S(1,25)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num)]))
336+
results_m["STAAR-S(1,1)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+(apc_num+1))]))
337+
results_m["STAAR-B(1,25)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+2*(apc_num+1))]))
338+
results_m["STAAR-B(1,1)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+3*(apc_num+1))]))
339+
results_m["STAAR-A(1,25)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+4*(apc_num+1))]))
340+
results_m["STAAR-A(1,1)"] <- CCT(as.numeric(results_m[6:length(results_m)][c(1:apc_num+5*(apc_num+1))]))
341+
}
342+
}
343+
}
344+
345+
192346
###### whole-genome results
193347
# plof
194348
save(results_plof_genome,file=paste0(output_path,"plof.Rdata"))

‎README.md

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
11
[![R build status](https://github.com/xihaoli/STAARpipelineSummary/workflows/R-CMD-check/badge.svg)](https://github.com/xihaoli/STAARpipelineSummary/actions)
2-
[![Build status](https://ci.appveyor.com/api/projects/status/glc569adbr7ar11j/branch/main?svg=true)](https://ci.appveyor.com/project/xihaoli/staarpipelinesummary/branch/main)
32
[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
43

54
# STAARpipelineSummary
@@ -30,7 +29,7 @@ Please see the <a href="docs/STAARpipelineSummary_manual.pdf">**STAARpipelineSum
3029
## Data Availability
3130
The whole-genome functional annotation data assembled from a variety of sources and the precomputed annotation principal components are available at the [Functional Annotation of Variant - Online Resource (FAVOR)](https://favor.genohub.org) site and [FAVOR Essential Database](https://doi.org/10.7910/DVN/1VGTJI).
3231
## Version
33-
The current version is 0.9.7 (March 23, 2024).
32+
The current version is 0.9.7 (July 25, 2024).
3433
## Citation
3534
If you use **STAARpipeline** and **STAARpipelineSummary** for your work, please cite:
3635

0 commit comments

Comments
 (0)