Mattbranch

R/locuszoom.R

@@ -24,7 +24,12 @@ optional <- c("flanking_region"=500,
               "track_file"=NA,
               "track_file_type"="window",
               "track_label"="",
-              "track_threshold"=5e-8)
+              "track_threshold"=5e-8,
+              "signif_level"=5e-8,
+              "gene_rows"=4,
+              "title"=TRUE,
+              "variant_label_file"=NA,
+              "variant_line_color"="transparent")
 config <- setConfigDefaults(config, required, optional)
 print(config)
 
@@ -47,28 +52,44 @@ assoc <- getobj(assocfile)
 if (config["locus_type"] == "variant") {
     stopifnot("variant.id" %in% names(locus))
     variant <- locus$variant.id
-    flank <- as.numeric(config["flanking_region"]) * 1000
     var.pos <- assoc$pos[assoc$variant.id == variant]
-    start <- var.pos - flank
-    if (start < 1) start <- 1
-    end <- var.pos + flank
-
-    lz.name <- paste0("chr", var.chr, ":", var.pos)
-    ld.region <- paste0("--refsnp \"", lz.name, "\"", " --flank ", config["flanking_region"], "kb")
-    prefix <- paste0(config["out_prefix"], "_var", variant, "_ld_", pop)
+
+    if("rsid" %in% names(locus) && !is.na(locus$rsid)){
+      lz.name <- locus$rsid
+    }else{
+      lz.name <- paste0("chr", var.chr, ":", var.pos)
+    }
+
+    if(all(c("variant.id", "start", "end") %in% names(locus))){
+      start <- locus$start
+      end <- locus$end
+      ld.region <- paste("--refsnp", lz.name, "--chr", var.chr, "--start", start, "--end", end)
+    }else{
+      flank <- as.numeric(config["flanking_region"]) * 1000
+      start <- var.pos - flank
+      if (start < 1) start <- 1
+      end <- var.pos + flank
+      ld.region <- paste0("--refsnp \"", lz.name, "\"", " --flank ", config["flanking_region"], "kb")
+    }   
+
     freq <- assoc$freq[assoc$variant.id == variant]
     maf <- min(freq, 1-freq)
     mac <- assoc$MAC[assoc$variant.id == variant]
-    title <- paste(lz.name, "- MAF:", formatC(maf, digits=3), "- MAC:", mac)
+
+    prefix <- paste0(config["out_prefix"], "_var", variant, "_ld_", pop)
 
 } else if (config["locus_type"] == "region") {
     stopifnot(all(c("start", "end") %in% names(locus)))
     start <- locus$start
     end <- locus$end
 
     ld.region <- paste("--chr", var.chr, "--start", start, "--end", end)
+
     prefix <- paste0(config["out_prefix"], "_ld_", pop)
-    title <- ""
+}
+
+if("locus_name" %in% names(locus)){
+  prefix <- paste(prefix, locus$locus_name, sep = "_")
 }
 
 ## construct METAL-format file
@@ -77,20 +98,18 @@ assoc <- assoc %>%
     select(variant.id, chr, pos, ends_with("pval"))
 names(assoc)[4] <- "pval"
 
-##### NOTE: i added the following two lines of code #####
 ## remove duplicate chr:pos rows, removing the variant with the less significant (higher) pvalue
 assoc <- assoc[order(assoc$pval,decreasing=FALSE),]
 assoc <- assoc[!duplicated(assoc$pos),]
 assoc <- assoc[order(assoc$pos),]
-#####
 
 assoc.filename <- tempfile()
 writeMETAL(assoc, file=assoc.filename)
 
+
 # LD
 if (pop != "TOPMED") {
     ld.cmd <- paste("--pop", pop, "--source 1000G_Nov2014")
-    ld.title <- paste("LD: 1000G", pop)
 } else {
     if (!is.na(config["ld_sample_include"])) {
         sample.id <- getobj(config["ld_sample_include"])
@@ -108,13 +127,16 @@ if (pop != "TOPMED") {
     }
     gdsfile <- insertChromString(config["gds_file"], var.chr)
     ld <- calculateLD(gdsfile, variant.id=assoc$variant.id, ref.var=ref.var, sample.id=sample.id)
+
     ld.filename <- tempfile()
-    writeLD(assoc, ld, ref.var, file=ld.filename)
+    if("rsid" %in% names(locus) && !is.na(locus$rsid)){
+      writeLD(assoc, ld, ref.var, file=ld.filename, rsid = locus$rsid)
+    }else{
+      writeLD(assoc, ld, ref.var, file=ld.filename, rsid = NULL)
+    }
 
     ld.cmd <- paste("--ld", ld.filename)
-    ld.title <- "LD: TOPMed"
 }
-title <- if (title == "") ld.title else paste(ld.title, title, sep=" - ")
 
 ## construct BED track file
 if (!is.na(config["track_file"])) {
@@ -128,6 +150,28 @@ if (!is.na(config["track_file"])) {
     track.cmd <- ""
 }
 
+
+## Plot Title
+if(config["title"]){
+  if(config["locus_type"] == "variant"){
+    title <- ifelse("locus_name" %in% names(locus), locus$locus_name, lz.name)
+    maf.title <- paste("- MAF:", formatC(maf, digits=3), "- MAC:", mac)
+  }else if(config["locus_type"] == "region"){
+    title <- ifelse("locus_name" %in% names(locus), locus$locus_name, "")
+    maf.title <- ""
+  }
+  if(pop != "TOPMED"){
+    ld.title <- paste("- LD: 1000G", pop)
+  }else{
+    ld.title <- "- LD: TOPMed"
+  }
+  title <- paste(title, ld.title, maf.title)
+}else{
+  title <- NULL
+}
+
+
+
 command <- paste("locuszoom",
                  "theme=publication",
                  "--cache None",
@@ -141,9 +185,22 @@ command <- paste("locuszoom",
                  ld.cmd,
                  ld.region,
                  "--prefix ", prefix,
-                 paste0("title=\"", title, "\""),
-                 paste0("signifLine=\"", -log10(5e-8), "\" signifLineColor=\"gray\" signifLineWidth=\"2\""),
-                 "ylab=\"-log10(p-value) from single variant test\"")
+                 paste0("signifLine=\"", -log10(as.numeric(config["signif_level"])), "\""), 
+                 "signifLineColor=\"gray\"", 
+                 "signifLineWidth=\"2\"",
+                 "ylab=\"-log10(p-value)\"",
+                 paste0("refsnpLineColor=\"", config["variant_line_color"], "\""),
+                 "refsnpLineWidth=\"2\"",
+                 "refsnpTextSize=\"1\"",
+                 paste0("rfrows=\"", as.numeric(config["gene_rows"]), "\""))
+
+if(!is.null(title)){
+  command <- paste(command, paste0("title=\"", title, "\""))
+}
+
+if(!is.na(config["variant_label_file"])){
+  command <- paste(command, "--denote-markers-file", config["variant_label_file"])
+}
 
 cat(paste(command, "\n"))
 system(command)

TopmedPipeline/R/locuszoom.R

@@ -63,8 +63,20 @@ calculateLD <- function(gdsfile, variant.id, ref.var=NULL, sample.id=NULL) {
         geno.ref <- geno[,as.character(ref.var)]
     }
 
-    # drop unnecessary dimensions (if ref.var is a single variant)
-    r <- drop(cor(geno, geno.ref, use="pairwise.complete.obs"))
+    # check to avoid long vectors
+    nr <- as.numeric(nrow(geno))
+    nc <- as.numeric(ncol(geno))
+    nblock <- ceiling(nr*nc/2^31)
+    if(nblock > 1){
+        # break into blocks to compute correlation if needed
+        blocks <- unname(split(1:nc, cut(1:nc, nblock)))
+        r <- unlist(lapply(blocks, function(b) {
+            drop(cor(geno[,b], geno.ref, use="pairwise.complete.obs"))
+        }))
+    }else{
+        # drop unnecessary dimensions (if ref.var is a single variant)
+        r <- drop(cor(geno, geno.ref, use="pairwise.complete.obs"))
+    }
     r^2
 }
 
@@ -78,10 +90,17 @@ calculateLD <- function(gdsfile, variant.id, ref.var=NULL, sample.id=NULL) {
 #' 
 #' @importFrom dplyr "%>%" mutate_ select_
 #' @export
-writeLD <- function(x, ld, ref.var, file) {
+writeLD <- function(x, ld, ref.var, file, rsid = NULL) {
     x <- x %>%
         mutate_(snp1=~paste0("chr", chr, ":", pos))
-    x$snp2 <- x$snp1[x$variant.id == ref.var]
+
+    if(!is.null(rsid)){
+        x$snp1[x$variant.id == ref.var] <- rsid
+        x$snp2 <- rsid
+    }else{
+        x$snp2 <- x$snp1[x$variant.id == ref.var]
+    }
+
     x$dprime <- 0
     x$rsquare <- ld
     x <- x %>%