From 0dff9806320a0376a26c4b9859debd8b12e10059 Mon Sep 17 00:00:00 2001 From: "jinghuazhao@github.com" Date: Thu, 19 Dec 2024 18:01:17 +0000 Subject: [PATCH] Deployed 87a0ee6 with MkDocs version: 1.5.3 --- index.html | 2 +- misc/coloc.R | 29 +++++++++++++----------- misc/json.sh | 60 ++++++++++++++++++++++++++++++++++++++++++++++++- sitemap.xml.gz | Bin 127 -> 127 bytes 4 files changed, 76 insertions(+), 15 deletions(-) diff --git a/index.html b/index.html index 26b44a2..37f206e 100644 --- a/index.html +++ b/index.html @@ -296,5 +296,5 @@ diff --git a/misc/coloc.R b/misc/coloc.R index 336805a..579a7d9 100644 --- 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 index cc49d6c..dac561a 100644 --- 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 index 1868238b339dd598ca432a159616567db3b611b6..84b08e628d29192f8e242b982bf1d7a37667281a 100644 GIT binary patch delta 13 Ucmb=gXP58h;K+?nnaExN02_n^00000 delta 13 Ucmb=gXP58h;1Kjlp2%JS02uiLkN^Mx