From 5a37268bf07872438dc3533fcd940788b875a406 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Fri, 5 Jun 2026 10:48:51 -0400 Subject: [PATCH 01/10] data.table-ified rest of gwas_format --- R/gwas_format.R | 238 +++++++++++++++++++++++++----------------------- 1 file changed, 124 insertions(+), 114 deletions(-) diff --git a/R/gwas_format.R b/R/gwas_format.R index 138243a..0b3bab5 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -25,108 +25,117 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, sample_size, allele_freq, output_file, compute_pval = TRUE){ + # --- make data.table --- + setDT(X) + + # --- check for missing inputs --- if(missing(snp) | missing(beta_hat) | missing(se) | missing(A1) | missing(A2)){ stop("snp, beta_hat, se, A1, and A2 are required.\n") } - if(missing(chrom)){ - X <- mutate(X, chrom = NA) - chrom <- "chrom" - }else if(is.na(chrom)){ - X <- mutate(X, chrom = NA) + + if(missing(chrom) || is.na(chrom)){ + X[, chrom := NA] chrom <- "chrom" } - if(missing(pos)){ - X <- mutate(X, pos = NA_integer_) - pos <- "pos" - }else if(is.na(pos)){ - X <- mutate(X, pos = NA_integer_) + if(missing(pos) || is.na(pos)){ + X <- [, pos := NA_integer_] pos <- "pos" } - if(missing(p_value)){ - X <- mutate(X, p_value = NA_real_) - p_value <- "p_value" - p_val_missing <- TRUE - }else if(is.na(p_value)){ - X <- mutate(X, p_value = NA_real_) + if(missing(p_value) || is.na(p_value)){ + X <- [, p_value := NA_real_] p_value <- "p_value" p_val_missing <- TRUE }else{ p_val_missing <- FALSE } - if(missing(sample_size)){ - X <- mutate(X, sample_size = NA_real_) - sample_size <- "sample_size" - }else if(is.na(sample_size)){ - X <- mutate(X, sample_size = NA_real_) + if(missing(sample_size) || is.na(sample_size)){ + X <- [, sample_size := NA_real_] sample_size <- "sample_size" }else if(is.numeric(sample_size)){ - X <- mutate(X, sample_size = sample_size) + X <- [, sample_size := NA_real_] sample_size <- "sample_size" } - if(missing(allele_freq)){ - X <- mutate(X, af = NA_real_) - allele_freq <- "af" - }else if(is.na(allele_freq)){ - X <- mutate(X, af = NA_real_) + if(missing(allele_freq) || is.na(allele_freq)){ + X <- [, af := NA_real_] allele_freq <- "af" } + + # --- keep columns we want and rename --- + old_cols <- c(chrom, pos, snp, A1, A2, beta_hat, se, p_value, sample_size, allele_freq) - keep_cols <- c(chrom, pos, snp, A1, A2, beta_hat, se, p_value, sample_size, allele_freq) - - X <- X %>% - select(all_of(keep_cols))%>% - rename(snp = snp, - beta_hat =beta_hat, - se = se, - A1 = A1, - A2 = A2, - chrom = chrom, - pos = pos, - p_value = p_value, - sample_size = sample_size, - allele_freq = allele_freq) %>% - mutate(A1 = toupper(A1), - A2 = toupper(A2)) + new_cols <- c( + "chrom", "pos", "snp", "A1", "A2", + "beta_hat", "se", "p_value", "sample_size", "allele_freq" + ) - if(p_val_missing & compute_pval){ - X <- X %>% mutate(p_value = 2*pnorm(-abs(beta_hat/se))) + # checks + stopifnot(all(old_cols %chin% names(X))) + + # drop columns not needed, by reference + drop_cols <- setdiff(names(X), old_cols) + if (length(drop_cols)) { + X[, (drop_cols) := NULL] } + # reorder columns + setcolorder(X, old_cols) + + # rename columns, by reference + setnames(X, old = old_cols, new = new_cols) + + # uppercase alleles, by reference + X[, `:=`( + A1 = toupper(A1), + A2 = toupper(A2) + )] + cat("There are ", nrow(X), " variants.\n") + # --- remove invalid snps --- #Duplicated variants - dup_vars <- X$snp[which(duplicated(X$snp))] - X <- X %>% filter(!snp %in% dup_vars) + dup_vars <- X[duplicated(snp), unique(snp)] + X <- X[!snp %in% dup_vars] cat("Removing ", length(dup_vars), " duplicated variants leaving ", nrow(X), "variants.\n") #Illegal alleles - illegal_vars <- X %>% - filter((!A1 %in% c("A", "C", "T", "G") | !A2 %in% c("A", "C", "T", "G") )) %>% - select(snp) + valid_alleles <- c("A", "C", "T", "G") + illegal_vars <- X[ + !A1 %chin% valid_alleles | !A2 %chin% valid_alleles, + snp + ] if(length(illegal_vars) > 0){ - X <- X %>% filter(!snp %in% illegal_vars$snp) + X <- X[!snp %in% illegal_vars] cat("Removing ", length(illegal_vars), " variants with illegal alleles leaving ", nrow(X), "variants.\n") }else{ cat("No variants have illegal alleles.\n") } #Ambiguous alleles - n <- nrow(X) - X <- remove_ambiguous(X, upper = TRUE) - cat("Removed ", n-nrow(X), " variants with ambiguous strand.\n") + # when filtering rows, need to assign to var to keep result + X <- remove_ambiguous(X) + cat("Removed ", nrow(X), " variants with ambiguous strand.\n") - # make X into a data.table for I can do the new align_beta on it. ideally this would be for everything - setDT(X) - + # --- compute pval ---- + # ask jean do we always want to compute even if not missing + if(p_val_missing & compute_pval){ + X[, beta_hat := as.numeric(beta_hat)] + X[, se := as.numeric(se)] + + X[, p_value := fifelse( + !is.na(beta_hat) & !is.na(se) & se != 0, + 2 * pnorm(-abs(beta_hat / se)), + NA_real_ + )] + cat("Computed p-value\n") + + # --- harmonize alleles --- cat("Flipping strand and effect allele so A1 is always A\n") align_beta(X) - # data table syntax - X <- X[, .(chrom, pos, snp, A1, A2, beta_hat, se, p_value, sample_size, allele_freq)] - + # --- write out results --- if(!missing(output_file)){ cat("Writing out ", nrow(X), " variants to file.\n") # changed from path= to file= @@ -145,76 +154,77 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, # return(dat) # } -remove_ambiguous <- function(X, upper=TRUE){ - if(upper){ - X <- X %>% dplyr::filter(!(A1 == "G" & A2 == "C") & - !(A1 == "C" & A2 == "G") & - !(A1 == "A" & A2 == "T") & - !(A1 == "T" & A2 == "A")) - return(X) +remove_ambiguous <- function(X) { + stopifnot(data.table::is.data.table(X)) + stopifnot(all(c("A1", "A2") %chin% names(X))) + + ambig_pairs <- data.table(A1 = c("G", "C", "A", "T", "g", "c", "a", "t"), + A2 = c("C", "G", "T", "A", "c", "g", "t", "a")) + + idx_ambig <- X[ambig_pairs, on = .(A1, A2), which = TRUE, nomatch = NULL] + + if (!length(idx_ambig)) { + return(invisible(X)) + } else { + return(invisible(X[-idx_ambig])) } - X <- X %>% filter(!(A1 == "g" & A2 == "c") & - !(A1 == "c" & A2 == "g") & - !(A1 == "a" & A2 == "t") & - !(A1 == "t" & A2 == "a")) - return(X) } # flip signs and strands so that allele 1 is always A # now modifies X in-place w/ data table for speed and memory savings # believe beta_hat and af are always the names assigned in gwas_format. but kept for consistency w/ old func -align_beta <- function(X, upper = TRUE, - beta_col = "beta_hat", - af_col = "af") { +align_beta <- function(X, beta_col = "beta_hat", af_col = "allele_freq") { + # --- checks --- stopifnot(is.data.table(X)) - stopifnot(all(c("A1", "A2") %in% names(X))) - stopifnot(beta_col %in% names(X)) - + stopifnot(all(c("A1", "A2") %chin% names(X))) + if (!is.character(X[["A1"]])) X[, A1 := as.character(A1)] + if (!is.character(X[["A2"]])) X[, A2 := as.character(A2)] + stopifnot(beta_col %chin% names(X)) + + # --- setup --- flp <- c("A"="T","G"="C","T"="A","C"="G", "a"="t","t"="a","c"="g","g"="c") - af_present <- af_col %in% names(X) - if (!af_present) { - X[, (af_col) := NA_real_] # create af col temporarily so code can be uniform - } + af_present <- af_col %chin% names(X) - # flip strands if we have Ts to get As - X[, flip_strands_flag := - if (upper) (A1 == "T" | A2 == "T") else (A1 == "t" | A2 == "t")] - X[, `:=`( - flipped_A1 = fifelse(flip_strands_flag, flp[A1], A1), - flipped_A2 = fifelse(flip_strands_flag, flp[A2], A2) - )] + # make beta and af numeric if needed + if (!is.numeric(X[[beta_col]])) { + X[, (beta_col) := as.numeric(get(beta_col))] + } - # flag that is true if the A1 is A (we want) - X[, A1_A_flag := flipped_A1 %chin% c("A","a")] + if (af_present && !is.numeric(X[[af_col]])) { + X[, (af_col) := as.numeric(get(af_col))] + } - # swap A1 and A2 if A1 is not A - X[, `:=`( - A1 = fifelse(A1_A_flag, flipped_A1, flipped_A2), - A2 = fifelse(A1_A_flag, flipped_A2, flipped_A1) - )] + # --- flipping --- + # flip strands if we have Ts to get As + idx_flip_strands <- X[, which(A1 %chin% c("T", "t") | A2 %chin% c("T", "t"))] - # flip beta hats if the A1 is not A - # sd means subset of data (select just beta column) - X[, (beta_col) := { - b <- .SD[[1L]] - if (!is.numeric(b)) b <- as.numeric(b) - fifelse(A1_A_flag, b, -b) - }, .SDcols = beta_col] - - # flip af if the A1 is not A - X[, (af_col) := { - p <- .SD[[1L]] - if (!is.numeric(p)) p <- as.numeric(p) - fifelse(A1_A_flag, p, 1 - p) - }, .SDcols = af_col] - - # delete columns we don't need anymore - X[, c("flip_strands_flag", "flipped_A1", "flipped_A2", "A1_A_flag") := NULL] - if (!af_present) X[, (af_col) := NULL] + if (length(idx_flip_strands)) { + X[idx_flip_strands, `:=`( + A1 = unname(flp[A1]), + A2 = unname(flp[A2]) + )] + } + + # we want A1 as A + idx_swap_alleles <- X[, which(!(A1 %chin% c("A", "a")))] + + if (length(idx_swap_alleles)) { + X[idx_swap_alleles, `:=`( + A1 = A2, + A2 = A1 + )] + + # flip beta and af if A1 was not A + X[idx_swap_alleles, (beta_col) := -get(beta_col)] + if (af_present){ + X[idx_swap_alleles, (af_col) := 1 - get(af_col)] + } + } + # --- return --- # since these are in-place mods, we can call func w/o assignment to new var. but nobody wants to see the whole table return(invisible(X)) } From d24912fecfdf79a8405104abbb5592d99e101790 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Fri, 5 Jun 2026 10:55:12 -0400 Subject: [PATCH 02/10] fix ambig alleles n --- R/gwas_format.R | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/R/gwas_format.R b/R/gwas_format.R index 0b3bab5..f86086a 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -97,7 +97,7 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, # --- remove invalid snps --- #Duplicated variants dup_vars <- X[duplicated(snp), unique(snp)] - X <- X[!snp %in% dup_vars] + X <- X[!(snp %in% dup_vars)] cat("Removing ", length(dup_vars), " duplicated variants leaving ", nrow(X), "variants.\n") #Illegal alleles @@ -106,8 +106,8 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, !A1 %chin% valid_alleles | !A2 %chin% valid_alleles, snp ] - if(length(illegal_vars) > 0){ - X <- X[!snp %in% illegal_vars] + if(length(illegal_vars)){ + X <- X[!(snp %in% illegal_vars)] cat("Removing ", length(illegal_vars), " variants with illegal alleles leaving ", nrow(X), "variants.\n") }else{ cat("No variants have illegal alleles.\n") @@ -115,8 +115,9 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, #Ambiguous alleles # when filtering rows, need to assign to var to keep result + n <- nrow(X) X <- remove_ambiguous(X) - cat("Removed ", nrow(X), " variants with ambiguous strand.\n") + cat("Removed ", n-nrow(X), " variants with ambiguous strand.\n") # --- compute pval ---- # ask jean do we always want to compute even if not missing From 6f8ba2febdc5e0180bc6a3933dfa3110e82636de Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Wed, 24 Jun 2026 10:33:30 -0400 Subject: [PATCH 03/10] made gwas_format have option to return all snps --- R/.gwas_format.R.swp | Bin 0 -> 24576 bytes R/gwas_format.R | 24 ++++++++++++++++++++++-- 2 files changed, 22 insertions(+), 2 deletions(-) create mode 100644 R/.gwas_format.R.swp diff --git a/R/.gwas_format.R.swp b/R/.gwas_format.R.swp new file mode 100644 index 0000000000000000000000000000000000000000..532afdb164c9a767c5b1e7f1f8dfd5aa274dbdc3 GIT binary patch literal 24576 zcmeHP--{$i9dAuyyu_SQQBY96nY~~S%#Z9U5|5lZ?q-jb&2hQgGYcwnz1`I_eb>{~ z*Hyi9o8&NvAQ(h^(~A-k-^EuE6_Phm!8}BWM$!BOi6-}ui2C_fbmjIFsE)v8xcl5pG;QMs6{pT(-L%VJOVRX*-s+Mmdx=%3yk3@0*D{H1@cKbi(= z7|UW3=-#=-{4V4ea17ie1Er3;_uqBl6Kkucsjsx|p%1+O!c7Wu*&G9o0mp!2z%k$$ za11yG90QJl|3wCLew+7I82--6_)&F#$I|__DqN>amY&~Qg>R|u@zV2uR@XrlKVEuX z=exa%|MSxGD(G&G0mp!2z%k$$a11yG90QI4$ADwNG2j?*4BP|*LF9QW$oG3*@ZA#eg%0d5EW{%+5E3wROuI`9Qx6W9Pg1{?!k$L7GRz;A#b0Y3n~ z2YeUEfQNuj0v`e12i$lU>;ql{egRwu4uBn?4;%;X0FD8_!v@Dsfgb`-1K$U*HR62; z_#kjA@G>?zo(7%*t^?l&A|M1l4IBslhK-CDffs<^0zU(u1)c$(0L}pq16|-A;7;HM zHdS5&UIm^9o&%l*o&?SSUEnR$&zk_()31SF0ck2mK}rGcM7#8Wmkm=9YkHtd_EK4{ zhEbVwj4CJo?un?7`7lm{QQPy7=Sv-$#n8!9B=&>U4`mi5I*}P?ZFgwzfKCxaeFZrd zp0U^~rPRXjouSXrY8wewj6fP}Z&Q*f9b}?QJf65|f{bY^2AH1Z{c@k7gN+&j=kQaqL&60iAV0 z8dY7HO_|gvmxx4k2qRkNNsJUeNb_qUbZgy=i6Y_e6F>)itGB| zLWe%9p4n((jEoND614-z>%dQdH59Zm4bcM6tyu`njp|f>Vob6HWN*GKL|)e!O$ANS zG$_a-xVdo#_G}{j*~)Abi*0NLu`L~1eV5>$L2T-GFf2rn4s9L#lSHW`8!_i_RBHe| z{x~J7mjw&4?Z@VR)3DvBYJ?YcHXbZ1za~;NkRYg0BthIZ>v#-oj3+diylc#?f!1TVRf zdD4(*!z>e-tyG0}?zc<}S!@~HUf8?Q>2eQ^)>$Tnh-{zC9n!oKvm838uG3mfGfs(S zXt25vTC=5G;jsaQL03XJO<2pr@R64ogQ|f8J@v*8#t9U_W|#z*1x0}ldOFy{uv(>~ zLoZos2ph&pinbrh2?Ao6qh~@7&yHf~J~U_3^A|b@?u-t<4BJ%~RQXBAYaEv3x{S*# z_NxZp9*q^?JhW)a^` zFUt1mnj6DP6F*I887gh`3*$3XvHm z^i*1t;Vp`M(W6n12K{yT63NWE3y=97+UcL>NTyHXqiA!vO%W}P%q%#Lq%jF6Trm%- zH5pQB4I{`mqMGt;1G0J-Rm}4Cs4nwzd94+by<}9%Qq^-8Ld%$%c@qn)jr31L37!o4 zJ8j%ihD=}@r-I_3LrFy{Gi6Gu;6OaXvks$Xp2QLtoT&CJR2fg>Xb&8*Uql!~lT=w3 zt6kXEoGy&7StSY|yN)?ZMx|6J%@x=1 zH1u3YWwzd5qLWN$lyG$-6m*!=%nBxP#7vRrmYVy zV8ViYI??dxb(vvJrH{-49|~B6ExPq#*u?A+ z?P1?4lY0`|vHbGL>=7wpq+?G{@)j3wVZ$7|wOwt$Ry*2uzy6w8{|~XA{|2w=vHr*F zXmdS<@bx;rigP!|fMdWh;23ZWI0hU8jseGjW56-s7;p?Y2L3w(hac>!Vpq0s(yBhr zjbp`pAk&;{-o$Bb{Byk#A8W?h0&}37&o$$4_VQV2b1J}`{Na-)sR;I)`S@%&GS)|+ z?Abeh|Nj%L^Ixvk|2S%Hu0K|Bf581E;FkdJ349H>2Dr;H;23ZWI0hU8jseGjW56-s z7;p?Y1{?#9fw#*5p9=8!V_{WZg^zr|Kcgo6z}fx^pkD7UCr229%#K0*tMKYq(S%=v JR~rf5zW_7pphf@y literal 0 HcmV?d00001 diff --git a/R/gwas_format.R b/R/gwas_format.R index f86086a..fedc85d 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -23,7 +23,7 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, chrom, pos, p_value, sample_size, allele_freq, - output_file, compute_pval = TRUE){ + output_file, compute_pval = TRUE, return_og_snps = FALSE){ # --- make data.table --- setDT(X) @@ -95,7 +95,13 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, cat("There are ", nrow(X), " variants.\n") # --- remove invalid snps --- + # before filtering begins, label rows with indices, so we can join back to get all snps if they want + # use ids instead of snp names bc may not be unique + X[, row_id := .I] + og_snp_index <- X[, .(row_id, snp)] + #Duplicated variants + # drop ALL instances of a variants which ever appears > 1x dup_vars <- X[duplicated(snp), unique(snp)] X <- X[!(snp %in% dup_vars)] cat("Removing ", length(dup_vars), " duplicated variants leaving ", nrow(X), "variants.\n") @@ -137,6 +143,21 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, align_beta(X) # --- write out results --- + if (return_og_snps){ + # label all snps that made it thru filtering as kept + X[, pass_filt := TRUE] + # make matrix with all snps given to function + X_full <- X[ + original_index, + on = .(row_id, snp) + ] + # make NAs in kept from dropped rows into false + X_full[is.na(pass_filt), pass_filt := FALSE] + + X <- X_full + + print(head(X)) + } if(!missing(output_file)){ cat("Writing out ", nrow(X), " variants to file.\n") # changed from path= to file= @@ -145,7 +166,6 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, } cat("Returning ", nrow(X), " variants.\n") return(X) - } From 5d98b3053cb7a5c50408f5c750a2bf19482ec104 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Wed, 24 Jun 2026 11:32:29 -0400 Subject: [PATCH 04/10] fixed missing closing brace --- R/.gwas_format.R.swp | Bin 24576 -> 0 bytes R/gwas_format.R | 13 +++++++------ man/gwas_format.Rd | 3 ++- 3 files changed, 9 insertions(+), 7 deletions(-) delete mode 100644 R/.gwas_format.R.swp diff --git a/R/.gwas_format.R.swp b/R/.gwas_format.R.swp deleted file mode 100644 index 532afdb164c9a767c5b1e7f1f8dfd5aa274dbdc3..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 24576 zcmeHP--{$i9dAuyyu_SQQBY96nY~~S%#Z9U5|5lZ?q-jb&2hQgGYcwnz1`I_eb>{~ z*Hyi9o8&NvAQ(h^(~A-k-^EuE6_Phm!8}BWM$!BOi6-}ui2C_fbmjIFsE)v8xcl5pG;QMs6{pT(-L%VJOVRX*-s+Mmdx=%3yk3@0*D{H1@cKbi(= z7|UW3=-#=-{4V4ea17ie1Er3;_uqBl6Kkucsjsx|p%1+O!c7Wu*&G9o0mp!2z%k$$ za11yG90QJl|3wCLew+7I82--6_)&F#$I|__DqN>amY&~Qg>R|u@zV2uR@XrlKVEuX z=exa%|MSxGD(G&G0mp!2z%k$$a11yG90QI4$ADwNG2j?*4BP|*LF9QW$oG3*@ZA#eg%0d5EW{%+5E3wROuI`9Qx6W9Pg1{?!k$L7GRz;A#b0Y3n~ z2YeUEfQNuj0v`e12i$lU>;ql{egRwu4uBn?4;%;X0FD8_!v@Dsfgb`-1K$U*HR62; z_#kjA@G>?zo(7%*t^?l&A|M1l4IBslhK-CDffs<^0zU(u1)c$(0L}pq16|-A;7;HM zHdS5&UIm^9o&%l*o&?SSUEnR$&zk_()31SF0ck2mK}rGcM7#8Wmkm=9YkHtd_EK4{ zhEbVwj4CJo?un?7`7lm{QQPy7=Sv-$#n8!9B=&>U4`mi5I*}P?ZFgwzfKCxaeFZrd zp0U^~rPRXjouSXrY8wewj6fP}Z&Q*f9b}?QJf65|f{bY^2AH1Z{c@k7gN+&j=kQaqL&60iAV0 z8dY7HO_|gvmxx4k2qRkNNsJUeNb_qUbZgy=i6Y_e6F>)itGB| zLWe%9p4n((jEoND614-z>%dQdH59Zm4bcM6tyu`njp|f>Vob6HWN*GKL|)e!O$ANS zG$_a-xVdo#_G}{j*~)Abi*0NLu`L~1eV5>$L2T-GFf2rn4s9L#lSHW`8!_i_RBHe| z{x~J7mjw&4?Z@VR)3DvBYJ?YcHXbZ1za~;NkRYg0BthIZ>v#-oj3+diylc#?f!1TVRf zdD4(*!z>e-tyG0}?zc<}S!@~HUf8?Q>2eQ^)>$Tnh-{zC9n!oKvm838uG3mfGfs(S zXt25vTC=5G;jsaQL03XJO<2pr@R64ogQ|f8J@v*8#t9U_W|#z*1x0}ldOFy{uv(>~ zLoZos2ph&pinbrh2?Ao6qh~@7&yHf~J~U_3^A|b@?u-t<4BJ%~RQXBAYaEv3x{S*# z_NxZp9*q^?JhW)a^` zFUt1mnj6DP6F*I887gh`3*$3XvHm z^i*1t;Vp`M(W6n12K{yT63NWE3y=97+UcL>NTyHXqiA!vO%W}P%q%#Lq%jF6Trm%- zH5pQB4I{`mqMGt;1G0J-Rm}4Cs4nwzd94+by<}9%Qq^-8Ld%$%c@qn)jr31L37!o4 zJ8j%ihD=}@r-I_3LrFy{Gi6Gu;6OaXvks$Xp2QLtoT&CJR2fg>Xb&8*Uql!~lT=w3 zt6kXEoGy&7StSY|yN)?ZMx|6J%@x=1 zH1u3YWwzd5qLWN$lyG$-6m*!=%nBxP#7vRrmYVy zV8ViYI??dxb(vvJrH{-49|~B6ExPq#*u?A+ z?P1?4lY0`|vHbGL>=7wpq+?G{@)j3wVZ$7|wOwt$Ry*2uzy6w8{|~XA{|2w=vHr*F zXmdS<@bx;rigP!|fMdWh;23ZWI0hU8jseGjW56-s7;p?Y2L3w(hac>!Vpq0s(yBhr zjbp`pAk&;{-o$Bb{Byk#A8W?h0&}37&o$$4_VQV2b1J}`{Na-)sR;I)`S@%&GS)|+ z?Abeh|Nj%L^Ixvk|2S%Hu0K|Bf581E;FkdJ349H>2Dr;H;23ZWI0hU8jseGjW56-s z7;p?Y1{?#9fw#*5p9=8!V_{WZg^zr|Kcgo6z}fx^pkD7UCr229%#K0*tMKYq(S%=v JR~rf5zW_7pphf@y diff --git a/R/gwas_format.R b/R/gwas_format.R index fedc85d..65a6bbd 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -38,12 +38,12 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, chrom <- "chrom" } if(missing(pos) || is.na(pos)){ - X <- [, pos := NA_integer_] + X[, pos := NA_integer_] pos <- "pos" } if(missing(p_value) || is.na(p_value)){ - X <- [, p_value := NA_real_] + X[, p_value := NA_real_] p_value <- "p_value" p_val_missing <- TRUE }else{ @@ -51,15 +51,15 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, } if(missing(sample_size) || is.na(sample_size)){ - X <- [, sample_size := NA_real_] + X[, sample_size := NA_real_] sample_size <- "sample_size" }else if(is.numeric(sample_size)){ - X <- [, sample_size := NA_real_] + X[, sample_size := NA_real_] sample_size <- "sample_size" } if(missing(allele_freq) || is.na(allele_freq)){ - X <- [, af := NA_real_] + X[, af := NA_real_] allele_freq <- "af" } @@ -137,7 +137,8 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, NA_real_ )] cat("Computed p-value\n") - + } + # --- harmonize alleles --- cat("Flipping strand and effect allele so A1 is always A\n") align_beta(X) diff --git a/man/gwas_format.Rd b/man/gwas_format.Rd index 3f802ba..4996770 100644 --- a/man/gwas_format.Rd +++ b/man/gwas_format.Rd @@ -17,7 +17,8 @@ gwas_format( sample_size, allele_freq, output_file, - compute_pval = TRUE + compute_pval = TRUE, + return_og_snps = FALSE ) } \arguments{ From 9a8d5fb7e98d83dceaebbd14038031ae9f7546d0 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Wed, 24 Jun 2026 11:45:34 -0400 Subject: [PATCH 05/10] fixed var name --- R/gwas_format.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/gwas_format.R b/R/gwas_format.R index 65a6bbd..46bf0bf 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -149,7 +149,7 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, X[, pass_filt := TRUE] # make matrix with all snps given to function X_full <- X[ - original_index, + og_snp_index, on = .(row_id, snp) ] # make NAs in kept from dropped rows into false From 9aacd51acc0c3474944c83d97c63d769f4f8efa6 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Wed, 24 Jun 2026 11:56:17 -0400 Subject: [PATCH 06/10] added print statements for debugging --- R/gwas_format.R | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/R/gwas_format.R b/R/gwas_format.R index 46bf0bf..ee9df70 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -104,7 +104,7 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, # drop ALL instances of a variants which ever appears > 1x dup_vars <- X[duplicated(snp), unique(snp)] X <- X[!(snp %in% dup_vars)] - cat("Removing ", length(dup_vars), " duplicated variants leaving ", nrow(X), "variants.\n") + cat("Removing", length(dup_vars), "duplicated variants leaving", nrow(X), "variants.\n") #Illegal alleles valid_alleles <- c("A", "C", "T", "G") @@ -114,7 +114,8 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, ] if(length(illegal_vars)){ X <- X[!(snp %in% illegal_vars)] - cat("Removing ", length(illegal_vars), " variants with illegal alleles leaving ", nrow(X), "variants.\n") + print(head(X[snp %in% illegal_vars])) + cat("Removing", length(illegal_vars), "variants with illegal alleles leaving", nrow(X), "variants.\n") }else{ cat("No variants have illegal alleles.\n") } @@ -123,7 +124,7 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, # when filtering rows, need to assign to var to keep result n <- nrow(X) X <- remove_ambiguous(X) - cat("Removed ", n-nrow(X), " variants with ambiguous strand.\n") + cat("Removing", n-nrow(X), " variants with ambiguous strand leaving", nrow(X), "variants.\n") # --- compute pval ---- # ask jean do we always want to compute even if not missing @@ -185,6 +186,8 @@ remove_ambiguous <- function(X) { idx_ambig <- X[ambig_pairs, on = .(A1, A2), which = TRUE, nomatch = NULL] + print(head(X[idx_ambig])) + if (!length(idx_ambig)) { return(invisible(X)) } else { From c835116cc4f4a4cd1aca3550386944da6e6a4c83 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Wed, 24 Jun 2026 12:03:06 -0400 Subject: [PATCH 07/10] corrected print statement --- R/gwas_format.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/gwas_format.R b/R/gwas_format.R index ee9df70..8eb2be9 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -113,8 +113,8 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, snp ] if(length(illegal_vars)){ - X <- X[!(snp %in% illegal_vars)] print(head(X[snp %in% illegal_vars])) + X <- X[!(snp %in% illegal_vars)] cat("Removing", length(illegal_vars), "variants with illegal alleles leaving", nrow(X), "variants.\n") }else{ cat("No variants have illegal alleles.\n") From 809f7bfcb55c8d890115bbe932bddc645aee87ab Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Wed, 24 Jun 2026 12:58:10 -0400 Subject: [PATCH 08/10] take away pass_filt, not needed bc estimate_R handles NAs --- R/gwas_format.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/gwas_format.R b/R/gwas_format.R index 8eb2be9..a780636 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -147,14 +147,14 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, # --- write out results --- if (return_og_snps){ # label all snps that made it thru filtering as kept - X[, pass_filt := TRUE] + #X[, pass_filt := TRUE] # make matrix with all snps given to function X_full <- X[ og_snp_index, on = .(row_id, snp) ] # make NAs in kept from dropped rows into false - X_full[is.na(pass_filt), pass_filt := FALSE] + #X_full[is.na(pass_filt), pass_filt := FALSE] X <- X_full From 7162580bc10c1db144815344c57138f2a5f83838 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Wed, 24 Jun 2026 14:59:08 -0400 Subject: [PATCH 09/10] correct matrix making in estimate_R if N is vector --- R/estimate_R.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/estimate_R.R b/R/estimate_R.R index dd0e07c..02c43a8 100644 --- a/R/estimate_R.R +++ b/R/estimate_R.R @@ -67,7 +67,7 @@ R_ldsc <- function(Z_hat, stopifnot(nrow(N) == J & identical(colnames(N),colnames(Z_hat))) }else if(class(N) == "numeric"){ stopifnot(identical(names(N),colnames(Z_hat))) - N <- matrix(rep(N, each = J), nrow = J) + N <- matrix(rep(N, each = J), nrow = J, ncol = M, dimnames = dimnames(Z_hat)) } if(is.null(comparisons)){ From 6006811bf192a6e06b3ee0768c36725694d285e1 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Wed, 15 Jul 2026 09:41:05 -0400 Subject: [PATCH 10/10] reduce user annoyance: remove debugging print statements and properly document gwas_format.R --- NAMESPACE | 11 ++++++++++- R/gwas_format.R | 12 ++++-------- R/package.R | 2 +- man/gwas_format.Rd | 6 +++++- 4 files changed, 20 insertions(+), 11 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 485f9f0..b02be59 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,7 +17,6 @@ export(min_norm) export(plot_factors) export(plot_factors_bars) export(snp_ldsc) -import(data.table) import(dplyr) import(flashier) import(ggplot2) @@ -25,3 +24,13 @@ import(lpSolve) import(purrr) import(readr) import(reshape2) +importFrom(data.table,"%chin%") +importFrom(data.table,":=") +importFrom(data.table,.I) +importFrom(data.table,data.table) +importFrom(data.table,fifelse) +importFrom(data.table,fwrite) +importFrom(data.table,is.data.table) +importFrom(data.table,setDT) +importFrom(data.table,setcolorder) +importFrom(data.table,setnames) diff --git a/R/gwas_format.R b/R/gwas_format.R index a780636..e43d308 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -10,14 +10,16 @@ #'@param pos Position column (optional) #'@param p_value p-value column (optional) #'@param sample_size Sample size column (optional) or an integer +#'@param allele_freq Effect allele frequency column (optional) #'@param compute_pval Logical, compute the p-value using a normal approximation if missing? Defaults to TRUE. #'@param output_file File to write out formatted data. If missing formatted data will be returned. +#'@param return_og_snps Option to return all input SNPs in original order. SNPs filtered out will have info set to NA #'@details This function will try to merge data sets X1 and X2 on the specified columns. Where #'necessary, it will flip the sign of effects so that the effect allele is the same in both #'data sets. It will remove variants with ambiguous alleles or where the alleles (G/C or A/T) or #'with alleles that do not match between data sets (e.g A/G in one data set and A/C in the other). #'It will not remove variants that are simply strand flipped between the two data sets (e. g. A/C in one data set, T/G in the other). -#'@return A data frame with columns chrom, pos, snp, A1, A2, beta_hat, se, p_value, and sample_size with all SNPs +#'@return A data frame with columns chrom, pos, snp, A1, A2, beta_hat, se, p_value, sample_size, and allele_freq with all SNPs #'aligned so that A is the effect allele. This is ready to be used with gwas_merge with formatted = TRUE. #'@export gwas_format <- function(X, snp, beta_hat, se, A1, A2, @@ -113,7 +115,6 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, snp ] if(length(illegal_vars)){ - print(head(X[snp %in% illegal_vars])) X <- X[!(snp %in% illegal_vars)] cat("Removing", length(illegal_vars), "variants with illegal alleles leaving", nrow(X), "variants.\n") }else{ @@ -157,13 +158,10 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, #X_full[is.na(pass_filt), pass_filt := FALSE] X <- X_full - - print(head(X)) } if(!missing(output_file)){ cat("Writing out ", nrow(X), " variants to file.\n") - # changed from path= to file= - write_tsv(X, file = output_file) + fwrite(X, file = output_file, sep="\t", na = "NA") return(NULL) } cat("Returning ", nrow(X), " variants.\n") @@ -186,8 +184,6 @@ remove_ambiguous <- function(X) { idx_ambig <- X[ambig_pairs, on = .(A1, A2), which = TRUE, nomatch = NULL] - print(head(X[idx_ambig])) - if (!length(idx_ambig)) { return(invisible(X)) } else { diff --git a/R/package.R b/R/package.R index d121965..af36761 100644 --- a/R/package.R +++ b/R/package.R @@ -7,6 +7,6 @@ #' @import flashier #' @import dplyr readr reshape2 lpSolve #' @import purrr ggplot2 -#' @import data.table +#' @importFrom data.table setDT setcolorder setnames data.table is.data.table fifelse fwrite %chin% .I := #' @name GFA NULL diff --git a/man/gwas_format.Rd b/man/gwas_format.Rd index 4996770..6759e7d 100644 --- a/man/gwas_format.Rd +++ b/man/gwas_format.Rd @@ -42,12 +42,16 @@ gwas_format( \item{sample_size}{Sample size column (optional) or an integer} +\item{allele_freq}{Effect allele frequency column (optional)} + \item{output_file}{File to write out formatted data. If missing formatted data will be returned.} \item{compute_pval}{Logical, compute the p-value using a normal approximation if missing? Defaults to TRUE.} + +\item{return_og_snps}{Option to return all input SNPs in original order. SNPs filtered out will have info set to NA} } \value{ -A data frame with columns chrom, pos, snp, A1, A2, beta_hat, se, p_value, and sample_size with all SNPs +A data frame with columns chrom, pos, snp, A1, A2, beta_hat, se, p_value, sample_size, and allele_freq with all SNPs aligned so that A is the effect allele. This is ready to be used with gwas_merge with formatted = TRUE. } \description{