APA Calling

APA Calling#

Quantifies alternative polyadenylation from RNA-seq coverage with DaPars2, producing per-sample PDUI values.

Overview#

Most genes have more than one polyadenylation site, and which one is used shifts the length of the 3’ UTR - changing which regulatory elements survive in the transcript. Treating a gene as a single expression value hides that. This step quantifies alternative polyadenylation from RNA-seq coverage with DaPars2, giving a per-sample PDUI (percentage of distal polyA site usage) that can be scanned for QTLs like any other molecular phenotype.

The workflows run in order: UTR_reference builds the annotation DaPars2 needs, bam2tools converts alignments to coverage, APAconfig writes the configuration file and APAmain runs the quantification.

When to run it. After RNA-seq alignment, as an alternative or complement to expression and splicing quantification. The output feeds phenotype preprocessing like any other molecular phenotype matrix.

%preview ../../images/apa_calling.png

Input#

The workflow runs in four stages, and each stage consumes the output of the one before it. Only the first two stages need files you supply yourself.

--cwd output/apa (working directory all outputs are written under. Used by every stage. Defaults to output.)

Step 1, UTR_reference.

--hg-gtf output/apa/chr22.gtf (a gene annotation GTF matching the genome build used for alignment. The step reads transcript records from it and derives the 3’UTR intervals that DaPars2 will test. Built from the transcriptome-level gene feature file prepared by the reference data module.)

Step 2, bam2tools.

--bam-dir (a directory of transcriptome-aligned BAM files, one per sample. Sample names are taken from the BAM file stems, so those stems must match the identifiers you expect downstream.)
--bam-list (an explicit list of BAM paths; use instead of --bam-dir)

Step 3, APAconfig.

--bfile output/apa/wig (the directory of WIG tracks written by bam2tools. The tracks must end in .wig; renaming a .bedgraph is enough. The first line is treated as a header, so add any content there if the file has none. If the first column does not carry the chr prefix, set no_chr_prefix = T.)
--annotation output/apa/chr22_3UTR.bed (the 3’UTR BED written by UTR_reference)
--coverage-threshold 10 (minimum read depth a sample must reach in a region)
--least-pass-coverage-percentage 0.3 (fraction of samples that must clear that depth)

Step 4, APAmain.

--chrlist chr22 (the chromosomes to process)
--chr-prefix true (set to False if contig names in your BAMs lack the chr prefix)
--dapars-path code/SoS/molecular_phenotypes/calling/apa (location of the bundled DaPars2 scripts; change it only if you have moved them)

From the example data, showing the first rows (including the header where the format has one):

output/apa/chr22.gtf

chr22	ensembl	gene	10736171	10736283	.	-	.	...
chr22	ensembl	transcript	10736171	10736283	.	-	.	...
chr22	ensembl	exon	10736171	10736283	.	-	.	...

output/apa/chr22_3UTR.bed

chr22	32037111	32037995	ENST00000658411|LINC02558|chr22|+	0	+
chr22	46352525	46352841	ENST00000463785|TRMU|chr22|+	0	+
chr22	34175390	34175780	ENST00000668901|LINC01643|chr22|+	0	+

Output#

output/apa/chr22_3UTR.bed (UTR_reference) (the extracted 3’UTR intervals; the name field packs four values separated by |: transcript ID, gene symbol, chromosome and strand)

chr22	32037111	32037995	ENST00000658411|LINC02558|chr22|+	0	+
chr22	46352525	46352841	ENST00000463785|TRMU|chr22|+	0	+

output/apa/wig/SAMPLE_001.wig (bam2tools) (per-base coverage track, one per BAM)

track type=bedGraph
chr22	0	10527919	0
chr22	10527919	10527927	1

output/apa/wig/SAMPLE_001.flagstat (bam2tools) (total mapped read count used to normalise that sample. The step declares this output as <sample>.depth; the example data carries .flagstat, so the two disagree and the code is the authority.)

202561 + 0 in total (QC-passed reads + QC-failed reads)
200000 + 0 primary

output/apa/sample_mapping_files.txt (APAconfig) (maps each sample’s WIG track to its read count)

output/apa/wig/SAMPLE_001.wig	202561
output/apa/wig/SAMPLE_002.wig	201999

output/apa/sample_configuration_file.txt (APAconfig) (the DaPars2 run configuration, carrying the coverage thresholds)

Annotated_3UTR=output/apa/chr22.gencode_3UTR.bed
Aligned_Wig_files=output/apa/wig/SAMPLE_001.wig,output/apa/wig/SAMPLE_002.wig
Output_directory=output/apa/apa 
Output_result_file=Dapars_result
Least_pass_coverage_percentage=0.3
Coverage_threshold=10
Num_Threads=1
sequencing_depth_file=output/apa/sample_mapping_files.txt

apa_<chr>/Dapars_result_result_temp.<chr>.txt (APAmain) (the PDUI (Percentage of Distal polyA site Usage Index) estimates, one row per 3’UTR region and one PDUI column per sample. A PDUI near 1 means the distal polyA site dominates; lower values indicate 3’UTR shortening.)

Minimal Working Example#

Quantify alternative polyadenylation with DaPars2#

The four workflows chain: each consumes the output of the one before it, so run them in order.

Step 1. Build the 3’UTR reference#

UTR_reference converts the GTF to BED12 with gtf2bed12.py (from Li Lab’s Exon_Intron_Extractor), then extracts the 3’UTR intervals from the whole-genome BED with DaPars_Extract_Anno.py.

Timing: TBD (on toy dataset)

sos run pipeline/apa_calling.ipynb UTR_reference \
    --cwd output/apa \
    --hg-gtf output/apa/chr22.gtf

Step 2. Convert the BAM files to coverage#

bam2tools turns each transcriptome-aligned BAM into a per-base coverage track and a read-depth file, using bedtools or rsem-bam2wig for RSEM-based alignment.

Timing: TBD (on toy dataset)

sos run pipeline/apa_calling.ipynb bam2tools \
    --cwd output/apa \
    --bam-dir output/rnaseq/bam

Step 3. Write the DaPars2 configuration#

APAconfig reads each WIG track line by line to total the read coverage across chromosomes, then writes the sample mapping file and the DaPars2 configuration file carrying the coverage thresholds.

Timing: TBD (on toy dataset)

sos run pipeline/apa_calling.ipynb APAconfig \
    --cwd output/apa \
    --bfile output/apa/wig \
    --annotation output/apa/chr22_3UTR.bed

Step 4. Calculate the PDUI matrix#

APAmain runs Dapars2_Multi_Sample.py, which uses least squares to estimate long-isoform usage per chromosome. The bundled copy is modified from source to handle a formatting discrepancy in the WIG files. Missing values in the result are imputed by KNN with the R impute package.

Timing: TBD (on toy dataset)

sos run pipeline/apa_calling.ipynb APAmain \
    --cwd output/apa \
    --chrlist chr22 \
    --chr-prefix true \
    --dapars-path code/SoS/molecular_phenotypes/calling/apa

Command Interface#

sos run pipeline/apa_calling.ipynb -h

usage: sos run pipeline/apa_calling.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:
  UTR_reference
  bam2tools
  APAconfig
  APAmain

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --walltime 400h
  --mem 200G
  --ncore 16 (as int)
  --cwd output (as path)
                        the output directory for generated files
  --numThreads 8 (as int)
                        Number of threads
  --job-size 1 (as int)
  --container ''

Sections
  UTR_reference:
    Workflow Options:
      --hg-gtf VAL (as path, required)
                        gtf file
  bam2tools:
    Workflow Options:
      --bam-list  paths

                        Transcriptome/coordinate-sorted BAM files to compute
                        per-base coverage for
      --bam-dir . (as path)
                        Retained for CLI back-compatibility (unused by the R
                        port, which takes --bam-list)
  APAconfig:
    Workflow Options:
      --bfile VAL (as path, required)
      --annotation VAL (as path, required)
      --least-pass-coverage-percentage 0.3 (as float)
      --coverage-threshold 10 (as int)
  APAmain:
    Workflow Options:
      --chrlist VAL VAL ... (as type, required)
      --[no-]chr-prefix (default to True)
                        Whether the input wig/BAM chromosomes already carry the
                        'chr' prefix (True, default, matching the chr-prefixed
                        reference data). Set False for bare-named BAMs/wigs
                        (e.g. '22'): dapars.R then prepends 'chr' to the wig
                        chroms to match the annotation.
      --dapars-path code/SoS/molecular_phenotypes/calling/apa (as path)
                        Retained for CLI back-compatibility (unused by the R
                        port of DaPars2)

Workflow implementation#

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
parameter: walltime = '400h'
parameter: mem = '200G'
parameter: ncore = 16
# the output directory for generated files
parameter: cwd = path("output")
# Number of threads
parameter: numThreads = 8
parameter: job_size = 1
parameter: container = ''

UTR_reference#

Generate 3’UTR regions based on GTF

[UTR_reference]
# gtf file
parameter: hg_gtf = path
input: hg_gtf
output: f'{cwd}/{_input:bn}.bed',
        f'{cwd}/{_input:bn}.transcript_to_geneName.txt',
        f'{cwd}/{_input:bn}_3UTR.bed'
bash: expand = '${ }', container = container
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/dapars.R \
        --step gtf2bed12 --gtf ${_input} --out-bed ${_output[0]} --out-map ${_output[1]}
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/dapars.R \
        --step extract_anno --bed ${_output[0]} --symbol ${_output[1]} --output ${_output[2]}

bam2tools#

Generate WIG coverage and flagstat files from BAM files

[bam2tools]
# Transcriptome/coordinate-sorted BAM files to compute per-base coverage for
parameter: bam_list = paths
# Retained for CLI back-compatibility (unused by the R port, which takes --bam-list)
parameter: bam_dir = path('.')
input: bam_list, group_by = 1
output: f'{cwd}/{_input:bn}.wig', f'{cwd}/{_input:bn}.depth'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand = "${ }", container = container
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/dapars.R \
        --step coverage --bam ${_input} --output-wig ${_output[0]} --output-depth ${_output[1]}

APAconfig#

Generate the DaPars2 configuration and sample mapping files

[APAconfig]
parameter: bfile = path
parameter: annotation = path
parameter: least_pass_coverage_percentage = 0.3
parameter: coverage_threshold = 10
output: [f'{cwd}/sample_mapping_files.txt', f'{cwd}/sample_configuration_file.txt']
task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = ncore
bash: expand = "${ }", container = container
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/dapars.R \
        --step apa_config --wig-dir ${bfile} --annotation ${annotation} \
        --apa-cwd ${cwd} --coverage-threshold ${coverage_threshold} --num-threads ${numThreads} \
        --config ${_output[1]} --output-mapping ${_output[0]}

APAmain#

Run DaPars2 main to calculate PDUIs

[APAmain]
parameter: chrlist = list
# Whether the input wig/BAM chromosomes already carry the 'chr' prefix (True, default,
# matching the chr-prefixed reference data). Set False for bare-named BAMs/wigs (e.g. '22'):
# dapars.R then prepends 'chr' to the wig chroms to match the annotation.
parameter: chr_prefix = True
# Retained for CLI back-compatibility (unused by the R port of DaPars2)
parameter: dapars_path = path('code/SoS/molecular_phenotypes/calling/apa')
input: for_each = 'chrlist'
output: [f'{cwd}/apa_{x}/Dapars_result_result_temp.{x}.txt' for x in chrlist], group_by = 1
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = ncore
bash: expand = "${ }", container = container
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/dapars.R \
        --step apa_main --config ${cwd}/sample_configuration_file.txt --chr ${_chrlist} \
        --chr-prefix ${chr_prefix}

APA Calling

Contents

APA Calling#

Overview#

Input#

Output#

Minimal Working Example#

Quantify alternative polyadenylation with DaPars2#

Step 1. Build the 3’UTR reference#

Step 2. Convert the BAM files to coverage#

Step 3. Write the DaPars2 configuration#

Step 4. Calculate the PDUI matrix#

Command Interface#

Workflow implementation#

UTR_reference#

bam2tools#

APAconfig#

APAmain#