Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Phenotype Data Formatting

Splits a molecular phenotype matrix by chromosome or by region so the association scan can be parallelised.

Overview

Association testing and fine-mapping run per chromosome or per region, so the phenotype matrix has to be partitioned the same way before that work can be parallelised. These workflows split a phenotype BED by chromosome or by region, annotate features against TAD boundaries when region-based analysis is wanted, and accept GCT-format input as well as BED. Two further workflows trim inputs rather than split them: one subsets BAM files by coordinate, the other drops samples from a GCT matrix.

The pipeline’s author has flagged it as needing improvement, so treat the interface as unstable and re-check -h before relying on any option.

When to run it. After phenotype QC, normalisation and imputation, and before covariate preprocessing and the association scan.

Input

  • --phenoFile: the molecular phenotype data, required by every workflow. A tabix-indexed UCSC bed.gz for the BED workflows, a gct.gz for the GCT ones, and a list of BAM paths for bam_subsetting. Example tests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz:

    #chr   start     end       ID               SAMPLE_001         SAMPLE_002
chr22  10939387  10961337  ENSG00000283047  37.0768403309204   0.0

  • --chrom: the chromosomes to extract, taking one or more values. Used by phenotype_by_chrom and phenotype_by_chrom_gct.

  • --region-list: a 4-column index of the regions to analyse, giving chr, start, end and region name. Used by phenotype_by_region. Example input/reference_data/TAD/protocol_example_protein.enhanced_cis_chr22.bed:

    #chr   start  end       gene_id
chr22  0      18960000  ENSG00000131100_P36543

  • --TAD-list: the TAD boundary file for phenotype_annotate_by_tad, in the same 4-column shape. Example tests/fixtures/generalized_TADB/expected/TADB_enhanced_cis.bed.

  • --phenotype-per-tad: the minimum number of features a TAD must contain to be kept, 2 by default.

  • --region: coordinates passed through to samtools view by bam_subsetting, for example --region chr21 chr22.

  • --keep-samples: for gct_extract_samples, a file listing the sample IDs to retain, one per line. Example tests/fixtures/phenotype_formatting/keep_samples.txt:

    SAMPLE_001
SAMPLE_002
SAMPLE_003

  • --name: the stem of the output files, derived from --phenoFile by default.

  • --cwd: the directory outputs are written to.

Output

  • {cwd}/{name}.{chrom}.bed.gz - one phenotype BED per requested chromosome, with a tabix index, from phenotype_by_chrom.

  • {cwd}/{name}.phenotype_by_chrom_files.txt and {cwd}/{name}.phenotype_by_chrom_files.region_list.txt - the manifest of those per-chromosome files and the matching region list, both written by the collecting substep. Example output/phenotype_protein/protocol_example_protein.phenotype_by_chrom_files.region_list.txt:

    #chr   start     end       ID                      path
chr22  17592135  17628748  ENSG00000131100_P36543  output/phenotype_protein/protocol_example_protein.chr22.bed.gz

  • {cwd}/{region_list}_phenotype_by_region/{name}.{region}.bed.gz and {cwd}/{name}.phenotype_by_region_files.txt - the per-region equivalents from phenotype_by_region.

  • {cwd}/{phenoFile}.{TAD_list}.{n}_pheno_per_region.region_list - the TAD-grouped region list from phenotype_annotate_by_tad.

  • {cwd}/{name}.{chrom}.gct - one GCT per chromosome, from phenotype_by_chrom_gct.

  • {cwd}/{bam}.subsetted.bam - the coordinate-restricted BAM from bam_subsetting.

  • {phenoFile}.sample_matched.gct.gz - the sample-filtered GCT from gct_extract_samples, written beside its input rather than under --cwd.

Minimal Working Example

Partition a phenotype BED by chromosome

phenotype_by_chrom runs two substeps in sequence: _1 pulls the records for each requested chromosome out of the tabix-indexed BED and writes a per-chromosome bed.gz with its index, and _2 collects those into the *_files.txt and *.region_list.txt lists. This is the partitioning cis TensorQTL analysis needs. Only chr22 is present in the toy phenotype data.

Timing: TBD (on toy dataset)

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

Partition a GCT matrix by chromosome

phenotype_by_chrom_gct does the same split for RNASeQC-style GCT input rather than BED, writing one .gct per chromosome.

Timing: TBD (on toy dataset)

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

Partition a phenotype BED by region

phenotype_by_region splits on an explicit region list instead of whole chromosomes, giving one bed.gz per region under a directory named after the list. Like phenotype_by_chrom it runs two substeps in sequence: _1 writes the per-region files and _2 collects them into a _files.txt manifest.

Timing: TBD (on toy dataset)

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

Group features into TAD windows

phenotype_annotate_by_tad assigns each feature to a TAD and emits a region list, so that region-based analysis works on TAD windows rather than fixed cis-windows. --phenotype-per-tad sets the minimum features a TAD must hold to be kept.

Timing: TBD (on toy dataset)

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

Extract samples from a GCT matrix

gct_extract_samples keeps only the samples named in --keep-samples, which is how a GCT from RNASeQC is cut down to the cohort actually being analysed.

Timing: TBD (on toy dataset)

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

Subset BAM files by region

bam_subsetting runs samtools view over each BAM in --phenoFile and keeps only the given coordinates. No BAM ships with this repository, so the command below is a template rather than a runnable example.

Timing: TBD (on toy dataset)

sos run pipeline/phenotype_formatting.ipynb bam_subsetting \
    --cwd output/bam_subset \
    --phenoFile <path>/bam_file_list.txt \
    --region chr21 chr22

Command Interface

sos run pipeline/phenotype_formatting.ipynb -h
usage: sos run pipeline/phenotype_formatting.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  phenotype_by_chrom
  phenotype_annotate_by_tad
  phenotype_by_chrom_gct
  phenotype_by_region
  bam_subsetting
  gct_extract_samples

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        Work directory & output directory
  --container ''
                        The filename namefor output data
  --entrypoint ''
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --numThreads 20 (as int)
                        Number of threads
  --phenoFile VAL (as path, required)
                        Path to the input molecular phenotype data.
  --name  str(phenoFile).removesuffix('.gz').removesuffix('.bed')

                        name for the analysis output

Sections
  phenotype_by_chrom_1:
    Workflow Options:
      --chrom VAL VAL ... (as type, required)
                        list of chroms to extract
  phenotype_by_chrom_2:
  phenotype_annotate_by_tad:
    Workflow Options:
      --TAD-list VAL (as path, required)
      --phenotype-per-tad 2 (as int)
  phenotype_by_chrom_gct_1:
    Workflow Options:
      --chrom VAL VAL ... (as type, required)
                        list of chroms to extract
  phenotype_by_region_1:
    Workflow Options:
      --region-list VAL (as path, required)
                        An index text file with 4 columns specifying the chr,
                        start, end and name of regions to analyze
  phenotype_by_region_2:
  bam_subsetting:
    Workflow Options:
      --region VAL VAL ... (as type, required)
                        Input to `samtools view` coordinates, for example,
                        --region chr21 chr22
  gct_extract_samples:  Extract samples from expression data generated by
                        RNASeQC
    Workflow Options:
      --keep-samples VAL (as path, required)

Workflow implementation

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
import os
# Work directory & output directory
parameter: cwd = path("output")
# The filename namefor output data
parameter: container = ''
import re
parameter: entrypoint= ""
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 20
# Path to the input molecular phenotype data.
parameter: phenoFile = path
# name for the analysis output
parameter: name = str(phenoFile).removesuffix('.gz').removesuffix('.bed')
[phenotype_by_chrom_1]
# list of chroms to extract 
parameter: chrom = list
chrom = list(set(chrom))
# Path to the input molecular phenotype data.
input: phenoFile, for_each = "chrom"
output: f'{cwd}/{name}.{_chrom}.bed.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output:bn}'  
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R --step phenotype_by_chrom_1 \
        --cwd "${cwd}" \
        --phenoFile "${_input}" \
        --name "${name}" \
        --chrom "${_chrom}" \
        --output "${_output}" \
        --numThreads ${numThreads}
[phenotype_by_chrom_2]
# Path to the input molecular phenotype data.
input: group_by = "all"
output: f'{cwd}/{name}.{step_name[:-2]}_files.txt',f'{cwd}/{name}.{step_name[:-2]}_files.region_list.txt'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R --step phenotype_by_chrom_2 \
        --cwd "${cwd}" \
        --phenoFile "${phenoFile}" \
        --inputs ${_input} \
        --output-files "${_output[0]}" \
        --output-region-list "${_output[1]}"
[phenotype_annotate_by_tad]
parameter: TAD_list = path
parameter: phenotype_per_tad = 2 # This is the minimum number of epigenomics marker for a tadb to be considered having a functions.
input: phenoFile,TAD_list
output: f'{cwd}/{_input[0]:b}.{_input[1]:b}.{phenotype_per_tad}_pheno_per_region.region_list'
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R \
        --step phenotype_annotate_by_tad \
        --cwd "${cwd}" \
        --phenoFile "${_input[0]}" \
        --TAD_list "${_input[1]}" \
        --phenotype-per-tad ${phenotype_per_tad} \
        --output "${_output}" \
        --numThreads ${numThreads}

[phenotype_by_chrom_gct_1]
# list of chroms to extract 
parameter: chrom = list
chrom = list(set(chrom))
# Path to the input molecular phenotype data.
input: phenoFile, for_each = "chrom"
output: f'{cwd:a}/{name}.{_chrom}.gct'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output:bn}'  
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R --step phenotype_by_chrom_gct_1 \
        --cwd "${cwd}" \
        --phenoFile "${_input}" \
        --name "${name}" \
        --chrom "${_chrom}" \
        --output "${_output}" \
        --numThreads ${numThreads}
[phenotype_by_region_1]
# An index text file with 4 columns specifying the chr, start, end and name of regions to analyze
parameter: region_list = path
regions = [x.strip().split() for x in open(region_list).readlines() if x.strip() and not x.strip().startswith('#')]
# Path to the input molecular phenotype data.
input: phenoFile, for_each = "regions"
output: f'{cwd}/{region_list:bn}_phenotype_by_region/{name}.{_regions[3]}.bed.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output:bn}'  
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R --step phenotype_by_region_1 \
        --cwd "${cwd}" \
        --phenoFile "${_input}" \
        --name "${name}" \
        --region "${_regions[0]}" "${_regions[1]}" "${_regions[2]}" "${_regions[3]}" \
        --output "${_output}" \
        --numThreads ${numThreads}
[phenotype_by_region_2]
input: group_by = "all"
output: f'{cwd}/{name}.{step_name[:-2]}_files.txt'
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R --step phenotype_by_region_2 \
        --cwd "${cwd}" \
        --inputs ${_input} \
        --output "${_output}"
[bam_subsetting]
# Input to `samtools view` coordinates, for example, --region chr21 chr22
parameter: region = list
# Path to the input molecular phenotype data.
parameter: phenoFile = paths
input: phenoFile , group_by = 1
output: f'{cwd}/{_input:bn}.subsetted.bam'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container=container, entrypoint=entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R --step bam_subsetting \
        --phenoFile "${_input}" \
        --region ${' '.join(region)} \
        --output "${_output}"
# Extract samples from expression data generated by RNASeQC
[gct_extract_samples]
parameter: keep_samples = path
input: phenoFile
output: f'{_input[0]:nn}.sample_matched.gct.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "${ }", stderr = f'{_output:nn}.stderr', stdout = f'{_output:nn}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R \
        --step gct_extract_samples \
        --cwd "${cwd}" \
        --phenoFile "${_input[0]}" \
        --keep-samples "${keep_samples}" \
        --output "${_output}" \
        --numThreads ${numThreads}
Data Preprocessing
Phenotype data imputation
Data Preprocessing
Covariate data preprocessing