rajlabmssm / echodata Goto Github PK

tidyr (NA -> 1.2.1) [CRAN]
Installing 1 packages: tidyr
Error: Failed to install 'echodata' from local:
  type 'binary' is not supported on this platform
Execution halted

Fetching echolocatoR Fine-mapping Portal study metadata.
+ 2 datasets remain after filtering.
+ Searching for multi_finemap files...
OK (HTTP 200).+ 0 unique files identified.
+ Downloading 0 files...
+ Returning local file paths.

ggtree    (NA -> 3.7.0) [CRAN]
orthogene (NA -> 1.5.0) [CRAN]

Nalls_top_SNPs
setnames(Nalls_top_SNPs,"BP","SNP_pos")

top_SNPs <- import_topSNPs(
  # topSS = "~/Desktop/Fine_Mapping/Data/GWAS/Nalls23andMe_2019/Nalls2019_TableS2.xlsx",
  topSS = Nalls_top_SNPs,
  chrom_col = "CHR", 
  position_col = "SNP_pos", 
  snp_col="SNP",
  pval_col="P, all studies", 
  effect_col="Beta, all studies", 
  gene_col="Nearest Gene", 
  locus_col = "Nearest Gene",
  grouping_vars = c("Locus"),
  remove_variants = "rs34637584")

#' Generate a named list of [e]gene-locus pairs
#'
#' @export
gene_locus_list <- function(top_SNPs){
  setNames(top_SNPs$Locus, top_SNPs$Gene)
}

make_topSNPs = function(data,
                        build = "hg19",
                        write.out = T,
                        custom_gene = NULL,
                        .metadata_file = metadata) {
  # I need to do this in advance to enter vautils function
  data = data %>%
    dplyr::rename(rsid = `#SNP`,
                  chromosome = CHR,
                  position = POS)

  snps_mapped = vautils::find_nearest_gene(as.data.frame(data),
                                           build=build,
                                           collapse=F,
                                           snp = "rsid",
                                           flanking=1000)


  if(any(snps_mapped$distance == "intergenic")) {
    snps_mapped = snps_mapped %>%
      mutate(distance = recode(distance, "intergenic" = "0"))
  }

  snps_mapped_filt = snps_mapped %>%
    mutate(distance = abs(as.numeric(distance))) %>%
    arrange(distance) %>%
    group_by(rsid) %>%
    filter(row_number() == 1) %>%
    ungroup()

  if (!is.null(custom_gene)) {
     custom_gene_filt = snps_mapped %>%
       dplyr::filter(GENE %in% custom_gene)
     # Remove the picked gene from vautils
     snps_mapped_filt = snps_mapped_filt %>%
       dplyr::filter(!((chromosome == custom_gene_filt$chromosome) & (position == custom_gene_filt$position)))

     snps_mapped_filt = rbind(snps_mapped_filt, custom_gene_filt)
  }

  snps_mapped_filt = data.frame(Locus = snps_mapped_filt$GENE,
                     Gene = snps_mapped_filt$GENE,
                     CHR = snps_mapped_filt$chromosome,
                     POS = snps_mapped_filt$position,
                     SNP = snps_mapped_filt$rsid)


  tmp_2 = data %>%
    dplyr::filter(rsid %in% snps_mapped_filt$SNP) %>%
    dplyr::select(SNP = rsid, Effect, P = Pval)
  mydf = snps_mapped_filt %>% inner_join(tmp_2)


  if (write.out) {
    fwrite(mydf, paste(.metadata_file[1],"top_SNPs.txt", sep = "/"),
           col.names = T,
           row.names = F,
           sep ="\t",
           quote = F)
  }
  return(mydf)
}

    top_SNPs <- dplyr::mutate(top_SNPs,
                              Locus=paste0("locus_chr",CHR,"_",SNP),
                              Gene=paste0("locus_chr",CHR,"_",SNP))

==> devtools::test()

ℹ Testing echodata
✔ | F W  S  OK | Context
⠏ |          0 | assign_lead_SNP                            + leadSNP missing. Assigning new one by min p-value.
✔ |          4 | assign_lead_SNP
✔ |          1 | column_dictionary                          
⠏ |          0 | construct_colmap                           Assuming fullSS_path summary stats have already been processed with MungeSumstats.
Assuming fullSS_path summary stats have already been processed with MungeSumstats.
✔ |          3 | construct_colmap
⠏ |          0 | coords_to_rsids                            Searching for RSIDs using: hg19
Loading required namespace: SNPlocs.Hsapiens.dbSNP144.GRCh37
5 / 5 SNPs annotated with RSIDs.
✔ |          5 | coords_to_rsids [6.6s]                     
✔ |          3 | createDT                                   
⠏ |          0 | dt_to_granges                              Converting dat to GRanges object.
dat is already a GRanges object.
✔ |          6 | dt_to_granges
⠏ |          0 | example_fullSS                             Writing file to ==> /var/folders/rd/rbc_wrdj4k3djf3brk6z0_dc0000gp/T//RtmpgfuhvU/nalls2019.fullSS_subset.tsv
Writing file to ==> /var/folders/rd/rbc_wrdj4k3djf3brk6z0_dc0000gp/T//RtmpgfuhvU/kunkle2019.fullSS_subset.tsv
✔ |          4 | example_fullSS                             
⠏ |          0 | extract_loci                               + Mapping colnames from MungeSumstats ==> echolocatoR
BIN1
CLU
CR1*
EPHA1
MS4A2
PICALM
PTK2B
SPI1*
   Locus SNPs
1   BIN1   55
2    CLU   17
3   CR1*   22
4  EPHA1    5
5  MS4A2  210
6 PICALM   61
7  PTK2B   17
8  SPI1*  149
✔ |          2 | extract_loci
⠏ |          0 | fillNA_CS_PP                               + Filling NAs in CS cols with 0.
+ Filling NAs in PP cols with 0.
✔ |          4 | fillNA_CS_PP
⠏ |          0 | filter_snps                                FILTER:: Filtering by SNP features.
+ FILTER:: Post-filtered data: 2913 x 26
FILTER:: Filtering by SNP features.
+ FILTER:: Post-filtered data: 6216 x 26
FILTER:: Filtering by SNP features.
+ FILTER:: Removing specified variants: rs4541502,rs34559912,rs4389574,rs35519415,rs3756246
+ FILTER:: Post-filtered data: 2908 x 26
FILTER:: Filtering by SNP features.
+ FILTER:: Post-filtered data: 2913 x 26
FILTER:: Filtering by SNP features.
BIOMART:: Trimming data to only include SNPs within gene coordinates.
BIOMART:: 3154 SNPs left after trimming.
+ leadSNP missing. Assigning new one by min p-value.
+ FILTER:: Post-filtered data: 1610 x 26
⠼ |          5 | filter_snps                                FILTER:: Filtering by SNP features.
BIOMART:: Trimming data to only include SNPs within gene coordinates.
BIOMART:: 0 SNPs left after trimming.
+ leadSNP missing. Assigning new one by min p-value.
+ FILTER:: Post-filtered data: 0 x 26
⠴ |          6 | filter_snps                                FILTER:: Filtering by SNP features.
BIOMART:: Trimming data to only include SNPs within gene coordinates.
BIOMART:: 3313 SNPs left after trimming.
+ FILTER:: Post-filtered data: 1727 x 26
✔ |          7 | filter_snps [1.3s]                         
⠏ |          0 | find_consensus_snps_no_polyfun             Identifying UCS and Consensus SNPs without PolyFun
Identifying Consensus SNPs...
+ support_thresh = 2
+ Calculating mean Posterior Probability (mean.PP)...
+ 3 fine-mapping methods used.
+ 7 Credible Set SNPs identified.
+ 1 Consensus SNPs identified.
✔ |          3 | find_consensus_snps_no_polyfun
⠏ |          0 | find_consensus_snps                        Identifying Consensus SNPs...
+ support_thresh = 2
+ Calculating mean Posterior Probability (mean.PP)...
+ 4 fine-mapping methods used.
+ 7 Credible Set SNPs identified.
+ 3 Consensus SNPs identified.
✔ |          2 | find_consensus_snps
✔ |          4 | find_top_consensus                         
⠏ |          0 | gene_locus_list                            WARNING: Dropping 97 loci not present in topSNPs: 
 - ASXL3_ASXL3
 - BAG3_BAG3
 - BIN3_BIN3
 - BRIP1_BRIP1
 - BST1_BST1
 - C5orf24_C5orf24
 - CAB39L_CAB39L
 - CAMK2D_CAMK2D
 - CASC16_CASC16
 - CD19_CD19
 - CHD9_CHD9
 - CHRNB1_CHRNB1
 - CLCN3_CLCN3
 - CNTN1_CNTN1
 - CRHR1_CRHR1
 - CRLS1_CRLS1
 - CTSB_CTSB
 - DDRGK1_DDRGK1
 - DLG2_DLG2
 - DNAH17_DNAH17
 - DYRK1A_DYRK1A
 - ELOVL7_ELOVL7
 - FAM171A2_FAM171A2
 - FAM47E_FAM47E
 - FAM47E-STBD1_FAM47E-STBD1
 - FAM49B_FAM49B
 - FBRSL1_FBRSL1
 - FCGR2A_FCGR2A
 - FGD4_FGD4
 - FGF20_FGF20
 - FYN_FYN
 - GAK_GAK
 - GALC_GALC
 - GBAP1_GBAP1
 - GBF1_GBF1
 - GCH1_GCH1
 - GPNMB_GPNMB
 - GS1-124K5.11_GS1-124K5.11
 - GXYLT1_GXYLT1
 - HIP1R_HIP1R
 - HLA-DRB5_HLA-DRB5
 - IGSF9B_IGSF9B
 - INPP5F_INPP5F
 - IP6K2_IP6K2
 - ITGA8_ITGA8
 - ITPKB_ITPKB
 - KANSL1_KANSL1
 - KCNIP3_KCNIP3
 - KCNS3_KCNS3
 - KPNA1_KPNA1
 - KRTCAP2_KRTCAP2
 - LCORL_LCORL
 - LINC00693_LINC00693
 - LOC100131289_LOC100131289
 - LRRK2_LRRK2
 - MAP3K14_MAP3K14
 - MAP4K4_MAP4K4
 - MAPT-AS1_MAPT-AS1
 - MBNL2_MBNL2
 - MCCC1_MCCC1
 - MED12L_MED12L
 - MEX3C_MEX3C
 - MIPOL1_MIPOL1
 - NOD2_NOD2
 - NSF_NSF
 - NUCKS1_NUCKS1
 - PAM_PAM
 - PMVK_PMVK
 - RAB29_RAB29
 - RETREG3_RETREG3
 - RIMS1_RIMS1
 - RIT2_RIT2
 - RNF141_RNF141
 - RPS12_RPS12
 - RPS6KL1_RPS6KL1
 - SATB1_SATB1
 - SCAF11_SCAF11
 - SCARB2_SCARB2
 - SEMA4A_SEMA4A
 - SETD1A_SETD1A
 - SH3GL2_SH3GL2
 - SIPA1L2_SIPA1L2
 - SLC44A4_SLC44A4
 - SNCA_SNCA
 - SPPL2B_SPPL2B
 - SPTSSB_SPTSSB
 - STK39_STK39
 - SYT17_SYT17
 - TMEM163_TMEM163
 - TMEM175_TMEM175
 - TRIM40_TRIM40
 - UBAP2_UBAP2
 - UBTF_UBTF
 - VAMP4_VAMP4
 - VPS13C_VPS13C
 - WNT3_WNT3
 - ZNF608_ZNF608
Reassigning gene-specific locus names.
Reassigning gene-specific locus names.
✔ |          3 | gene_locus_list
✔ |          1 | get_data                                   
✔ |          3 | get_header                                 
⠏ |          0 | get_Nalls2019_loci                         Writing fine-mapped locus data ==> /var/folders/rd/rbc_wrdj4k3djf3brk6z0_dc0000gp/T//RtmpgfuhvU/Nalls23andMe_2019/BST1.multi_finemap.csv.gz
Writing fine-mapped locus data ==> /var/folders/rd/rbc_wrdj4k3djf3brk6z0_dc0000gp/T//RtmpgfuhvU/Nalls23andMe_2019/LRRK2.multi_finemap.csv.gz
Writing fine-mapped locus data ==> /var/folders/rd/rbc_wrdj4k3djf3brk6z0_dc0000gp/T//RtmpgfuhvU/Nalls23andMe_2019/MEX3C.multi_finemap.csv.gz
✔ |          2 | get_Nalls2019_loci
⠏ |          0 | get_sample_size                            Preparing sample size column (N).
Computing effective sample size using the LDSC method:
 Neff = (N_CAS+N_CON) * (N_CAS/(N_CAS+N_CON)) / mean((N_CAS/(N_CAS+N_CON))[(N_CAS+N_CON)==max(N_CAS+N_CON)]))
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Computing effective sample size using the LDSC method:
 Neff = (N_CAS+N_CON) * (N_CAS/(N_CAS+N_CON)) / mean((N_CAS/(N_CAS+N_CON))[(N_CAS+N_CON)==max(N_CAS+N_CON)]))
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Computing effective sample size using the LDSC method:
 Neff = (N_CAS+N_CON) * (N_CAS/(N_CAS+N_CON)) / mean((N_CAS/(N_CAS+N_CON))[(N_CAS+N_CON)==max(N_CAS+N_CON)]))
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Computing effective sample size using the LDSC method:
 Neff = (N_CAS+N_CON) * (N_CAS/(N_CAS+N_CON)) / mean((N_CAS/(N_CAS+N_CON))[(N_CAS+N_CON)==max(N_CAS+N_CON)]))
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Computing effective sample size using the LDSC method:
 Neff = (N_CAS+N_CON) * (N_CAS/(N_CAS+N_CON)) / mean((N_CAS/(N_CAS+N_CON))[(N_CAS+N_CON)==max(N_CAS+N_CON)]))
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Computing effective sample size using the GIANT method:
 Neff = 2 / (1/N_CAS + 1/N_CON)
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Computing effective sample size using the METAL method:
 Neff = 4 / (1/N_CAS + 1/N_CON)
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Computing sample size using the sum method:
 N = N_CAS + N_CON
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Assigning N=10000 for all SNPs.
+ Mapping colnames from MungeSumstats ==> echolocatoR
Preparing sample size column (N).
Numeric vector supplied to compute_n.
Preparing sample size column (N).
Numeric vector supplied to compute_n.
✔ |         13 | get_sample_size
⠏ |          0 | get_SNPgroup_counts                        All loci (75) : 
           Total.SNPs          nom.sig.GWAS              sig.GWAS 
              4948.16                924.68                 82.36 
                   CS             Consensus          topConsensus 
                 7.88                  2.69                  1.47 
topConsensus.leadGWAS 
                 0.41 
Loci with at least one Consensus SNP (69) : 
           Total.SNPs          nom.sig.GWAS              sig.GWAS 
              5019.07                911.41                 84.28 
                   CS             Consensus          topConsensus 
                 7.77                  2.93                  1.59 
topConsensus.leadGWAS 
                 0.45 
✔ |          1 | get_SNPgroup_counts                        
⠏ |          0 | granges_to_bed                             Converting dat to GRanges object.
Converting 1 GRanges object to separate BED files.
Saving BED file ==> /var/folders/rd/rbc_wrdj4k3djf3brk6z0_dc0000gp/T//RtmpgfuhvU/granges1.bed.txt
✔ |          1 | granges_to_bed
⠏ |          0 | import_topSNPs                             Renaming column: BP ==> POS
Renaming column: P, all studies ==> P
Renaming column: Beta, all studies ==> Effect
Renaming column: Nearest Gene ==> Locus
Renaming column: QTL Nominated Gene (nearest QTL) ==> Gene
[1] "+ Assigning Gene and Locus independently."
Standardising column headers.
First line of summary statistics file: 
SNP	CHR	POS	Locus	Gene	Effect allele	Other allele	Effect allele frequency	Effect	SE, all studies	P	P, COJO, all studies	P, random effects, all studies	P, Conditional 23AndMe only	P, 23AndMe only	I2, all studies	Freq1, previous studies	Beta, previous studies	StdErr, previous studies	P, previous studies	I2, previous studies	Freq1, new studies	Beta, new studies	StdErr, new studies	P, new studies	I2, new studies	Passes pooled 23andMe QC	Known GWAS locus within 1MB	Failed final filtering and QC	Locus within 250KB	Locus Number	
Returning unmapped column names without making them uppercase.
+ Mapping colnames from MungeSumstats ==> echolocatoR
⠋ |          1 | import_topSNPs                             Renaming column: BP ==> POS
Renaming column: P, all studies ==> P
Renaming column: Beta, all studies ==> Effect
Renaming column: Nearest Gene ==> Locus
Renaming column: QTL Nominated Gene (nearest QTL) ==> Gene
[1] "+ Assigning Gene and Locus independently."
Standardising column headers.
First line of summary statistics file: 
SNP	CHR	POS	Locus	Gene	Effect allele	Other allele	Effect allele frequency	Effect	SE, all studies	P	P, COJO, all studies	P, random effects, all studies	P, Conditional 23AndMe only	P, 23AndMe only	I2, all studies	Freq1, previous studies	Beta, previous studies	StdErr, previous studies	P, previous studies	I2, previous studies	Freq1, new studies	Beta, new studies	StdErr, new studies	P, new studies	I2, new studies	Passes pooled 23andMe QC	Known GWAS locus within 1MB	Failed final filtering and QC	Locus within 250KB	Locus Number	
Returning unmapped column names without making them uppercase.
+ Mapping colnames from MungeSumstats ==> echolocatoR
✔ |          6 | import_topSNPs [1.6s]                      
⠏ |          0 | is_empty                                   Subset file looks good.
No SNPs identified within the summary stats file that met your criterion. :o
✔ |          6 | is_empty
✔ |          2 | is_ggbio                                   
✔ |          2 | is_ggplot                                  
✔ |          2 | is_url                                     
⠏ |          0 | limit_snps                                 Limiting to only 500 SNPs.
+ Reduced number of SNPs: 6,216 ==> 500
✔ |          2 | limit_snps
⠏ |          0 | melt_finemapping_results                   Melting PP and CS from 5 fine-mapping methods.
✔ |          4 | melt_finemapping_results
⠏ |          0 | merge_finemapping_results                  + Gathering all fine-mapping results from storage...
+ 3 multi-finemap files found.
+ Removing duplicate Multi-finemap files per locus.
+ Importing results... RtmpgfuhvU
+ Importing results... RtmpgfuhvU
+ Importing results... RtmpgfuhvU
Identifying Consensus SNPs...
+ support_thresh = 2
+ Calculating mean Posterior Probability (mean.PP)...
+ 4 fine-mapping methods used.
+ 20 Credible Set SNPs identified.
+ 9 Consensus SNPs identified.
+ Saving merged results ==> /var/folders/rd/rbc_wrdj4k3djf3brk6z0_dc0000gp/T//RtmpgfuhvU/file13a997e13e3bdmerged_results.csv.gz
✔ |          4 | merge_finemapping_results                  
✔ |          5 | merge_robust                               
⠏ |          0 | mungesumstats                              + Mapping colnames from MungeSumstats ==> echolocatoR
✔ |          5 | mungesumstats
⠏ |          0 | parquet                                    Loading required namespace: arrow
Failed with error:  'bad value'
Writing parquet file ==> /var/folders/rd/rbc_wrdj4k3djf3brk6z0_dc0000gp/T//RtmpgfuhvU/file13a994bcccdf5.parquet

 *** caught segfault ***
address 0x0, cause 'invalid permissions'

Traceback:
 1: compute__GetFunctionNames()
 2: list_compute_functions()
 3: create_binding_cache()
 4: fun(libname, pkgname)
 5: doTryCatch(return(expr), name, parentenv, handler)
 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
 7: tryCatchList(expr, classes, parentenv, handlers)
 8: tryCatch(fun(libname, pkgname), error = identity)
 9: runHook(".onLoad", env, package.lib, package)
10: loadNamespace(x)
11: echodata::write_parquet(dat)
12: eval(code, test_env)
13: eval(code, test_env)
14: withCallingHandlers({    eval(code, test_env)    if (!handled && !is.null(test)) {        skip_empty()    }}, expectation = handle_expectation, skip = handle_skip, warning = handle_warning,     message = handle_message, error = handle_error)
15: doTryCatch(return(expr), name, parentenv, handler)
16: tryCatchOne(expr, names, parentenv, handlers[[1L]])
17: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
18: doTryCatch(return(expr), name, parentenv, handler)
19: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]),     names[nh], parentenv, handlers[[nh]])
20: tryCatchList(expr, classes, parentenv, handlers)
21: tryCatch(withCallingHandlers({    eval(code, test_env)    if (!handled && !is.null(test)) {        skip_empty()    }}, expectation = handle_expectation, skip = handle_skip, warning = handle_warning,     message = handle_message, error = handle_error), error = handle_fatal,     skip = function(e) {    })
22: test_code(desc, code, env = parent.frame(), default_reporter = local_interactive_reporter())
23: test_that("read/write_parquet works", {    dat <- echodata::BST1    path <- echodata::write_parquet(dat)    dat2 <- echodata::read_parquet(path = path)})
24: eval(code, test_env)
25: eval(code, test_env)
26: withCallingHandlers({    eval(code, test_env)    if (!handled && !is.null(test)) {        skip_empty()    }}, expectation = handle_expectation, skip = handle_skip, warning = handle_warning,     message = handle_message, error = handle_error)
27: doTryCatch(return(expr), name, parentenv, handler)
28: tryCatchOne(expr, names, parentenv, handlers[[1L]])
29: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
30: doTryCatch(return(expr), name, parentenv, handler)
31: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]),     names[nh], parentenv, handlers[[nh]])
32: tryCatchList(expr, classes, parentenv, handlers)
33: tryCatch(withCallingHandlers({    eval(code, test_env)    if (!handled && !is.null(test)) {        skip_empty()    }}, expectation = handle_expectation, skip = handle_skip, warning = handle_warning,     message = handle_message, error = handle_error), error = handle_fatal,     skip = function(e) {    })
34: test_code(test = NULL, code = exprs, env = env, default_reporter = StopReporter$new())
35: source_file(path, env = env(env), desc = desc, error_call = error_call)
36: FUN(X[[i]], ...)
37: lapply(test_paths, test_one_file, env = env, desc = desc, error_call = error_call)
38: doTryCatch(return(expr), name, parentenv, handler)
39: tryCatchOne(expr, names, parentenv, handlers[[1L]])
40: tryCatchList(expr, classes, parentenv, handlers)
41: tryCatch(code, testthat_abort_reporter = function(cnd) {    cat(conditionMessage(cnd), "\n")    NULL})
42: with_reporter(reporters$multi, lapply(test_paths, test_one_file,     env = env, desc = desc, error_call = error_call))
43: test_files_serial(test_dir = test_dir, test_package = test_package,     test_paths = test_paths, load_helpers = load_helpers, reporter = reporter,     env = env, stop_on_failure = stop_on_failure, stop_on_warning = stop_on_warning,     desc = desc, load_package = load_package, error_call = error_call)
44: test_files(test_dir = path, test_paths = test_paths, test_package = package,     reporter = reporter, load_helpers = load_helpers, env = env,     stop_on_failure = stop_on_failure, stop_on_warning = stop_on_warning,     load_package = load_package, parallel = parallel)
45: test_dir(test_path, package = package, reporter = reporter, ...,     load_package = load_package)
46: testthat::test_local(pkg$path, filter = filter, stop_on_failure = stop_on_failure,     load_package = load_package, ...)
47: devtools::test()
An irrecoverable exception occurred. R is aborting now ...

Exited with status 139.