IsoComp: Comparing isoform composition between cohorts using high-quality long-read RNA-seq

Table of contents

1. Contributors
2. Introduction
3. Aim
4. Workflow
5. Running the pipeline
6. Methods

Contributors:

• Yutong Qiu (Carnegie Mellon)
• Chia Sin Liew (University of Nebraska-Lincoln)
• Chase Mateusiak (Washington University)
• Rupesh Kesharwani (Baylor College of Medicine)
• Bida Gu (University of Southern California)
• Muhammad Sohail Raza (Beijing Institute of Genomics, Chinese Academy of Sciences/China National Center for Bioinformation)
• Evan Biederstedt (HMS)
• Umran Yaman (UK Dementia Research Institute, University College London)
• Abdullah Al Nahid (Shahjalal University of Science and Technology)
• Trinh Tat (Houston Methodist Research Institute)
• Sejal Modha (Theolytics Limited)
• Jędrzej Kubica (University of Warsaw)

Cite

Introduction

Transcriptomic profiling has gained traction over the past few decades, but its progress has been hindered by short-read sequencing, particularly in tasks such as inferring alternative splicing, allelic imbalance, and isoform variation due to read length and required assembly.

The potential of long-read sequencing lies in its ability to overcome the inherent limitations of short-reads. Tools like Isoseq3 [link: https://www.pacb.com/products-and-services/applications/rna-sequencing/] offer high-quality, polished, assembled full-length isoforms. This advancement allows us to identify alternatively spliced isoforms and detect gene fusions. Further, with the introduction of HiFi sequencing, the error rates have significantly decreased in third-generation sequencing long reads.

Aim

The aim of this project is to algorithmically characterize the "unique" (differing) isoforms between any number of samples using high-quality assembled isoforms.

Workflow

Running the pipeline

Installation

pip install isocomp==0.3.0

For guidelines run:

isocomp --help

Step 1. Create windows

isocomp create_windows -i sample1.gtf sample2.gtf sample3.gtf -f transcript -o clustered_file.gtf

Step 2. Find unique isoforms across multiple samples

isocomp find_unique_isoforms -a clustered_file.gtf -f fasta_map.csv

File fasta_map.csv:

source,fasta NA24385.filtered,BCM-data-HG002-All2Samples-hg38-Results/NA24385_HG002/MMSQANTI3Filter/NA24385.filtered.fasta NA26105.filtered,BCM-data-HG002-All2Samples-hg38-Results/NA26105_GM26105/MMSQANTI3Filter/NA26105.filtered.fasta

Example output

For each isoform that is unique to at least one sample, we provide information about the read and the similarity between that isoform and the most similar isoform within the same window.

The last column describes the normalized edit distance and the CIGAR string.

win_chr win_start win_end total_isoform isoform_name sample_from sample_compared_to mapped_start isoform_sequence selected_alignments NC_060925.1 255178 288416 4 PB.6.2 HG004 HG002 255173 GGATTATCCGGAGCCAAGGTCCGCTCGGGTGAGTGCCCTCCGCTTTTT 0.02_HG002_PB.6.2_3=6I1=3I1286=11I NC_060925.1 255178 288416 4 PB.6.2 HG004 HG005 255173 GGATTATCCGGAGCCAAGGTCCGCTCGGGTGAGTGCCCTCCGCTTTTTG 0.02_HG002_PB.6.2_3=6I1=3I1286=11

Detailed project overview

https://github.com/collaborativebioinformatics/isocomp/blob/main/FinalPresentationBCMHackathon_12Oct2022.pdf

Methods

Overview of Methods

The core challenge, referred to as the "Isoform set comparison problem," involves identifying distinct isoforms between two sets of samples.

A direct approach to solving this problem is through sequence matching between the sets of isoforms. However, this becomes time-consuming given the size of isoform sets in the human genome (consider a number significantly larger than 10,000, for instance).

We recognize that an all-against-all alignment across complete isoform sets isn't necessary. Instead, the focus is on comparing isoforms aligned to the same genomic regions. Genomic windows containing at least one isoform from any sample are extracted. The isoform sets are then subdivided into smaller subsets based on their origin in these extracted regions.

For each pair of samples under comparison, intersections are made between subsets of isoforms within each genomic window. This process identifies isoforms shared by both samples and isoforms unique to each sample.

For each unique isoform S from sample A, a deeper analysis is conducted on the differences between S and other isoforms from sample B within the same genomic window.

Aligning isoforms to the reference genome

For each individual sample, we initially prefix the sample name to the FASTA sequences in the finalized corrected FASTA output from SQANTI. This step ensures the uniqueness of all sequence names. Subsequently, we employ the minimap2 aligner (v2.24-r1122) to align the renamed FASTA sequences against the human Telomere-to-Telomere genome assembly of the CHM13 cell line (T2T-CHM13v2.0; RefSeq - GCF_009914755.1). The resultant alignment is presented in a SAM file, which is then converted into BAM format and sorted using samtools (v1.15.1; Danecek et al, 2021).

Segmentation of Isoforms into Subsets

Regions from the CHM13v2.0 genome that intersect with at least one isoform from any given sample are extracted. We initially determine the average coverage of isoforms per base using samtools mpileup (version?, Danecek et al, 2021). Subsequently, we identify and extract the 20,042 annotated protein-coding gene regions from the reference genome. To create windows, we merge these regions where overlaps occur. Further refinement is applied by filtering windows to those displaying per-base coverage greater than 0.05, resulting in a final set of 11,936 windows.

Apart from the annotated gene regions, each sample encompasses over 100,000 isoforms (Table 1) that align with intron regions. These isoforms, often considered novel, hold potential relevance to the observed phenotypes. To account for these, we divide the genome into 100-base-pair windows and retain those exhibiting per-base coverage exceeding 0.05.

Following this, the gene-related windows and the 100-base-pair windows are merged to form a comprehensive set of windows aligning with any isoform.

Intersecting subsets of isoforms

For every isoform S within the subset of sample A, we conduct precise string matching against all isoforms within the subset of sample B. If no isoform in sample B, within the same genomic window, precisely matches S, we classify S as unique.

Comparing unique isoforms with other isoforms

For each unique isoform U, we employ the Needleman-Wunch alignment method to compare U with other isoforms within the identical genomic window. The comparison is quantified through the percentage of matched bases in U.

Annotating the differences between unique isoforms and the other sequences

Differences among isoforms are categorized into [TODO] SNPs (<5bp), large-scale variants (>5bp), gene fusions, distinct exon utilization, and entirely novel sequences. These categories build upon those used by SQANTI to annotate disparities between sample isoforms and the reference transcriptome. Notably, we extend SQANTI's categories by incorporating SNPs and large-scale variants.

Iso-Seq analysis

Isoseq3 (v3.2.2) generated HQ (Full-length high quality) transcripts [Table 1] were mapped to GRCh38 (v33 p13) using Minimap2 long-read alignment tools 1. Basic statistics of the alignment for each sample [NA24385 /HG002, NA24143/HG004, and NA24631/HG005] are provided in Table 2. The cDNAcupcake workflow [https://github.com/Magdoll/cDNACupcake] was then executed to collapse redundant isoforms from the BAM file. Low-count isoforms (<10) were filtered out, as well as 5' degraded isoforms that might lack biological significance. Subsequently, SQANTI3 [2] was employed to generate the final corrected fasta [Table 3a] transcripts and GTF [Table 3b] files, along with isoform classification reports. External databases, including the reference data set of transcription start sites (refTSS), a list of polyA motifs, tappAS-annotation, and Genecode hg38 annotation, were utilized. Finally, IsoAnnotLite (v2.7.3) analysis was conducted to annotate the GTF file from SQANTI3.

Differences between isoforms are categorized into [TODO] SNPs (<5bp), large-scale variants (>5bp), gene fusion, different exon usage, and completely novel sequences. These categories build upon those used by SQANTI to annotate disparities between sample isoforms and the reference transcriptome. Note that we extend the categories provided by SQANTI by adding SNPs and large-scale variants.

DEPENDENCIES

python >=3.9

If you're working on ada, you'll need to update the old, crusty version of python to something more modern and exciting.

The easy way (untested, but should work):

Install miniconda and create a conda env with python 3.9

The manual method (source) (tested, works):

``` ssh ... # your username login to ada

mkdir /home/${USER}/.local

use your favorite text editor. no need to be vim

vim /home/${USER}/.bashrc

add the following to the end (or where ever)

export PATH=/home/$USER/.local/bin:$PATH

logout of the current session and log back in

exit ssh ... (your username, etc)

Download a more current version of python

wget https://www.python.org/ftp/python/3.9.15/Python-3.9.15.tgz

unpack

tar xfp Python-3.9.15.tgz

remove the tarball

rm Python-3.9.15.tgz

cd into the Python package dir, configure and make

cd Python-3.9.15/

./configure --prefix=/home/${USER}/.local --exec_prefix=/home/${USER}/.local --enable-optimizations

make # this takes some time

make altinstall

the following should point at a python in your /home/$USER/.local/bin dir

which python3.9

optional, but convenient

ln -s /home/$USER/.local/bin/python3.9 /home/$USER/.local/bin/python

Download the pip installer

curl https://bootstrap.pypa.io/get-pip.py -o get-pip.py

install pip

python3.9 get-pip.py

confirm that pip is where you think it is

which pip # location should be in your .local

at this point, you can do:

pip install poetry

and continue with the development install below

```

Github Codespace for Development

To use codespaces for development purposes, do the following:

1. fork the repo
2. switch to the 'develop' branch
   • NOTE: if you plan to code/add a feature, create a branch from the 'develop' branch. Switch to it, and then continue on with the steps below.
3. click the green 'code' button. But, rather than copying the https or ssh link, click the tab that says "Codespace"
4. click the button that says "create codespace on develop". Go make some tea -- it takes ~5 minutes or so to set up the environment. But, once it is set up, you will have a fully functioning vscode environment with all the dependencies installed. Start running the tests, set some breakpoints, take a look around!

Development

Install poetry and consider setting the configuration such that virtual environments for a given projects are installed in that project directory.

Next, I like working on a fork rather than the actual repository of record. I set my git remotes so that origin points to my fork, and upstream points to the 'upstream' repository.

bash ➜ isocomp git:(develop) ✗ git remote -v origin https://github.com/cmatKhan/isocomp.git (fetch) origin https://github.com/cmatKhan/isocomp.git (push) upstream https://github.com/collaborativebioinformatics/isocomp.git (fetch) upstream https://github.com/collaborativebioinformatics/isocomp.git (push)

On your machine, cd into your local repository, git checkout the development branch, and make sure it is up-to-date with the upstream (ie the original) repository.

NOTE: if you branch, in general make sure you branch off the develop repo, not main!

Then (assuming poetry is installed already), do:

bash $ poetry install

This will install the virtual environment with the dependencies (and the dependencies' dependencies) listed in the pyproject.toml.

Adding dependencies

To add a development dependency (eg, mkdocs is not something a user needs), use poetry add -D <dependency> this is equivalent to pip installing into your virtual environment with the added benefit that the dependency is tracked in the pyproject.toml.

To add a deployment dependency, just omit the -D flag.

Writing code

Do this first!

bash $ pip install -e .

This is an 'editable install' and means that any change you make in your code is immediately available in your environment. NOTE: If you happen to see a Logging Error when you run the pip install -e . command, you can ignore it.

If you use vscode, this is a useful plugin which will automatically generate docstrings for you. Default docstring format is google, which is what the scripts we currently have use. This is an example of what a google formatted docstring looks like:

```python def getallwindows(genedf:pd.DataFrame, bpdf:pd.DataFrame) -> pd.DataFrame: """From gene boundaries and 100 bp nonzero coverage windows, produce a merged window df

Args:
    gene_df (pd.DataFrame): one window per gene, > 0.05 avg coverage
    bp_df (pd.DataFrame): one window per 100 bp, > 0.05 avg coverage

Returns:
    pd.DataFrame: merged windows df
"""
...

```

In the function definition, the type hints of the arguments (eg gene_df:pdDataFrame) are not required, but if you include them, autoDocs will automatically generate the data types in the docstring skeleton, also, which is nice. The -> <datatype> at the end of the function definition is the return data type.

Tests

Unit tests can be written into the src/tests directory. There is an example in src/tests/test_isocomp.py. There are a couple other examples of tests -- ie for logging and error handling -- here, too.

Build

It's good to intermittently build the package as you go. To do so, use poetry build which will create a .whl and .tar.gz (dist is already included in the gitignore). You can 'distribute' these files to others -- they can be installed with pip or conda -- or use them to install the software outside of your current virtual environment.

Documentation

If you would like to write documentation (ie not docstrings, but long form letters to your adoring users), then this can be done in markdown or jupyter notebooks (already added as a dev dependency) in the docs directory. Add the markdown/notebook document to the nav section in the mkdocs.yml and it will be added to the menu of the documentation site. Use mkdocs serve locally to see what the documentation looks like. mkdocs build will build the site in a directory called site, which is in the .gitignore already. Like poetry build it is a good diea to do mkdocs build intermittently as you write documentation. Eventually, we'll use mkdocs gh-deploy to deploy the site to github pages. Maybe if we get fancy, we'll set up the github actions to build the package on mac,windows and linux OSes on every push to develop, and rebuild the docs and push the package to pypi on every push to main.

Computational Resources / Operation

On a 16 core cloud instance with three GTF files (replicates) from Genome in a Bottle subject HG002, find_unique_isoforms took ~15 minutes and ~3 to 5 GB of RAM