Skip to content

fixes for pecotmr #1418

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
gaow merged 4 commits into from
Sep 2, 2026
Merged

fixes for pecotmr #1418

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
60bce4a
fixes for pecotmr
danielnachun Aug 31, 2026
104d5cb
debug OSTYPE
danielnachun Aug 31, 2026
4ddec57
more fixes
danielnachun Sep 1, 2026
d9066a3
remove debug
danielnachun Sep 1, 2026
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
95 changes: 5 additions & 90 deletions code/SoS/data_preprocessing/phenotype/gene_annotation.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
"\n",
"Coordinates come from a collapsed gene-model GTF, following the [GTEx pipeline](https://github.com/broadinstitute/gtex-pipeline/blob/master/qtl/src/eqtl_prepare_expression.py) convention, so the annotation used here matches the one used to build the GTF in the first place. Feature IDs that the GTF does not contain are dropped rather than guessed at, which is why the row count of the output can be lower than that of the input.\n",
"\n",
"Which route applies depends on how the features are named. Gene-level expression matrices match the GTF directly. Protein matrices carry `gene_id|UniProt` composite IDs and need an auxiliary mapping to reach a gene. Leafcutter intron clusters are not genes at all: they are genomic intervals that must first be mapped onto the genes they overlap, which also produces the phenotype-group file that grouped multiple-testing correction later consumes. Psichomics PSI output is already grouped by gene symbol, so it only needs that symbol translated to an ENSEMBL ID.\n",
"Which route applies depends on how the features are named. Gene-level expression matrices match the GTF directly. Protein matrices carry `gene_id|UniProt` composite IDs and need an auxiliary mapping to reach a gene. Leafcutter intron clusters are not genes at all: they are genomic intervals that must first be mapped onto the genes they overlap, which also produces the phenotype-group file that grouped multiple-testing correction later consumes.\n",
"\n",
"An older `biomaRt` route is retained for reference. It places a single-base TSS interval rather than the full gene body and is obsolete; new analyses should use the GTF-based route.\n",
"\n",
Expand Down Expand Up @@ -79,7 +79,7 @@
"- `<phenoFile stem>.region_list.txt` -- one row per retained feature giving its coordinates, ID and the path of the `bed.gz` it came from; QTL modules read this to enumerate cis regions.\n",
"- `<phenoFile stem>.gene_list.tsv` -- the feature IDs that were successfully matched to the GTF.\n",
"- `*.leafcutter.clusters_to_genes.txt` and `*.exon_list` -- from `map_leafcutter_cluster_to_gene`: the cluster-to-gene assignment and the exon table it was derived from.\n",
"- `*.phenotype_group.txt` -- from the leafcutter and psichomics isoform workflows: the isoform-to-gene grouping used by grouped multiple-testing correction.\n",
"- `*.phenotype_group.txt` -- from the leafcutter isoform workflow: the isoform-to-gene grouping used by grouped multiple-testing correction.\n",
"\n",
"`tests/fixtures/gene_annotation/expected/protocol_example.rnaseq.bed.bed.gz`:\n",
"\n",
Expand Down Expand Up @@ -313,47 +313,6 @@
" --map-stra site\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"kernel": "SoS"
},
"source": [
"### Psichomics isoforms"
]
},
{
"cell_type": "markdown",
"metadata": {
"kernel": "SoS"
},
"source": [
"For psichomics splicing quantifications, each event ID ends in `_<gene_id>`; the gene ID is matched against the GTF to obtain coordinates. The output is a coordinate-annotated phenotype BED plus a phenotype-group file."
]
},
{
"cell_type": "markdown",
"metadata": {
"kernel": "SoS"
},
"source": [
"**Timing**: ~2 min (on toy dataset)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"kernel": "Bash"
},
"outputs": [],
"source": [
"sos run pipeline/gene_annotation.ipynb annotate_psichomics_isoforms \\\n",
" --cwd output/gene_annotation \\\n",
" --phenoFile tests/fixtures/gene_annotation/protocol_example.psichomics.phenotype.tsv \\\n",
" --coordinate-annotation <path/to/Homo_sapiens.GRCh38.103.chr.gtf>\n"
]
},
{
"cell_type": "markdown",
"metadata": {
Expand Down Expand Up @@ -476,7 +435,6 @@
" annotate_coord_biomart\n",
" map_leafcutter_cluster_to_gene\n",
" annotate_leafcutter_isoforms\n",
" annotate_psichomics_isoforms\n",
"\n",
"Global Workflow Options:\n",
" --modular-script-dir code/script (as path)\n",
Expand Down Expand Up @@ -523,9 +481,7 @@
" map_leafcutter_cluster_to_gene:\n",
" Workflow Options:\n",
" --intron-count VAL (as path, required)\n",
" Extract the code in case psichromatic needs to be\n",
" processed the same way PheoFile in this step is the\n",
" intron_count file\n",
" PhenoFile in this step is the intron_count file\n",
" --coordinate-annotation VAL (as path, required)\n",
" gtf annotation with exons (the full gene model), used to\n",
" map intron splice sites to genes\n",
Expand All @@ -545,11 +501,6 @@
" --sample-participant-lookup . (as path)\n",
" --coordinate-annotation VAL (as path, required)\n",
" gtf annotation with exons (the full gene model)\n",
" annotate_psichomics_isoforms:\n",
" Workflow Options:\n",
" --sample-participant-lookup . (as path)\n",
" --coordinate-annotation VAL (as path, required)\n",
" gtf annotation with exons (the full gene model)\n",
"```"
]
},
Expand Down Expand Up @@ -704,8 +655,7 @@
"outputs": [],
"source": [
"[map_leafcutter_cluster_to_gene]\n",
"## Extract the code in case psichromatic needs to be processed the same way\n",
"## PheoFile in this step is the intron_count file\n",
"## PhenoFile in this step is the intron_count file\n",
"parameter: intron_count = path\n",
"# gtf annotation with exons (the full gene model), used to map intron splice sites to genes\n",
"parameter: coordinate_annotation = path\n",
Expand Down Expand Up @@ -754,41 +704,6 @@
" --output-bed \"${_output[0]}\" \\\n",
" --output-phenotype-group \"${_output[1]}\""
]
},
{
"cell_type": "markdown",
"metadata": {
"kernel": "SoS"
},
"source": [
"### Processing of psichomics output\n",
"It occurs that the psichomatic by default grouped the isoforms by gene name, so only thing needs to be done is to extract this information and potentially renamed the gene symbol into ENSG ID"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"kernel": "SoS"
},
"outputs": [],
"source": [
"[annotate_psichomics_isoforms]\n",
"parameter: sample_participant_lookup = path()\n",
"# gtf annotation with exons (the full gene model)\n",
"parameter: coordinate_annotation = path\n",
"input: phenoFile, coordinate_annotation\n",
"output: f'{cwd:a}/{_input[0]:bn}.formated.bed.gz', f'{cwd:a}/{_input[0]:bn}.phenotype_group.txt'\n",
"task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output[0]:bn}' \n",
"bash: expand= \"${ }\", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint = entrypoint\n",
" Rscript ${modular_script_dir}/data_preprocessing/phenotype/gene_annotation.R \\\n",
" --step annotate_psichomics_isoforms \\\n",
" --phenoFile \"${_input[0]}\" \\\n",
" --annotation-gtf \"${_input[1]}\" \\\n",
" --sample-participant-lookup \"${sample_participant_lookup}\" \\\n",
" --output-bed \"${_output[0]}\" \\\n",
" --output-phenotype-group \"${_output[1]}\""
]
}
],
"metadata": {
Expand Down Expand Up @@ -827,4 +742,4 @@
},
"nbformat": 4,
"nbformat_minor": 4
}
}
73 changes: 22 additions & 51 deletions code/SoS/data_preprocessing/phenotype_preprocessing.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@
"\n",
"This mini-protocol walks through how molecular-phenotype matrices are prepared for downstream xQTL analysis. Each step calls a workflow from [`phenotype_imputation.ipynb`](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_imputation.html), [`gene_annotation.ipynb`](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/gene_annotation.html), or [`phenotype_formatting.ipynb`](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html).\n",
"\n",
"The commands are selectable routes rather than one mandatory chain. Impute only matrices containing missing values. Choose the annotation workflow matching gene, protein, LeafCutter, or psichomics identifiers, then choose the formatting workflow required by the downstream analysis unit. GCT sample extraction and BAM subsetting are independent utilities.\n"
"The commands are selectable routes rather than one mandatory chain. Impute only matrices containing missing values. Choose the annotation workflow matching gene, protein, or LeafCutter identifiers, then choose the formatting workflow required by the downstream analysis unit. GCT sample extraction and BAM subsetting are independent utilities.\n"
],
"id": "416d15c6"
},
Expand All @@ -48,21 +48,20 @@
"source": [
"## Steps\n",
"\n",
"Choose a route before running commands; the 12 commands are not one mandatory chain.\n",
"Choose a route before running commands; the 11 commands are not one mandatory chain.\n",
"\n",
"| **Analysis goal** | **Commands to run, in order** | **Inputs** |\n",
"| --- | --- | --- |\n",
"| Gene-expression or protein phenotype by chromosome | 2 → 7 | `tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz`; `input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf` |\n",
"| Phenotype with missing values by chromosome | 1 → 2 → 7 | `tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz`; gene-coordinate GTF above |\n",
"| Gene-expression or protein phenotype by chromosome | 2 → 6 | `tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz`; `input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf` |\n",
"| Phenotype with missing values by chromosome | 1 → 2 → 6 | `tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz`; gene-coordinate GTF above |\n",
"| Retrieve gene coordinates from Ensembl BioMart | 3 | `tests/fixtures/gene_annotation/protocol_example.rnaseq.gene_ID.tsv`; internet access |\n",
"| LeafCutter clusters mapped to genes | 4 | `tests/fixtures/gene_annotation/protocol_example.leafcutter.phenotype.bed.gz`; `tests/fixtures/gene_annotation/protocol_example.leafcutter.intron_count.tsv`; `input/reference_data/Homo_sapiens.GRCh38.103.chr.gtf` |\n",
"| LeafCutter isoforms annotated for QTL analysis | 5 | `tests/fixtures/gene_annotation/protocol_example.leafcutter.phenotype.bed.gz`; `tests/fixtures/gene_annotation/protocol_example.leafcutter.intron_count.tsv`; `input/reference_data/Homo_sapiens.GRCh38.103.chr.gtf` |\n",
"| psichomics isoforms annotated for QTL analysis | 6 | `tests/fixtures/gene_annotation/protocol_example.psichomics.phenotype.tsv`; `input/reference_data/Homo_sapiens.GRCh38.103.chr.gtf` |\n",
"| Partition a GCT matrix by chromosome | 8 | `input/rnaseq/protocol_example.rnaseq.gene_tpm.gct.gz` |\n",
"| Partition a BED phenotype by predefined regions | 2 → 9 | `tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz`; gene-coordinate GTF above; `input/reference_data/TAD/protocol_example_protein.enhanced_cis_chr22.bed` |\n",
"| Define TAD-based phenotype regions | 2 → 10 | `tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz`; gene-coordinate GTF above; `tests/fixtures/generalized_TADB/expected/TADB_enhanced_cis.bed` |\n",
"| Extract selected samples from a GCT matrix | 11 | `input/rnaseq/protocol_example.rnaseq.gene_tpm.gct.gz`; `tests/fixtures/phenotype_formatting/keep_samples.txt` |\n",
"| Subset BAM files to selected genomic regions | 12 | `input/rnaseq/bam_file_list.txt`; referenced BAM files; no example BAM is currently bundled |\n",
"| Partition a GCT matrix by chromosome | 7 | `input/rnaseq/protocol_example.rnaseq.gene_tpm.gct.gz` |\n",
"| Partition a BED phenotype by predefined regions | 2 → 8 | `tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz`; gene-coordinate GTF above; `input/reference_data/TAD/protocol_example_protein.enhanced_cis_chr22.bed` |\n",
"| Define TAD-based phenotype regions | 2 → 9 | `tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz`; gene-coordinate GTF above; `tests/fixtures/generalized_TADB/expected/TADB_enhanced_cis.bed` |\n",
"| Extract selected samples from a GCT matrix | 10 | `input/rnaseq/protocol_example.rnaseq.gene_tpm.gct.gz`; `tests/fixtures/phenotype_formatting/keep_samples.txt` |\n",
"| Subset BAM files to selected genomic regions | 11 | `input/rnaseq/bam_file_list.txt`; referenced BAM files; no example BAM is currently bundled |\n",
"\n",
"Run only the row matching the intended analysis goal, following its commands in numerical order. Other imputation methods are available through `phenotype_imputation.ipynb` and its Command Interface."
],
Expand Down Expand Up @@ -214,34 +213,7 @@
"kernel": "SoS"
},
"source": [
"### 6. [Annotate psichomics isoforms](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/gene_annotation.html)\n",
"\n",
"**What it does:** Add genomic and gene annotations to psichomics-derived splicing phenotypes.\n"
],
"id": "acfd2814"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"kernel": "Bash"
},
"outputs": [],
"source": [
"sos run pipeline/gene_annotation.ipynb annotate_psichomics_isoforms \\\n",
" --cwd output/gene_annotation \\\n",
" --phenoFile tests/fixtures/gene_annotation/protocol_example.psichomics.phenotype.tsv \\\n",
" --coordinate-annotation <path/to/Homo_sapiens.GRCh38.103.chr.gtf>\n"
],
"id": "e6d20f6d"
},
{
"cell_type": "markdown",
"metadata": {
"kernel": "SoS"
},
"source": [
"### 7. [Partition a BED phenotype by chromosome](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"### 6. [Partition a BED phenotype by chromosome](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"\n",
"**What it does:** Split a coordinate-annotated BED phenotype into chromosome-specific files.\n"
],
Expand Down Expand Up @@ -269,7 +241,7 @@
"kernel": "SoS"
},
"source": [
"### 8. [Partition a GCT phenotype by chromosome](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"### 7. [Partition a GCT phenotype by chromosome](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"\n",
"**What it does:** Split a coordinate-aware GCT matrix into chromosome-specific GCT files.\n"
],
Expand All @@ -296,7 +268,7 @@
"kernel": "SoS"
},
"source": [
"### 9. [Partition a phenotype by predefined regions](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"### 8. [Partition a phenotype by predefined regions](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"\n",
"**What it does:** Extract phenotype features falling within each region in a supplied region list.\n"
],
Expand All @@ -323,7 +295,7 @@
"kernel": "SoS"
},
"source": [
"### 10. [Define TAD-based phenotype regions](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"### 9. [Define TAD-based phenotype regions](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"\n",
"**What it does:** Assign phenotype features to TAD windows and generate a region list for downstream analysis.\n"
],
Expand Down Expand Up @@ -351,7 +323,7 @@
"kernel": "SoS"
},
"source": [
"### 11. [Extract selected samples from a GCT matrix](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"### 10. [Extract selected samples from a GCT matrix](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"\n",
"**What it does:** Retain only samples listed in a supplied keep file.\n"
],
Expand All @@ -378,7 +350,7 @@
"kernel": "SoS"
},
"source": [
"### 12. [Subset BAM files by genomic region](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"### 11. [Subset BAM files by genomic region](https://statfungen.github.io/xqtl-protocol/code/data_preprocessing/phenotype/phenotype_formatting.html)\n",
"\n",
"**What it does:** Extract selected chromosomes or regions from every BAM listed in the input manifest.\n"
],
Expand Down Expand Up @@ -414,13 +386,12 @@
"| Step 3 | `output/gene_annotation/protocol_example.rnaseq.gene_ID.bed.gz` |\n",
"| Step 4 | `output/gene_annotation/protocol_example.leafcutter.phenotype.bed.gz.exon_list`; `output/gene_annotation/protocol_example.leafcutter.phenotype.bed.gz.leafcutter.clusters_to_genes.txt` |\n",
"| Step 5 | `tests/fixtures/gene_annotation/expected/protocol_example.leafcutter.phenotype.bed.formated.bed.gz`; `output/gene_annotation/protocol_example.leafcutter.phenotype.phenotype_group.txt` |\n",
"| Step 6 | `tests/fixtures/gene_annotation/expected/protocol_example.psichomics.phenotype.formated.bed.gz`; `tests/fixtures/gene_annotation/expected/protocol_example.psichomics.phenotype.phenotype_group.txt` |\n",
"| Step 7 | `output/phenotype_uf/protocol_example.genotype.chr22.bed.gz`; `tests/fixtures/phenotype_formatting/expected/protocol_example.phenotype_by_chrom_files.txt`; its region list |\n",
"| Step 8 | `output/phenotype_gct/protocol_example.rnaseq.gene_tpm.chr21.gct`; `output/phenotype_gct/protocol_example.rnaseq.gene_tpm.chr22.gct` |\n",
"| Step 9 | `output/phenotype_by_region/protocol_example_protein.enhanced_cis_chr22_phenotype_by_region/*.bed.gz`; `output/phenotype_by_region/*.phenotype_by_region_files.txt` |\n",
"| Step 10 | `output/phenotype_by_region/*_pheno_per_region.region_list`; TAD-annotated region list under `output/phenotype_by_region/` |\n",
"| Step 11 | `output/phenotype_gct/protocol_example.rnaseq.gene_tpm.sample_matched.gct.gz` |\n",
"| Step 12 | `output/bam_subset/*.subsetted.bam` |"
"| Step 6 | `output/phenotype_uf/protocol_example.genotype.chr22.bed.gz`; `tests/fixtures/phenotype_formatting/expected/protocol_example.phenotype_by_chrom_files.txt`; its region list |\n",
"| Step 7 | `output/phenotype_gct/protocol_example.rnaseq.gene_tpm.chr21.gct`; `output/phenotype_gct/protocol_example.rnaseq.gene_tpm.chr22.gct` |\n",
"| Step 8 | `output/phenotype_by_region/protocol_example_protein.enhanced_cis_chr22_phenotype_by_region/*.bed.gz`; `output/phenotype_by_region/*.phenotype_by_region_files.txt` |\n",
"| Step 9 | `output/phenotype_by_region/*_pheno_per_region.region_list`; TAD-annotated region list under `output/phenotype_by_region/` |\n",
"| Step 10 | `output/phenotype_gct/protocol_example.rnaseq.gene_tpm.sample_matched.gct.gz` |\n",
"| Step 11 | `output/bam_subset/*.subsetted.bam` |"
],
"id": "48f58267"
},
Expand Down Expand Up @@ -501,4 +472,4 @@
},
"nbformat": 4,
"nbformat_minor": 5
}
}
Loading
Loading