Prepare molecular-phenotype matrices for xQTL analysis by imputing missing values when needed, adding genomic annotations, and formatting data for chromosome-, region-, TAD-, or sample-specific analyses.
Miniprotocol Timing¶
This is the total duration for the standard toy-data route; module-specific timings appear on their respective pages.
Timing: <12 min (on the toy dataset)
Overview¶
This mini-protocol walks through how molecular-phenotype matrices are prepared for downstream xQTL analysis. Each step calls a workflow from phenotype_imputation.ipynb, gene_annotation.ipynb, or phenotype_formatting.ipynb.
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.
Steps¶
Choose a route before running commands; the 11 commands are not one mandatory chain.
| Analysis goal | Commands to run, in order | Inputs |
|---|---|---|
| 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 |
| Phenotype with missing values by chromosome | 1 → 2 → 6 | tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz; gene-coordinate GTF above |
| Retrieve gene coordinates from Ensembl BioMart | 3 | tests/fixtures/gene_annotation/protocol_example.rnaseq.gene_ID.tsv; internet access |
| 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 |
| 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 |
| Partition a GCT matrix by chromosome | 7 | input/rnaseq/protocol_example.rnaseq.gene_tpm.gct.gz |
| 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 |
| 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 |
| Extract selected samples from a GCT matrix | 10 | input/rnaseq/protocol_example.rnaseq.gene_tpm.gct.gz; tests/fixtures/phenotype_formatting/keep_samples.txt |
| Subset BAM files to selected genomic regions | 11 | input/rnaseq/bam_file_list.txt; referenced BAM files; no example BAM is currently bundled |
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.
1. Impute missing phenotype values¶
What it does: Use generalized empirical Bayes matrix factorization to complete the example protein matrix.
sos run pipeline/phenotype_imputation.ipynb gEBMF \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf \
--num_factor 30
2. Add genomic coordinates to gene or protein phenotypes¶
What it does: Join phenotype identifiers to a supplied coordinate annotation and write a coordinate-aware BED matrix.
sos run pipeline/gene_annotation.ipynb annotate_coord \
--cwd output/gene_annotation \
--phenoFile tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz \
--coordinate-annotation tests/fixtures/gene_annotation/Homo_sapiens.GRCh38.103.collapse_only.gene.chr22.gtf.gz \
--phenotype-id-column gene_id
3. Retrieve gene coordinates from Ensembl BioMart¶
What it does: Query the selected Ensembl release when a local coordinate annotation is unavailable.
sos run pipeline/gene_annotation.ipynb annotate_coord_biomart \
--cwd output/gene_annotation \
--phenoFile tests/fixtures/gene_annotation/protocol_example.rnaseq.gene_ID.tsv \
--ensembl-version 115
4. Map LeafCutter clusters to genes¶
What it does: Assign LeafCutter clusters to genes using splice-site overlap with the gene annotation.
sos run pipeline/gene_annotation.ipynb map_leafcutter_cluster_to_gene \
--cwd output/gene_annotation \
--phenoFile tests/fixtures/gene_annotation/protocol_example.leafcutter.phenotype.bed.gz \
--intron-count tests/fixtures/gene_annotation/protocol_example.leafcutter.intron_count.tsv \
--coordinate-annotation <path/to/Homo_sapiens.GRCh38.103.chr.gtf> \
--map-stra site
5. Annotate LeafCutter isoforms¶
What it does: Convert LeafCutter intron-cluster phenotypes into annotated isoform features for QTL analysis.
sos run pipeline/gene_annotation.ipynb annotate_leafcutter_isoforms \
--cwd output/gene_annotation \
--phenoFile tests/fixtures/gene_annotation/protocol_example.leafcutter.phenotype.bed.gz \
--intron-count tests/fixtures/gene_annotation/protocol_example.leafcutter.intron_count.tsv \
--coordinate-annotation <path/to/Homo_sapiens.GRCh38.103.chr.gtf> \
--map-stra site
6. Partition a BED phenotype by chromosome¶
What it does: Split a coordinate-annotated BED phenotype into chromosome-specific files.
sos run pipeline/phenotype_formatting.ipynb phenotype_by_chrom \
--cwd output/phenotype_uf \
--phenoFile tests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz \
--name protocol_example \
--chrom chr22
7. Partition a GCT phenotype by chromosome¶
What it does: Split a coordinate-aware GCT matrix into chromosome-specific GCT files.
sos run pipeline/phenotype_formatting.ipynb phenotype_by_chrom_gct \
--cwd output/phenotype_gct \
--phenoFile output/phenotype/phenotype_by_chrom_for_cis/protocol_example.rnaseq.gene_tpm.gct.gz \
--chrom chr21 chr22
8. Partition a phenotype by predefined regions¶
What it does: Extract phenotype features falling within each region in a supplied region list.
sos run pipeline/phenotype_formatting.ipynb phenotype_by_region \
--cwd output/phenotype_by_region \
--phenoFile tests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz \
--region-list output/phenotype/phenotype_by_chrom_for_cis/protocol_example_protein.enhanced_cis_chr22.bed
9. Define TAD-based phenotype regions¶
What it does: Assign phenotype features to TAD windows and generate a region list for downstream analysis.
sos run pipeline/phenotype_formatting.ipynb phenotype_annotate_by_tad \
--cwd output/phenotype_by_region \
--phenoFile tests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz \
--TAD-list tests/fixtures/generalized_TADB/expected/TADB_enhanced_cis.bed \
--phenotype-per-tad 2
10. Extract selected samples from a GCT matrix¶
What it does: Retain only samples listed in a supplied keep file.
sos run pipeline/phenotype_formatting.ipynb gct_extract_samples \
--cwd output/phenotype_gct \
--phenoFile output/phenotype/phenotype_by_chrom_for_cis/protocol_example.rnaseq.gene_tpm.gct.gz \
--keep-samples tests/fixtures/phenotype_formatting/keep_samples.txt
11. Subset BAM files by genomic region¶
What it does: Extract selected chromosomes or regions from every BAM listed in the input manifest.
sos run pipeline/phenotype_formatting.ipynb bam_subsetting \
--cwd output/bam_subset \
--phenoFile output/phenotype/phenotype_by_chrom_for_cis/bam_file_list.txt \
--region chr21 chr22
Output¶
| Route or step | Products and relative paths |
|---|---|
| Step 1 | output/phenotype_imputation_uf/protocol_example.protein.missing.bed.imputed.bed.gz |
| Step 2 | tests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz; tests/fixtures/gene_annotation/expected/protocol_example.rnaseq.bed.region_list.txt |
| Step 3 | output/gene_annotation/protocol_example.rnaseq.gene_ID.bed.gz |
| 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 |
| 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 |
| 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 |
| Step 7 | output/phenotype_gct/protocol_example.rnaseq.gene_tpm.chr21.gct; output/phenotype_gct/protocol_example.rnaseq.gene_tpm.chr22.gct |
| 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 |
| Step 9 | output/phenotype_by_region/*_pheno_per_region.region_list; TAD-annotated region list under output/phenotype_by_region/ |
| Step 10 | output/phenotype_gct/protocol_example.rnaseq.gene_tpm.sample_matched.gct.gz |
| Step 11 | output/bam_subset/*.subsetted.bam |
Anticipated Results¶
The selected route produces a molecular-phenotype matrix with the genomic annotation and layout required by its downstream xQTL analysis. Standard gene-expression routes typically end with chromosome-specific BED files; splicing routes end with gene- or isoform-annotated phenotypes; region and TAD routes end with per-region files and their manifest.
Continue with covariate preprocessing and then the appropriate association-testing workflow.
Command Interface¶
List the workflows and parameters available in each module used by this mini-protocol.
sos run pipeline/phenotype_imputation.ipynb -h
sos run pipeline/gene_annotation.ipynb -h
sos run pipeline/phenotype_formatting.ipynb -h