Deployed 87a0ee6 with MkDocs version: 1.5.3

jinghuazhao · Dec 19, 2024 · 0dff980 · 0dff980
1 parent 0c7cef6
commit 0dff980
Show file tree

Hide file tree

Showing 4 changed files with 76 additions and 15 deletions.
diff --git a/index.html b/index.html
@@ -296,5 +296,5 @@ <h4 class="modal-title" id="keyboardModalLabel">Keyboard Shortcuts</h4>
 
 <!--
 MkDocs version : 1.5.3
-Build Date UTC : 2024-12-18 22:26:25.266596+00:00
+Build Date UTC : 2024-12-19 18:01:17.146975+00:00
 -->
diff --git a/misc/coloc.R b/misc/coloc.R
@@ -85,7 +85,6 @@ gtex <- function(gwas_stats_hg38,ensGene,region38)
   invisible(sapply(1:49, function(i) {
     f <- file.path(analysis, "coloc", "GTEx", "sumstats", paste0(prot, "-", names(result_filtered)[i], ".gz"))
     write.table(result_filtered[[i]], file = gzfile(f, "w"), col.names = TRUE, row.names = FALSE, quote = FALSE, sep = "\t")
-    close(gz_connection)
   }))
   purrr::map_df(result_filtered, ~run_coloc(., gwas_stats_hg38), .id = "qtl_id")
 }
@@ -106,10 +105,9 @@ ge <- function(gwas_stats_hg38,ensGene,region38)
   result_list <- result_list[!unlist(purrr::map(result_list, is.null))]
   result_filtered <- purrr::map(result_list[lapply(result_list,nrow)!=0],
                                 ~dplyr::filter(., !is.na(se)))
-  invisible(sapply(1:49, function(i) {
+  invisible(sapply(1:length(result_filtered), function(i) {
     f <- file.path(analysis, "coloc", "eQTLCatalogue", "sumstats", paste0(prot, "-", names(result_filtered)[i], ".gz"))
     write.table(result_filtered[[i]], file = gzfile(f, "w"), col.names = TRUE, row.names = FALSE, quote = FALSE, sep = "\t")
-    close(gz_connection)
   }))
   purrr::map_df(result_filtered, ~run_coloc(., gwas_stats_hg38), .id = "unique_id")
 }
@@ -166,27 +164,25 @@ all_coloc <- function(prot,chr,ensGene,region37,region38,out)
 
 single_run <- function(r, batch="GTEx")
 {
-  sentinel <- sentinels[r,]
-  chr <- with(sentinel,geneChrom)
+  ss <- subset(pQTLdata::caprion,Protein==paste0(prot,"_HUMAN"))
+  ensGene <- ss[["ensGenes"]]
   ensRegion37 <- with(sentinel,
                       {
                         start <- geneStart-M
                         if (start<0) start <- 0
                         end <- geneEnd+M
-                        paste0(chr,":",start,"-",end)
+                        paste0(geneChrom,":",start,"-",end)
                       })
-  ss <- subset(pQTLdata::caprion,Protein==paste0(sentinel[["prot"]],"_HUMAN"))
-  ensGene <- ss[["ensGenes"]]
-  x <- with(sentinel,list(chr=geneChrom,start=geneStart,end=geneEnd))
+  x <- list(chr=geneChrom,start=geneStart,end=geneEnd)
   lr <- liftRegion(x)
   ensRegion38 <- with(lr,paste0(chr,":",start-M,"-",end+M))
-  cat(chr,ensGene,ensRegion37,ensRegion38,"\n")
-  f <- file.path(analysis,"coloc",batch,with(sentinel,paste0(prot,"-",SNP)))
+  cat(geneChrom,ensGene,ensRegion37,ensRegion38,"\n")
+  f <- file.path(analysis,"coloc",batch,paste0(prot,"-",snp))
   if (batch=="GTEx")
   {
-    gtex_coloc(sentinel[["prot"]],sentinel[["SNP"]],chr,ensGene,ensRegion37,ensRegion38,f)
+    gtex_coloc(prot,snp,geneChrom,ensGene,ensRegion37,ensRegion38,f)
   } else {
-    ge_coloc(sentinel[["prot"]],chr,ensGene,ensRegion37,ensRegion38,f)
+    ge_coloc(prot,geneChrom,ensGene,ensRegion37,ensRegion38,f)
   }
 }
 
@@ -282,6 +278,13 @@ f <- file.path(find.package("pQTLtools"),"eQTL-Catalogue","hg19ToHg38.over.chain
 chain <- rtracklayer::import.chain(f)
 
 r <- as.integer(Sys.getenv("r"))
+sentinel <- sentinels[r,]
+prot <- sentinel[["prot"]]
+snp <- sentinel[["SNP"]]
+geneChrom <- sentinel[["geneChrom"]]
+geneStart <- sentinel[["geneStart"]]
+geneEnd <- sentinel[["geneEnd"]]
+
 single_run(r)
 single_run(r,batch="eQTLCatalogue")
 

diff --git a/misc/json.sh b/misc/json.sh
@@ -5,7 +5,7 @@ export analysis=~/Caprion/analysis
 export pre_qc_data=/rds/project/rds-MkfvQMuSUxk/interval/caprion_proteomics
 export suffix=_dr
 
-module load ceuadmin/htslib
+module load ceuadmin/htslib ceuadmin/R
 
 function gz()
 {
@@ -83,6 +83,64 @@ function lz_json()
   ' | gzip -f > ${analysis}/json/gz/top_hits.json.gz
 }
 
+function gtex()
+{
+  awk 'NR>1{print $1,$2,$3,$4}' ${analysis}/coloc/GTEx.tsv | \
+  parallel -C ' ' '
+  export prot={1}
+  export rsid={2}
+  export snpid={3}
+  export tissue={4}
+  Rscript -e "
+    suppressMessages(library(dplyr))
+    suppressMessages(library(jsonlite))
+    analysis <- Sys.getenv(\"analysis\")
+    prot <- Sys.getenv(\"prot\")
+    rsid <- Sys.getenv(\"rsid\")
+    snpid <- Sys.getenv(\"snpid\")
+    tissue <- Sys.getenv(\"tissue\")
+    if (!dir.exists(file.path(analysis,\"json\",\"GTEx\"))) dir.create(file.path(analysis,\"json\",\"GTEx\"))
+    f <- file.path(analysis,\"coloc\",\"GTEx\",\"sumstats\")
+    sumstats <- read.delim(file.path(f,paste0(prot,\"-\",tissue,\".gz\"))) %>%
+                dplyr::mutate(log_pvalue=-log10(pvalue),ref_allele=ref,alt_allele=alt) %>%
+                dplyr::select(chromosome,position,variant,ref_allele,alt_allele,log_pvalue,beta,se)
+    j <- gzfile(file.path(analysis,\"json\",\"GTEx\",paste0(prot,\"-\",tissue,\".json.gz\")))
+    sink(j)
+    print(jsonlite::toJSON(list(ppid=paste0(prot),data=sumstats),auto_unbox=TRUE,pretty=FALSE))
+    sink()
+  "
+  '
+}
+
+function eQTLCatalogue()
+{
+  awk 'NR>1{print $1,$2,$3,$4}' ${analysis}/coloc/eQTLCatalogue.tsv | \
+  parallel -C ' ' '
+  export prot={1}
+  export rsid={2}
+  export snpid={3}
+  export tissue={4}
+  Rscript -e "
+    suppressMessages(library(dplyr))
+    suppressMessages(library(jsonlite))
+    analysis <- Sys.getenv(\"analysis\")
+    prot <- Sys.getenv(\"prot\")
+    rsid <- Sys.getenv(\"rsid\")
+    snpid <- Sys.getenv(\"snpid\")
+    tissue <- Sys.getenv(\"tissue\")
+    if (!dir.exists(file.path(analysis,\"json\",\"eQTLCatalogue\"))) dir.create(file.path(analysis,\"json\",\"eQTLCatalogue\"))
+    f <- file.path(analysis,\"coloc\",\"eQTLCatalogue\",\"sumstats\")
+    sumstats <- read.delim(file.path(f,paste0(prot,\"-\",tissue,\".gz\"))) %>%
+                dplyr::mutate(log_pvalue=-log10(pvalue),ref_allele=ref,alt_allele=alt) %>%
+                dplyr::select(chromosome,position,variant,ref_allele,alt_allele,log_pvalue,beta,se)
+    j <- gzfile(file.path(analysis,\"json\",\"eQTLCatalogue\",paste0(prot,\"-\",tissue,\".json.gz\")))
+    sink(j)
+    print(jsonlite::toJSON(list(ppid=paste0(prot),data=sumstats),auto_unbox=TRUE,pretty=FALSE))
+    sink()
+  "
+  '
+}
+
 function umich()
 {
   curl https://portaldev.sph.umich.edu/api/v1/statistic/single/ > single.json

diff --git a/sitemap.xml.gz b/sitemap.xml.gz