Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 32 additions & 30 deletions code/SoS/association_scan/TensorQTL/TensorQTL.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@
"data": {
"text/html": [
"<table class=\"dataframe\">\n",
"<caption>A data.table: 1 \u00d7 2</caption>\n",
"<caption>A data.table: 1 × 2</caption>\n",
"<thead>\n",
"\t<tr><th scope=col>#id</th><th scope=col>#dir</th></tr>\n",
"\t<tr><th scope=col>&lt;int&gt;</th><th scope=col>&lt;chr&gt;</th></tr>\n",
Expand All @@ -67,7 +67,7 @@
"</table>\n"
],
"text/latex": [
"A data.table: 1 \u00d7 2\n",
"A data.table: 1 × 2\n",
"\\begin{tabular}{ll}\n",
" \\#id & \\#dir\\\\\n",
" <int> & <chr>\\\\\n",
Expand All @@ -77,7 +77,7 @@
],
"text/markdown": [
"\n",
"A data.table: 1 \u00d7 2\n",
"A data.table: 1 × 2\n",
"\n",
"| #id &lt;int&gt; | #dir &lt;chr&gt; |\n",
"|---|---|\n",
Expand All @@ -96,7 +96,7 @@
"data": {
"text/html": [
"<table class=\"dataframe\">\n",
"<caption>A data.table: 5 \u00d7 8</caption>\n",
"<caption>A data.table: 5 × 8</caption>\n",
"<thead>\n",
"\t<tr><th scope=col>#chr</th><th scope=col>start</th><th scope=col>end</th><th scope=col>ID</th><th scope=col>SAMPLE_001</th><th scope=col>SAMPLE_002</th><th scope=col>SAMPLE_003</th><th scope=col>SAMPLE_004</th></tr>\n",
"\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
Expand All @@ -111,7 +111,7 @@
"</table>\n"
],
"text/latex": [
"A data.table: 5 \u00d7 8\n",
"A data.table: 5 × 8\n",
"\\begin{tabular}{llllllll}\n",
" \\#chr & start & end & ID & SAMPLE\\_001 & SAMPLE\\_002 & SAMPLE\\_003 & SAMPLE\\_004\\\\\n",
" <chr> & <int> & <int> & <chr> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
Expand All @@ -125,7 +125,7 @@
],
"text/markdown": [
"\n",
"A data.table: 5 \u00d7 8\n",
"A data.table: 5 × 8\n",
"\n",
"| #chr &lt;chr&gt; | start &lt;int&gt; | end &lt;int&gt; | ID &lt;chr&gt; | SAMPLE_001 &lt;dbl&gt; | SAMPLE_002 &lt;dbl&gt; | SAMPLE_003 &lt;dbl&gt; | SAMPLE_004 &lt;dbl&gt; |\n",
"|---|---|---|---|---|---|---|---|\n",
Expand Down Expand Up @@ -192,7 +192,7 @@
"data": {
"text/html": [
"<table class=\"dataframe\">\n",
"<caption>A data.table: 6 \u00d7 2</caption>\n",
"<caption>A data.table: 6 × 2</caption>\n",
"<thead>\n",
"\t<tr><th scope=col>#id</th><th scope=col>#path</th></tr>\n",
"\t<tr><th scope=col>&lt;int&gt;</th><th scope=col>&lt;chr&gt;</th></tr>\n",
Expand All @@ -208,7 +208,7 @@
"</table>\n"
],
"text/latex": [
"A data.table: 6 \u00d7 2\n",
"A data.table: 6 × 2\n",
"\\begin{tabular}{ll}\n",
" \\#id & \\#path\\\\\n",
" <int> & <chr>\\\\\n",
Expand All @@ -223,7 +223,7 @@
],
"text/markdown": [
"\n",
"A data.table: 6 \u00d7 2\n",
"A data.table: 6 × 2\n",
"\n",
"| #id &lt;int&gt; | #path &lt;chr&gt; |\n",
"|---|---|\n",
Expand Down Expand Up @@ -271,7 +271,7 @@
"data": {
"text/html": [
"<table class=\"dataframe\">\n",
"<caption>A matrix: 5 \u00d7 5 of type int</caption>\n",
"<caption>A matrix: 5 × 5 of type int</caption>\n",
"<thead>\n",
"\t<tr><th></th><th scope=col>SAMPLE_001</th><th scope=col>SAMPLE_002</th><th scope=col>SAMPLE_003</th><th scope=col>SAMPLE_004</th><th scope=col>SAMPLE_005</th></tr>\n",
"</thead>\n",
Expand All @@ -285,7 +285,7 @@
"</table>\n"
],
"text/latex": [
"A matrix: 5 \u00d7 5 of type int\n",
"A matrix: 5 × 5 of type int\n",
"\\begin{tabular}{r|lllll}\n",
" & SAMPLE\\_001 & SAMPLE\\_002 & SAMPLE\\_003 & SAMPLE\\_004 & SAMPLE\\_005\\\\\n",
"\\hline\n",
Expand All @@ -298,7 +298,7 @@
],
"text/markdown": [
"\n",
"A matrix: 5 \u00d7 5 of type int\n",
"A matrix: 5 × 5 of type int\n",
"\n",
"| <!--/--> | SAMPLE_001 | SAMPLE_002 | SAMPLE_003 | SAMPLE_004 | SAMPLE_005 |\n",
"|---|---|---|---|---|---|\n",
Expand Down Expand Up @@ -356,7 +356,7 @@
"data": {
"text/html": [
"<table class=\"dataframe\">\n",
"<caption>A data.table: 28 \u00d7 5</caption>\n",
"<caption>A data.table: 28 × 5</caption>\n",
"<thead>\n",
"\t<tr><th scope=col>#id</th><th scope=col>SAMPLE_001</th><th scope=col>SAMPLE_002</th><th scope=col>SAMPLE_003</th><th scope=col>SAMPLE_004</th></tr>\n",
"\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
Expand Down Expand Up @@ -394,7 +394,7 @@
"</table>\n"
],
"text/latex": [
"A data.table: 28 \u00d7 5\n",
"A data.table: 28 × 5\n",
"\\begin{tabular}{lllll}\n",
" \\#id & SAMPLE\\_001 & SAMPLE\\_002 & SAMPLE\\_003 & SAMPLE\\_004\\\\\n",
" <chr> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
Expand Down Expand Up @@ -431,7 +431,7 @@
],
"text/markdown": [
"\n",
"A data.table: 28 \u00d7 5\n",
"A data.table: 28 × 5\n",
"\n",
"| #id &lt;chr&gt; | SAMPLE_001 &lt;dbl&gt; | SAMPLE_002 &lt;dbl&gt; | SAMPLE_003 &lt;dbl&gt; | SAMPLE_004 &lt;dbl&gt; |\n",
"|---|---|---|---|---|\n",
Expand Down Expand Up @@ -559,7 +559,7 @@
"data": {
"text/html": [
"<table class=\"dataframe\">\n",
"<caption>A data.table: 6 \u00d7 4</caption>\n",
"<caption>A data.table: 6 × 4</caption>\n",
"<thead>\n",
"\t<tr><th scope=col>#chr</th><th scope=col>start</th><th scope=col>end</th><th scope=col>gene_id</th></tr>\n",
"\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;chr&gt;</th></tr>\n",
Expand All @@ -575,7 +575,7 @@
"</table>\n"
],
"text/latex": [
"A data.table: 6 \u00d7 4\n",
"A data.table: 6 × 4\n",
"\\begin{tabular}{llll}\n",
" \\#chr & start & end & gene\\_id\\\\\n",
" <chr> & <int> & <int> & <chr>\\\\\n",
Expand All @@ -590,7 +590,7 @@
],
"text/markdown": [
"\n",
"A data.table: 6 \u00d7 4\n",
"A data.table: 6 × 4\n",
"\n",
"| #chr &lt;chr&gt; | start &lt;int&gt; | end &lt;int&gt; | gene_id &lt;chr&gt; |\n",
"|---|---|---|---|\n",
Expand Down Expand Up @@ -737,12 +737,12 @@
"\n",
" **Computational strategy designed for trans analysis:**\n",
" \n",
" Trans analysis faces significant memory challenges as we calculate all associations between all molecular traits \u00d7 all genetic variants across the genome, creating a massive computational burden. To address this challenge, we implement a two-stage chromosome-based parallelization approach:\n",
" Trans analysis faces significant memory challenges as we calculate all associations between all molecular traits × all genetic variants across the genome, creating a massive computational burden. To address this challenge, we implement a two-stage chromosome-based parallelization approach:\n",
"\n",
" **Stage 1 (trans_1): Chromosome-based parallelization**\n",
" - Phenotype data is processed per chromosome (e.g., 22 separate jobs for autosomes)\n",
" - For each phenotype chromosome, we test associations against variants from all 22 chromosomes\n",
" - This creates phenotype_chr \u00d7 genotype_chr combinations (e.g., phenotype chr1 vs genotype chr1-22); Garbage was collected between each chromosome combination caculation to release memory\n",
" - This creates phenotype_chr × genotype_chr combinations (e.g., phenotype chr1 vs genotype chr1-22); Garbage was collected between each chromosome combination caculation to release memory\n",
" - Results are combined across all chromosome combinations and saved as compressed files\n",
"\n",
" **Stage 2 (trans_2): Significance filtering**\n",
Expand Down Expand Up @@ -1117,8 +1117,11 @@
" ${\"--keep-sample \" + str(keep_sample) if keep_sample.is_file() else \"\"} \\\n",
" ${\"--customized-cis-windows \" + str(customized_cis_windows) if customized_cis_windows.is_file() else \"\"} \\\n",
" ${\"--phenotype-group \" + str(phenotype_group) if phenotype_group.is_file() else \"\"} \\\n",
" --skip-postprocess \\\n",
" --numThreads ${numThreads}\n"
" --numThreads ${numThreads}\n",
"\n",
" # Add per-gene Storey q-values to the raw nominal table (in place; re-bgzip+tabix)\n",
" Rscript ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.R \\\n",
" nominal_qvalues \"${_output['nominal']}\" \"${interaction}\"\n"
]
},
{
Expand All @@ -1138,13 +1141,8 @@
"input_files = [str(x) for x in _input[\"regional\"]]\n",
"task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
"bash: expand= \"${ }\", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', container = container, entrypoint = entrypoint\n",
" python3 ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.py \\\n",
" --step cis_postprocess \\\n",
" --cwd \"${cwd}\" \\\n",
" --output-tsv \"${_output[0]}\" \\\n",
" --output-summary \"${_output[1]}\" \\\n",
" --regional-files ${\" \".join(input_files)} \\\n",
" --numThreads ${numThreads}\n"
" Rscript ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.R \\\n",
" regional_postprocess \"${_output[0]}\" \"${_output[1]}\" ${\" \".join(input_files)}\n"
]
},
{
Expand Down Expand Up @@ -1189,7 +1187,11 @@
" ${\"--region-list \" + str(region_list) + \" --region-list-phenotype-column \" + str(region_list_phenotype_column) if region_list.is_file() else \"\"} \\\n",
" ${\"--keep-sample \" + str(keep_sample) if keep_sample.is_file() else \"\"} \\\n",
" ${\"--customized-cis-windows \" + str(customized_cis_windows) if customized_cis_windows.is_file() else \"\"} \\\n",
" --numThreads ${numThreads}\n"
" --numThreads ${numThreads}\n",
"\n",
" # Add global Storey q-values to the raw trans table (in place; re-bgzip+tabix)\n",
" Rscript ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.R \\\n",
" trans_qvalues \"${_output['nominal']}\"\n"
]
},
{
Expand Down
File renamed without changes.
Original file line number Diff line number Diff line change
Expand Up @@ -143,6 +143,7 @@
"outputs": [],
"source": [
"[global]\n",
"parameter: modular_script_dir = path('code/script') # override with --modular-script-dir\n",
"import glob\n",
"import pandas as pd\n",
"## Path to work directory where output locates\n",
Expand Down Expand Up @@ -181,9 +182,6 @@
"import os\n",
"import re\n",
"\n",
"def group_by_region(lst, partition):\n",
" return partition\n",
"\n",
"bim_files = []\n",
"\n",
"for path in geno_list_paths:\n",
Expand All @@ -209,7 +207,7 @@
"\n",
"chr_info = regional_data['chr']\n",
"\n",
"input: regional_data[\"geno_list\"], group_by = lambda x: group_by_region(x, regional_data[\"geno_list\"]), group_with = \"chr_info\"\n",
"input: regional_data[\"geno_list\"], group_by = lambda x: regional_data[\"geno_list\"], group_with = \"chr_info\"\n",
"# a funtion to get geno prefixs as genos\n",
"\n",
"geno_list = regional_data['geno_list'][0]\n",
Expand All @@ -218,43 +216,9 @@
"\n",
"output: f\"{cwd}/{name}.{genos}.{_chr_info}.aligned.bim.gz\", f\"{cwd}/{name}.{genos}.{_chr_info}.aligned.bim.gz.tbi\"\n",
"task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
"R: expand = '${ }', stdout = f\"{_output[0]:nn}.stdout\", stderr = f\"{_output[0]:nn}.stderr\", container = container\n",
" library(tidyverse) \n",
" library(data.table)\n",
" \n",
" merge_genotype_data <- function(df1, df2, all = TRUE) {\n",
" setDT(df1)\n",
" setDT(df2)\n",
" df1[, key := paste(chr, pos, pmin(alt, ref), pmax(alt, ref))]\n",
" df2[, key := paste(chr, pos, pmin(alt, ref), pmax(alt, ref))]\n",
" df2[df1, on = \"key\", flip := i.alt == ref & i.ref == alt, by = .EACHI]\n",
" df2[flip == TRUE, c(\"alt\", \"ref\") := .(ref, alt)]\n",
" if (all) {\n",
" df_combined <- unique(rbindlist(list(df1[, .(chr, pos, alt, ref)], df2[, .(chr, pos, alt, ref)])), by = c(\"chr\", \"pos\", \"alt\", \"ref\"))\n",
" } else {\n",
" df_combined <- df2[, .(chr, pos, alt, ref)]\n",
" }\n",
" return(df_combined)\n",
" }\n",
"\n",
" \n",
" bim_files = c(${\",\".join(['\"%s\"' % x.absolute() for x in _input])})\n",
" rosmap_bim <- fread(bim_files[1], header = FALSE) \n",
" colnames(rosmap_bim) <- c('chr', 'id', 'Mics', 'pos', 'alt', 'ref')\n",
" if(length(bim_files) > 1 ){\n",
" for(i in 2:length(bim_files)){\n",
" message('Aligning ', bim_files[i])\n",
" tmp_bim <- fread(bim_files[i], header = FALSE)\n",
" colnames(tmp_bim) <- c('chr', 'id', 'Mics', 'pos', 'alt', 'ref')\n",
" rosmap_bim <- merge_genotype_data(rosmap_bim, tmp_bim)\n",
" }\n",
" } else {\n",
" rosmap_bim <- rosmap_bim[, .(chr, pos, alt, ref)]\n",
" }\n",
" rosmap_bim <- rosmap_bim %>% arrange(pos)\n",
" fwrite(rosmap_bim, ${_output[0]:nr}, sep = '\\t', col.names = FALSE)\n",
" system(paste(\"bgzip \", ${_output[0]:nr}))\n",
" system(paste(\"tabix -s 1 -b 2 -e 2 \", ${_output[0]:r}))"
"bash: expand = '${ }', stdout = f\"{_output[0]:nn}.stdout\", stderr = f\"{_output[0]:nn}.stderr\", container = container\n",
" Rscript ${modular_script_dir}/misc/genotype_alignment.R \\\n",
" --output \"${_output[0]}\" ${_input:r}\n"
]
}
],
Expand Down
Loading
Loading