{"paper_id":"03398bf4-1b4d-4323-b706-742bb3f80bc6","body_text":"Biosurfer to track protein isoform variation \nBiosurfer for systematic tracking of regulatory mechanisms leading to protein isoform diversity 1 \n 2 \nMayank Murali 1, Jamie Saquing 2, Senbao Lu 6,7, Ziyang Gao 6,7, Ben Jordan 2, Zachary Peters Wakefield 3 \n8,9, Ana Fiszbein 8,9, David R. Cooper 2, Peter J. Castaldi 10,11, Dmitry Korkin 6,7, Gloria Sheynkman 2,3,4,5,*  4 \n 5 \n1 Broad Institute of MIT and Harvard University, Cambridge, MA, USA 6 \n2 Department of Molecular Physiology and Biological Physics, University of Virginia, Charlottesville, 7 \nVA, USA 8 \n3 Department of Biochemistry and Molecular Genetics, University of Virginia, Charlottesville, VA, USA 9 \n4 Center for Public Health Genomics, University of Virginia, Charlottesville, VA, USA 10 \n5 UVA Cancer Center, University of Virginia, Charlottesville, VA, USA 11 \n6 Bioinformatics and Computational Biology Program, Worcester Polytechnic Institute, Worcester, MA, 12 \nUSA 13 \n7 Computer Science Department, Worcester Polytechnic Institute, Worcester, MA, USA 14 \n8 Bioinformatics Program, Boston University, Boston, MA, USA 15 \n9 Department of Biology, Boston University, Boston, MA, USA 16 \n10 Channing Division of Network Medicine, Department of Medicine, Brigham and Women’s 17 \nHospital, Boston, MA, USA 18 \n11 Division of General Medicine and Primary Care, Department of Medicine, Brigham and Women’s 19 \nHospital, Boston, MA, USA 20 \n*Correspondence: gs9yr@virginia.edu 21 \n 22 \n 23 \n  24 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n2 \nABSTRACT 25 \nLong-read RNA sequencing has shed light on transcriptomic complexity, but questions remain about the 26 \nfunctionality of downstream protein products. We introduce Biosurfer, a computational approach for 27 \ncomparing protein isoforms, while systematically tracking the transcriptional, splicing, and translational 28 \nvariations that underlie differences in the sequences of the protein products. Using Biosurfer, we analyzed 29 \nthe differences in 32,799 pairs of GENCODE annotated protein isoforms, finding a majority (70%) of 30 \nvariable N-termini are due to the alternative transcription start sites, while only 9% arise from 5’ UTR 31 \nalternative splicing. Biosurfer’s detailed tracking of nucleotide-to-residue relationships helped reveal an 32 \nuncommonly tracked source of single amino acid residue changes arising from the codon splits at 33 \njunctions. For 17% of internal sequence changes, such split codon patterns lead to single residue 34 \ndifferences, termed “ragged codons”. Of variable C-termini, 72% involve splice- or intron retention-35 \ninduced reading frameshifts. We found an unusual pattern of reading frame changes, in which the first 36 \nframeshift is closely followed by a distinct second frameshift that restores the original frame, which we 37 \nterm a “snapback” frameshift. We analyzed long read RNA-seq-predicted proteome of a human cell line 38 \nand found similar trends as compared to our GENCODE analysis, with the exception of a higher 39 \nproportion of isoforms predicted to undergo nonsense-mediated decay. Biosurfer’s comprehensive 40 \ncharacterization of long-read RNA-seq datasets should accelerate insights of the functional role of protein 41 \nisoforms, providing mechanistic explanation of the origins of the proteomic diversity driven by the 42 \nalternative splicing. Biosurfer is available as a Python package at https://github.com/sheynkman-43 \nlab/biosurfer. 44 \n 45 \nKeywords: LRS Special Issue, Alternative splicing, long-read sequencing, protein isoforms, GENCODE, 46 \nprotein sequence, alternative transcriptional start site (altTSS), reading frame shift, intron retention, 47 \nnonsense mediated decay (NMD), open reading frame (ORF), transcriptional start site (TSS), 48 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n3 \ntranscriptional termination site (TSS), poison exon, Transcription Initiation Site (TIS), sequence 49 \nalignment 50 \n 51 \nINTRODUCTION 52 \nThrough the isoform diversifying mechanisms of alternative transcription, splicing, and 53 \npolyadenylation, nearly every human gene can produce multiple protein products, with ~20K genes 54 \ngiving rise to at least 180K annotated isoforms (Frankish et al. 2023). The pathway from gene to protein 55 \nis marked by several regulatory mechanisms that are highly tuned across development and cell states, 56 \nwith disruption of this regulation producing aberrant isoforms that lead to pathophysiological states such 57 \nas cancer and cardiovascular disease (Cooper et al. 2009). Hence, approaches are needed to systematically 58 \ncharacterize the upstream regulatory causes and functional impacts of such protein isoform sequence 59 \nchanges. 60 \nTranscript and, by extension, protein isoform diversity may now be globally characterized at 61 \ngreat depth for individual samples (Glinos et al. 2022; Reese et al. 2023). Transcript diversity can be 62 \nreadily characterized by long read RNA-Seq, which employs single molecule sequencing of individual 63 \ncDNA or RNA molecules to determine the sequence across the entire length of spliced transcripts (Sharon 64 \net al. 2013; Workman et al. 2019; Pardo/i1Palacios et al. 2021; Tian et al. 2021; Joglekar et al. 2023), with 65 \nplatforms from Oxford Nanopore and PacBio being most commonly used (Clarke et al. 2009; Eid et al. 66 \n2009). Long read RNA sequencing captures long range connectivity between multiple exons of a 67 \ntranscript and can reveal complex splice patterns unattainable by short read sequencing (De Paoli /i1Iseppi 68 \net al. 2021), including dependencies across distal splicing events (Anvar et al. 2018) and alternative 5’ 69 \nand 3’ transcript usage. Given this readily characteri zed complexity afforded by long read sequencing, a 70 \nnatural question is the extent to which such variations of the transcriptome lead to functional effects of the 71 \nproteome. Towards this goal, a necessary step is defining the potential proteome.  Both our group and 72 \nothers have reported methods in which long-read-derived transcript sequences serve as templates for 73 \npredicting full-length protein isoform sequences, thereby providing a global snapshot of potential protein 74 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n4 \nisoforms expressed in a particular biological cond ition (Miller et al. 2022; Veiga et al. 2022; Abood et al. 75 \n2023). 76 \nThe complexity of alternative splicing (AS)—which for practical purposes in this manuscript we 77 \ndefine here as all transcriptional variations, including alternative transcription start sites (TSS) and 78 \ntranscription termination sites (TTS)—can be observed and characterized at different levels: RNA 79 \ntranscript, open reading frames (ORFs), and finally protein sequences (Reixachs /i1Solé and Eyras 2022). 80 \nThe complex interplay between the changes occurring to AS variants and their ORF and protein products 81 \nare hard to characterize and quantify. Changes in mRNA sequence may lead to non-linear or traditionally 82 \nuntracked variations. For example, a subtle splicing event could lead to a reading frame shift and thus to 83 \nmore dramatic changes to the C-terminus of the protein than the originating small change at the RNA 84 \nlevel would suggest. Or, AS could occur at codon boundaries, leading to altered amino acid (AA) 85 \nidentities of codons that technically overlap in genome-space but are differentially “split” across exon-86 \nexon junctions. And complex interplay may also be observed between transcriptional variations and ORF 87 \nchoice, as alternative 5’ transcription or AS could lead to differentially availability of initiator codons, 88 \ndelimiting start codon choice co-translationally. 89 \nWe argue that rather than being an esoteric exercise, the ability to characterize all potential 90 \ninterplay of RNA-protein variation is critical for fully elucidating the transcript and proteomic diversity 91 \nencoded within long-read RNA-seq datasets. As long-read RNA-seq approaches are increasingly adopted 92 \nin large-scale studies of hundreds of samples (Glinos et al. 2022; Reese et al. 2023) and are maturing into 93 \nstable tools being adopted by the community (Pardo /i1Palacios et al. 2021), extracting all sources of 94 \nbiological molecular diversity is critical. Such interplay cannot be characterized by comparison of 95 \nisoforms using just one modality, such as transcript-focused annotation tools like SQANTI or Matt 96 \n(Tardaguila et al. 2018; Gohr and Irimia 2019), or conventional protein sequence alignment tools like 97 \nClustalW (Chenna et al. 2003) or BLAST (Altschul et al. 1990). R ecently, multi-modal comparisons have 98 \nbeen reported. For example, ORFanage is an approach for large-scale annotation of ORFs across 99 \npredicted transcripts in the CHESS database the main focus being optimal selection of ORFs based on 100 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n5 \nprotein-alignments (Varabyou et al. 2023). Other tools geared towards the phylogenetics community have 101 \ndeveloped frame-aware alignment in which the alignment scoring system is includes penalties for frame-102 \nshift-inducing gaps (Evans and Loose 2015; Ranwez et al. 2011; Jammali et al. 2022) . However, these 103 \ntools do not comprehensively elucidate the interplay between transcript and protein variation. 104 \nTo characterize how AS impacts protein sequence, several bioinformatic tools and databases have 105 \nbeen developed, such as VastDB, ASPicDB, ExonOntology, and DIGGER (Martelli et al. 2011; Tapial et 106 \nal. 2017; Tranchevent et al. 2017; Louadi et al. 2021). Tools such as tappAS and IsoTV annotate how 107 \nprotein isoform sequence and potential functional differences (de la Fuente et al. 2020; Annaldasula et al. 108 \n2021). For example, tappAS is a Java application for quantifying differential isoform usage but also to 109 \nfunctionally annotate such isoforms, using the module IsoAnnot (de la Fuente et al. 2020). IsoAnnot maps 110 \nprotein features (e.g., Pfam domains) across isoforms of a gene, and determines how splicing leads to 111 \npartial or full removal of protein features, indicating potential changes to molecular functions. Despite 112 \nexistence of these tools, of need is the ability to systematically capture all possible effects to protein 113 \nisoforms, with the accompanying information of the underlying complex RNA-protein relationships.  114 \nWe developed Biosurfer, a computational pipeline that tracks simultaneously the changes at all 115 \nthree levels, to understand the impact of alternative splicing on transcriptome, ORFeome, and proteome 116 \ndiversity. Biosurfer computes details not immediately apparent from genome annotation files or manual 117 \ninspection in genome browsers, such as how between isoforms of the same gene the frame of translation 118 \nand codon topology influences amino acid sequence identity changes. In order to accomplish this multi-119 \nlayered comparison, a genome is used as a “scaffold ” to exactly position all nucleotides, codons, AA 120 \nelements, and associate with each element local and context-dependent attributes. The resulting data 121 \nstructure of three distinct yet interlinked layers inform on the upstream biological mechanisms leading to 122 \nprotein sequence changes. To demonstrate the utility of Biosurfer, we globally characterize variation 123 \nacross protein isoforms in the reference human annotation (GENCODE) and protein isoforms predicted 124 \nfrom long-read RNA-seq of a human stem cell. This characterization includes comprehensive elaboration 125 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n6 \nof all sources of N-terminal, internal, and C-terminal protein variations observed, highlighting 126 \nmechanisms of RNA-protein interplay. 127 \n 128 \nMETHODS 129 \nBiosurfer package for isoform analysis and visualization 130 \nBiosurfer is a computational pipeline that performs a multi-layered comparison between a pair of 131 \nisoforms, in which differences at three different levels: RNA (nucleotides, nt, ORF (codon), and protein 132 \n(AA residues, AA) are simultaneously tracked. The developed data structure enables not only the 133 \ncomparison of AAs, but tracking frameshifts, patterns of codon splitting at junctions, and the attendant 134 \nupstream nucleotide and codon differences that explain AA changes. Such tracking aids in the systematic 135 \nannotation of explanatory mechanism(s) underlying AA residue changes, such as whether a substitution 136 \nof a stretch of AA residues is due to alternative splicing or a frameshift (Supplementary Figure S1). 137 \nThe Biosurfer pipeline is organized into the following three stages ( Figure 1). First, an SQLite 138 \ndatabase is populated with detailed information on each isoform at the transcript, ORF, and protein 139 \nproduct levels. For each isoform, the required inputs are (i) a transcript FASTA, (ii) a protein FASTA, 140 \nand (iii) a matching GTF with both Exon and CDS features. The inputs can either be extracted from the 141 \nreference annotations (e.g., GENCODE (Harrow et al. 2012)), or they can be user-defined ( e.g., predicted 142 \nprotein isoform sequences from long-read RNA-seq data (Miller et al. 2022)by using ORF callers such as 143 \nCPAT (Wang et al. 2013) or Transdecoder (Haas et al. 2013)). Second, mu ltilayered isoform-level 144 \nalignments are generated, and each alignment is represented as three key data structures: t-blocks, c-145 \nblocks, and p-blocks corresponding to the view of the aligned isoforms at the transcript, codon sequence, 146 \nand protein sequence levels, respectively ( Supplementary Figure S2 ). Third, all information is 147 \nsummarized in tabular format, to facilitate analysis of the transcription-level and codon-level mechanisms 148 \ndriving the proteomic diversity. In addition, a visual representation of such mechanisms integrated with 149 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n7 \nthe protein-level view of the isoforms can be output as png files.150 \n 151 \nFigure 1: Biosurfer for analysis and visualization of protein isoform sequence differences.  152 \nBiosurfer analyzes protein isoforms from reference annotations (e.g., GENCODE) or proteins predicted 153 \nfrom long-read RNA-seq data. Analysis initiates with the creation of an SQLite database is populated 154 \nwith isoform-relevant elements. Biosurfer performs a multi-layered comparison of transcript-, codon-, and 155 \nprotein-level differences between pairs of protein isoforms. Variable regions as well as their associated 156 \nannotations are output in tabular format and visualization files, which includes protein-relevant details 157 \nsuch as the frame of translation. Note that the terms “Match”, “Deletion”, etc. represent very different 158 \ncomparisons depending on the biological layer. GTF=Gene Transfer Format; iPSCs=induced pluripotent 159 \nstem cells. 160 \n 161 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n8 \nData structures for presenting transcript-, ORF-, and protein-level isoform alignments in Biosurfer 162 \nBiosurfer compares transcript-, ORF-, and protein-level sequences for each gene. One isoform is 163 \nselected as reference and the other isoforms are denoted as alternative. When aligning each alternative 164 \nisoform to the reference, the matched regions are referred to as matched blocks. The remaining regions on 165 \neach sequence that are not matched blocks are referred to as unmatched blocks. See Supplementary File 166 \n1 for step-by-step schema of the comparison process, which is summarized below. 167 \nTranscript blocks (t-blocks). T-blocks represent subsegments of the transcript sequence that are shared or 168 \nunique to the reference or alternative isoforms. Transcript-level differences are determined by analyzing 169 \nthe transcript-to-genome coordinates alignment (GTF file) provided as an input by the user. Specifically, 170 \nBiosurfer defines the aligned exonic regions that are shared or unique to each isoform. The resulting 171 \nranges, called t-blocks, are categorized as Match, Deletion, or Insertion t-blocks. Deletion or Insertion t-172 \nblocks are further annotated with the associated biological mechanism leading to the transcript nucleotide 173 \nchange, e.g., alternative transcriptional start site, splicing event, or polyadenylation. The alternative 174 \nsplicing events are then further classified into four basic types: retained intron, alternative donor, 175 \nalternative exon, or alternative acceptor. 176 \nCodon blocks (c-blocks). C-blocks represent a codon-centric layer defined through the comparison of the 177 \nprotein coding regions of transcripts, i.e., open reading frames (ORF), between two isoforms. The c-block 178 \ndata structure is the most complex, but critical, layer in Biosurfer that connects information between the 179 \ntranscript and protein layers. 180 \nFor ultimate granularity and precision, ORFs are compared based on the alignment of codons that 181 \noverlap in the genome space, in which one codon in the reference isoform is compared with another 182 \ncodon in the alternative isoform (see Supplemental Methods for additional details). First, codons across 183 \nthe two isoforms are “ paired” based on their mutual positions and base overlap: in a basic scenario, the 184 \ntwo codons are identical and their positions match ( Table1, Supplemental Figure S3A ); in a more 185 \ncomplex scenario, the codons are split and only partially overlapped (1 or 2 bases, see Table1, 186 \nSupplementary Fig S3B and C). In all other cases, the codons will be designated as “unpaired” ( Table1, 187 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n9 \nSupplementary Figure S3B).  To keep track of unpaired codons within the data structure, they are linked188 \nto a “placeholder” codon, which serves to maintain consistency of the c-block structures.  Once aligned,189 \neach single codon pair is categorized with respect to multiple attributes that explain their relationship.190 \nOverall, the paired and unpaired codons are classified into 9 categories based on their translation status,191 \nframeshift status, and codon topology (Table 1). 192 \n193 \nTable 1: Categories of overlapping codon pairs that are the basis for codon blocks (c-blocks). 194 \n 195 \nProtein blocks (p-blocks).  P-blocks represent a protein-centric layer defined through t- and c- block196 \nguided comparisons of two protein isoforms, thus focusing solely on the relationships between the AA197 \nresidues of the respective protein isoforms. This decision was motivated by the fact that in most cases,198 \nfunctional differences exhibited by an alternative protein isoform are attributable to differences in protein199 \nrather than nucleotide sequences between the alternative and reference isoforms , and these differences200 \ncould be conceptually decoupled from the knowledge of the upstream (ORF- or transcript-level) residue -201 \naltering mechanisms.  202 \ned \nd, \nip. \nus, \n \nck \nA \nes, \nin \nes \n-\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n10 \nThe p-blocks for a pair of isoforms are determined by translating c-blocks defined at the previous 203 \nlevel. The resulting data structure represents subsegments of contiguous stretches of AA residues from 204 \neach protein sequence. Each subsegment could be as short as a small translated portion of exon (even a 205 \nsingle residue, in the case of NAGNAG splicing), but it can also span multiple exons of a gene. Because 206 \nthe actual matching happened at the two previous levels (t-blocks and c-blocks), no alignment is required 207 \nat the p-block level. The resulting protein subsegments can be either (1) fully matched (100% sequence 208 \nidentity, across all residues in the subsegment), or (2) mismatched, where a subsegment in one protein 209 \nsequence will be aligned against a gapped region in the other sequence. Thus, a p-block represents either 210 \na matched pair of protein subsequences or a subsequence matched against a gapped region. P-blocks are 211 \nthen classified as Match, Insertion, Deletion, or Substitution changes. Substitution p-blocks must arise 212 \nfrom a combination of insertion/deletion/frameshift events found at the c-block level. 213 \nOverall, these p-block changes are agnostic to the upstream mechanisms, e.g., at the transcript or 214 \nORF levels; however, at the same time, the corresponding upstream mechanisms can be retrieved and 215 \nanalyzed within Biosurfer, unlike using traditional protein aligners.  216 \n 217 \nAnalysis of GENCODE isoforms using Biosurfer 218 \nThe principal and alternative isoforms were defined using APPRIS annotation of genes in GENCODE 219 \n(Rodriguez et al. 2013). First, the set of APPRIS isoforms is identified for each gene from the input 220 \ngenome data (GENCODE v42, basic annotation (Frankish et al. 2021)) by extracting the key transcript 221 \nfeatures, such as 'transcript_id', 'transcript_name', and the associated APPRIS ‘tag’ (Rodriguez et al. 222 \n2013). Second, the transcript's APPRIS status is determined as 'principal', 'alternative', or ‘none’, based on 223 \nfirst rank-ordering of transcripts based on APPRIS tag and setting the transcript with the highest APPRIS 224 \nvalue as ‘principle’ and all other transcripts as ‘alternative’. We did not consider further genes with only 225 \none annotated isoform and lacking an APPRIS tag (‘none’). Subsequently, we compiled a structured 226 \ndataset for each transcript, encompassing identifiers, gene associations, strand orientation, and APPRIS 227 \nstatus. During the annotation, the strand orientation is accounted for when necessary. 228 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n11 \n 229 \nLong read RNA-seq analysis of WTC-11 cell line 230 \nTotal RNA from WTC-11 cells was extracted using the RNeasy Kit (Qiagen) and analyzed on an 231 \nAgilent Bioanalyzer. We observed a RNA concentration of 30 ng/uL with the RNA Integrity Score (RIN) 232 \nof 9.9. As described previously, (Mehlferber et al. 2022) cDNA was synthesized from the extracted RNA 233 \nand the Iso-Seq Express Kit SMRT Bell Express Template prep kit 2.0 (Pacific Biosciences) was used on 234 \na Sequel II system to obtain long-read sequence information and output Circular Consensus (CCS) reads. 235 \nData is available at the Sequence Read Archive: SRR18130587 and previously published (de Souza et al. 236 \n2022). 237 \nWe analyzed the WTC-11 data with a proteogenomics Nextflow pipeline we previously 238 \ndeveloped (Miller et al. 2022) (Mehlferber et al. 2022). The output CCS reads from long-read sequencing 239 \nwere processed with Iso-Seq3 and SQANTI3 (version 1.3) for transcript isoform classification and quality 240 \nassessment. The CPAT(Wang et al. 2013) algorithm was used to predict Open Reading Frames (ORFs), 241 \nwhich were grouped into protein isoforms. 242 \n 243 \nBiosurfer is implemented as a Python package freely available at GitHub repository: 244 \nhttps://github.com/sheynkman-lab/biosurfer. The analysis code is available at 245 \nhttps://github.com/sheynkman-lab/biosurfer_analysis. All necessary input files and intermediate and final 246 \noutput files from Biosurfer analysis are uploaded to Zenodo at: https://zenodo.org/records/10822882.  247 \n 248 \nRESULTS 249 \nCharacterization of altered protein regions in the human annotation (GENCODE) 250 \nHere, we demonstrate the utility of Biosurfer through a genome-wide analysis of protein isoforms 251 \nin the GENCODE annotation (basic annotation, version 42 (Frankish et al. 2021)). We analyzed 35,083 252 \nreference-alternative protein isoform pairs across 11,815 genes ( Figure 2A ). Each pair consists of the 253 \nreference protein isoform for a gene—the APPRIS principal isoform (Rodriguez et al. 2013)—and an 254 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n12 \nalternative protein isoform. The number of isoform pairs correspond to the number of “alternative” 255 \nisoforms for a gene ( Supplementary Table S1). Genes lacking an APPRIS annotation (8,166 of 19,981) 256 \nwere excluded from the analysis (Supplementary Table S2). 257 \nGlobally, we found a total of 44,326 altered protein regions with an average of 1.3 altered regions 258 \nper isoform. Note that we are using the term altered protein region interchangeably with p-block (see 259 \nMethods). Altered protein regions are contiguous regions of altered protein sequence relative to the 260 \n“reference” (i.e., APPRIS principal) protein, which includes p-block Insertions, Deletions, and 261 \nSubstitutions. A majority (79%, 27,672 of 35,083) of isoform pairs contain a single altered region (Figure 262 \n2B). Still, 7,411 alternative isoforms contain two or more discontinuous altered regions. Notably, some 263 \nisoform pairs exhibited an extreme number of regions. For example, 14 regions are found for proteins of 264 \nDNAH14 (Reference: DNAH14-220, Alternative: DNAH14-211), which is explained by its extremely 265 \nlarge number of exons (86 exons for DNAH14-220). 266 \nAmong the 44,326 altered protein regions ( Figure 2C ), the median number of affected amino 267 \nacids (lost or gained due to a protein insertion, deletion, or substitution) is 49 AA, with the first and third 268 \nquartiles containing 21 AA and 128 AA, respectively. Among the altered regions, 14% (6,189) are 269 \ninsertions, 47% (20,780) are deletions, and 39% (17,357) are substitutions (Figure 2D ). Full annotations 270 \nfor these altered regions at the protein and codon-level, representing the Biosurfer output for protein-271 \nblocks (p-blocks) and codon-blocks (c-blocks) can be found in Supplementary Table S3, S4.  272 \nThe lengths of p-block insertions tend to be shorter than the length of deletions (p < 2.2e-16, 273 \nMann-Whitney U test) or substitutions (p < 2.2e-16, Mann-Whitney U test) ( Figure 2E-H ). Since the 274 \ndeletion or insertion status of a protein region is dependent on which isoform is denoted as the reference, 275 \nthis trend may reflect the tendency that longer isofor ms are more likely to be defined as the “reference” 276 \nisoform, which may be biologically driven or influenced by genome annotation guidelines (Rodriguez et 277 \nal. 2013; O’Leary et al. 2016; The UniProt Consortium 2017; Frankish et al. 2021; Varabyou et al. 2023). 278 \nWe found a similar trend for p-block substitutions; the lengths for substituted regions tended to be longer 279 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n13 \nfor the sequence affected in the reference ( Figure 2G ) versus the alternative ( Figure 2H ) isoform;280 \nalthough, for 3,263 (18% of cases), the affected regions in alternative isoforms can be longer (Figure 2I). 281 \n282 \nFigure 2:  Characterization of altered protein regions (Biosurfer p-blocks) across the GENCODE283 \nannotated human proteome. ( A) Schematic of genes with one or two alternative protein isoforms and284 \ndistribution of the number of alternative protein isoforms per gene. (B) Schematic of the altered protein285 \nregions (here, highlighted in pink and b lue), displayed relative to the underlying transcript structures. The286 \nproteome-wide distribution of the number of affected protein regions observed per alternative isoform.287 \n(C) Distribution of the length of altered protein regions across the annotated proteome. Differences 288 \n greater than 600 AAs are not included (2,347 cases, 5.3% of the data). These altered protein regions289 \ninclude cases of 1) deleted regions, 2) inserted regions, and 3) the region (in the reference isoform) in290 \nwhich one polypeptide subsegme nt is substituted for another. In other words, this distribution (C)291 \nrepresents an aggregation of the distributions shown in panels E- G. (D) Fraction of altered protein regions292 \nm; \n \n \nE \nnd \nin \nhe \nm. \nns \n in \nC) \nns \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n14 \naffected by insertions, deletions, and substitutions. (E-H) Distribution of the lengths of altered protein 293 \nregions for (E) deletions, (F) insertions, (G) substituted region in the reference isoform, and (H) 294 \nsubstituted region in the alternative isoform. (I) Comparison of the lengths of altered regions in the 295 \nreference versus alternative isoforms for substitutions. 296 \n 297 \n 298 \nAnalysis of N-terminal protein variations 299 \nVariations in the N-termini are found in 28% (12,504 of 44,326) of all possible reference-300 \nalternative isoform pairs, corresponding to 5,872 genes (Supplementary Table S5). We examined for the 301 \nN-terminal variations the explanatory mechanism, including alternative transcription start sites (TSSs), 302 \nalternative splicing, and alternative translation initiation sites (TISs). 303 \nAll cases of variable N-termini involve two initiation codons (AUGs), one upstream and one 304 \ndownstream, relative to the genome. A major category we first observed are those in which the N-305 \nterminus is different due to start codons that are mutually exclusively present across the two transcript 306 \nisoforms (Figure 3A ). Specifically, the start codon present in the reference isoform is absent from the 307 \nalternative isoform, and vice versa. These “mutually exclusive start codons” or MXS were observed for 308 \n3,123 reference-alternative isoform pairs ( Figure 3A). MXS may arise either from an alternative TSS or 309 \nfrom alternative splicing in the 5’ UTR. Strikingly, we found that nearly all (99%, 3,097 of 3,123) cases 310 \nare caused by alternative TSS usage ( Figure 3A, hatched region of the bar), with only a small, but non-311 \nzero, fraction (1%, 26 of 3,123) of MXS cases arise from splicing of the 5’ UTR, in which splicing 312 \nregulation is influencing N-terminal usage. An example of TSS-driven MXS for PRKACA is shown in 313 \nFigure 3B (pair, PRKACA-201 and PRKACA-202). 314 \nThe second category is when the upstream start codon is transcribed in only one of the two 315 \nisoforms, but the downstream start codon is present in both transcripts. We refer to this scenario as shared 316 \ndownstream starts (SDS), of which there were 6,878 cases ( Figure 3A). Like with cases of MXS, SDS 317 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n15 \narises primarily from alternative TSS usage (80% of  cases) versus 5’ UTR alternative splicing (20% of 318 \ncases). 319 \nMXS and SDS are common patterns underlying alterations of the N-termini or protein, driven by 320 \ndifferential availability of initiator codons in the mature transcript. We asked if there may be differences 321 \nin length of such N-terminal alterations for SDS versus MXS events. Measuring the differential length of 322 \nthe affected N-terminal regions between reference-alternative isoform pairs, we found that, on average, 323 \nSDS tends to affect a greater proportion of protein length, as compared to MXS ( Figure 3C , 324 \nSupplementary Figure S4; p = 2.8e-148, Mann-Whitney U test). The larger differences in length driven 325 \nby SDS could be explained by cases in which transcription is initiated from internal sites of the gene, 326 \ngiving rise to an ORF that corresponds to a subsequence of the ORF in the other isoform, theoretically 327 \nproducing a truncated C-terminal-containing subsequence of the full-length protein. 328 \nRecently, a mechanism related to SDS was described in which internal exons (not the 5’ most 329 \nexon, first exon, of a transcript) in one transcript can be immediately downstream of a DNA element of 330 \nnovel promoter activity and thus serve as the first transcribed exon in other isoforms (Fiszbein et al. 331 \n2022). Such so-called hybrid exons thus can operate as both sites of transcription initiation and alternative 332 \nsplicing, in effect, swapping their roles depending on the regulatory context. Of the cases of SDS, we 333 \nobserved that ~22% correspond to these hybrid exon swaps ( Figure 3D, Supplementary Table S6). The 334 \nfunctional consequences of hybrid exon usage are not well understood; however, one potential function 335 \ncould be the production of protein with a truncated N-terminus, which could remove signal peptide 336 \nsequences or binding domains (Kelemen et al. 2013).  337 \nIn addition to MXS and SDS, wherein start codon availability is controlled through differential 338 \ntranscription, we also observed many cases in which both upstream and downstream start codons co-339 \noccur in one or both transcript isoforms of a pair. In these cases, the choice of start codon may be 340 \ninfluenced by co-translational regulation, e.g., ribosome initiates translation at alternative initiation sites 341 \n(altTIS). 342 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n16 \nFor 2,054 cases, we found altTIS, which we classified as instances of a mutually shared start 343 \n(MSS) codon. We also found 449 cases in which the upstream start codon is present in both isoforms, but 344 \nthe downstream start codon is only present in one isoform. 345 \nIn such cases of shared upstream start (SUS) ( Figure 3A ), the explanatory co-translational 346 \nmechanism is not as clear, based on the ribosomal scanning model of translation initiation (Kozak 1978). 347 \nThe upstream start codon present in both isoforms would need to be bypassed by the ribosome only in the 348 \nalternative isoform. Therefore, some of the SUS annotations may need to be validated or may be 349 \nerroneous ORF calls, as early ORF prediction workflows attribute higher scores to longer ORFs (Wang et 350 \nal. 2013; Varabyou et al. 2023), under-annotating ORFs that utilize the upstream (annotated) start codon 351 \nbut is much shorter than the reference due to a reading frame shift (Wang et al. 2013). 352 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n17 \n 353 \n 354 \nFigure 3:  Analysis of mechanisms underlying variable N- terminal proteins across the GENCODE355 \nannotated human proteome.  (A) Distribution of the types of alternative N- terminal regions, classified356 \nbased on presence and translational status of the start codon. Hatches denote the fraction of alternative N -357 \nterminal regions associated with alternative transcription start sites, as opposed to 5’ UTR alternative358 \nsplicing. (B) Biosurfer output of altered N-terminal regions for PRKACA gene that undergo MXS (light359 \ngreen arrows) and SDS (dark green arrows). In the example of MXS, the yellow bars above and below 360 \nE \ned \n-\nve \nht \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n18 \n PRKACA-202 (alternate) transcript indicate the N-terminal ranges that differ between the reference and 361 \nalternative isoform. The yellow bar above PRKACA-202 shows the N-terminal protein sequence that is 362 \nspecific to the reference (PRKACA-201) and the bar below the transcript indicates the N-terminal region 363 \nspecific to the alternative isoform. The Biosurfer bars span the intronic lines between exons, but intronic 364 \nregions do not contribute to the protein sequence differences. In the example of SDS, the pink Biosurfer 365 \nbar above the transcript of PRKACA-203 represents the range of transcript sequence that is translated in 366 \nthe reference, but not translated in the alternative isoform. (C) Scatterplot of the length of affected N-367 \nterminal variation in the reference versus alternative, faceted by mutually exclusive starts (MXS) or 368 \nshared downstream start (SDS) status. An interesting case of an SDS leading to unique N-terminal 369 \nsequence in the reference is caused by usage of a different frame at the initiation of translation, with an 370 \nexample shown for isoforms of the gene FHOD3. (D) Fraction of shared downstream starts (SDS) caused 371 \nby hybrid exon swaps. 372 \n 373 \n 374 \nAnalysis of internal protein variations 375 \nThe internal regions of the protein isoforms account for 43% (19,263 of 44,326) of possible 376 \nreference-alternative isoform pairs, corresponding to 6,673 genes ( Supplementary Table S7 ). A large 377 \nmajority of these regions (80%, 15,444 of 19,263) are caused by single simple splicing events: exon 378 \nskipping, alternative acceptor, alternative donor, or  an in-frame retained intron. As expected, exon 379 \nskipping events are most numerous making up 9,505 of cases. In terms of the general effect on protein 380 \nsequence, most altered regions (70%, 13,453 of 19,263) lead to a deletion or removal of AA residues 381 \n(Figure 4A), and, again, in such cases, exon skipping is most common (51%, 6,908 of 13,453). 382 \nGoing beyond simple splicing events, we observed that 19% of variable internal regions (3,701 of 383 \n19,263) were associated with multiple events, which we refer to as compound, or linked, splicing events. 384 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n19 \nWe found that 56% (2,099 of 3,701) compound events involve multi-exon skipping, the rest being 385 \ncombinations of alternative donor/acceptor sites with exon skipping/inclusion (Figure 4B). 386 \nComplex nucleotide to amino acid relationships that affect internal protein sequence 387 \nPrevious studies of the impact splicing on proteins have typically focused on cases in which 388 \ndifferential splicing of transcript regions directly corresponds to changes in protein sequence 389 \n(Reixachs/i1Solé and Eyras 2022). However, in many instances, there is not a simple one-to-one 390 \nrelationship between nucleotides in a transcript and the corresponding amino acid identities in the 391 \nencoded protein. Such complexity alters the protein sequence in non-intuitive ways. Using the detailed 392 \ncodon tracking afforded by Biosurfer, we systematically characterized the protein-level impact of 393 \nvariations not commonly described: codons that span junctions and unusual reading frame shifts. 394 \nTo characterize differentially split codons, we examined all paired codons (see c-block section in 395 \nMethods) that are split across junctions and determined the identity of the associated AAs 396 \n(Supplemental Figure S3 ). Across all internal altered protein regions, we found 17% (3,213 of 19,263) 397 \nof regions that are flanked by one or more split codon pairs that encode different amino acid residues 398 \n(Figure 4C , also see Table 1 ). These so-called ragged codons affect a single residue and are always 399 \nadjacent to an altered protein region. Interestingly, while split codons are not particularly enriched by 400 \nsplice event type (Chi-square test: p-value = 5.54e-99), we found that ragged codons are more frequently 401 \nfound in protein insertions compared to deletions or substitutions (Chi-square test: p-value = 4.11e-106).  402 \nFor a majority of splice-driven frame shifts, the shifted frame is maintained to the end of the 403 \nprotein, leading to a protein variation affecting the C-terminus. Such frameshifted proteins could lead to 404 \ntruncated proteins or destabilization of the transcript via mechanisms such as nonsense mediated decay 405 \n(NMD). Surprisingly, we found an un commonly characterized pattern of  successive readin g frame shifts 406 \nthat exclusively affects the internal residues of a protein. In these cases, the alternative isoform’s reading 407 \nframe is shifted due to one splicing event, but then shifts back into register of the reference frame due to a 408 \nsecond, independent splicing event. We refer to these events as “snapback” frameshifts, as there is a 409 \nreturn back to the original reading frame. Snapback frameshifts have been previously observed, such as in 410 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n20 \nthe gene HSF4 , which produces internally frameshifted isoforms that have been demonstrated to exert411 \ndifferent regulatory effects (Tanabe et al. 1999), but the snapback phenomenon generally speaking has not412 \nbeen systematically described. Within GENCODE, we found 118 examples of suc h snapback isoforms413 \nacross 95 genes, including FAS , and PLEKHJ1 (Figure 4D, Supplementary Table S8). What is notable414 \nabout these cases is that the same underlying genomic sequence encodes different amino acid residues,415 \nand genetic mutations could lead to two different residue changes depending on the isoform. 416 \n417 \nFigure 4: Analysis of internal protein altered regions across proteins in GENCODE.  (A) Frequencies418 \nof the categories of the splicing mechanism underlyi ng internal protein sequence changes, split by their419 \nprotein-coding impact (Deletion, Insertion, Substitution). All sequence regions involve a reference -420 \nalternative isoform pair. (B) Frequency of compound splicing events across the altered internal regions.421 \n(C) Proportion of each internal protein region type for which there exists a split codon pattern near its 422 \n boundaries that would cause a single amino acid difference, or “ragged” codon. (D) Examples of423 \nsuccessive frameshifting (snapback” frameshift) that leads to an affected protein region that is wh olly424 \ninternal to the protein, for genes FAS and PLEKHJ1. 425 \n 426 \nAnalysis of C-terminal variations 427 \nert \not \nms \nle \nes, \nies \neir \n-\nns. \nof \nlly \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n21 \nC-terminal changes make up 28% (12,241 of 44,326) of reference-alternative pairs, 428 \ncorresponding to 6,138 genes (Supplementary Table S9). 429 \nTo break down the sources of C-termini variability, we found it useful to distinguish the most 430 \ndirect preceding cause of altered C-termini. In principle, all C-terminal changes must arise from an 431 \nupstream splicing event that influences the termination codon used (notwithstanding post-translational 432 \ncleavage events). However, such changes could be further classified. The altered C-terminus could arise 433 \nfrom alternative terminal (i.e., last) exons, each harboring a different stop codon, so that the splicing event 434 \nmore or less directly influences the stop codon availability (i.e., direct splice-driven events). In other 435 \ninstances, C-terminal changes could arise from a somewhat indirect relationship to the splice event, such 436 \nas when a splicing event causes a translational frameshift, in effect, “revealing” in the other frame a new 437 \nstop codon that is now decoded by the ribosome (i.e., frameshift-driven events) ( Supplementary Table 438 \nS9). 439 \nIn direct splice-driven events, the stop codon availability is dictated by the actively transcribed 440 \nregions that contain the stop codon. We find this scenario for 72% (8,877 of 12,241) of all C-terminal 441 \nvariations (Figure 5A). These variations can be further classified based on the pattern of splicing at the C 442 \nterminus. The first pattern involves an exon extension into the intron region, introducing a premature stop 443 \ncodon. These exon extensions introduce termination, or “EXIT”, make up 3,498 (39.4% of 8,877) cases 444 \n(Figure 5B). The second pattern involves usage of alternative terminal coding exons or “ATE”, making 445 \nup 3,301 (37.1% of 8,877) of cases ( Figure 5B ). Overall, EXIT and ATE changes lead to a shorter C-446 \ntermini in the alternative isoform (distribution shown in Figure 5C). 447 \nEXIT versus ATE reflect how different “modes” of spliceosome regulation could lead to distinct 448 \nC-terminal consequences. In EXIT, the reference-containing donor splice site fails to be spliced in the 449 \nalternative isoform, leading to partial or full intron retention. An example of EXIT is shown in the right 450 \npanel of Figure 5C for the pair, TDRD12-206 and TDRD102-202. In ATE, on the other hand, the 451 \nspliceosome catalyzes splicing at one of two splice site acceptor sites, influencing terminal exon identity 452 \nand thus stop codon used. A well-known pattern of ATE are poison exons, a mechanism by which 453 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n22 \ninclusion of an exon leads to a premature termination codon that either elicits nonsense mediated decay or 454 \ngenerates a truncated protein product (Carvill and Mefford 2020). Poison exons are evolutionarily 455 \nconserved and likely play a role in downregulating gene expression (Lareau et al. 2007). We found 550 456 \n(6.2% of total 8,877) cases of potential poison exons. Other ATE patterns include one in which the 457 \nalternative last exon of the alternative isoform resides in the UTR region of the reference isoform, 458 \nsuggesting that such sequences in the 3’ UTR could dually code both transcript and protein functional 459 \nelements. We found 819 (9% of 8,877) cases of such alternative last exon in UTR (ALE in UTR) ( Figure 460 \n5B). A third ATE variation, referred to as a cut-out splice terminal exon (COSTE), the 5’ end of the last 461 \nexon is shared, but the alternative isoform utilizes a splice site that skips over the remaining portion of the 462 \nlast exon in the reference, thereby creating a different last exon not found in the original reference. We 463 \nidentified 273 cases of this pattern (3% of 8,877) (Figure 5B and Supplementary Figure S5). 464 \nFrameshift-driven events influence stop codon usage somewhat indirectly through shifts in the 465 \ntranslational reading frame. In such cases, a splice-induced reading frame shift causes all downstream 466 \ncodons to be read in a different frame and stop codons are “revealed”, or decoded, by the ribosome. We 467 \nfound 3,364 (28% of 12,241) cases of frameshift-induced C-terminal changes ( Figure 5A, Examples are 468 \nshown in Figure 5D , pairs, GIPC3-201 and GIPC3-202, along with FANCM-201 and FANCM-231). 469 \nGenerally speaking, frameshifts lead to a dramatic shortening of the C-terminal region in the alternative 470 \nisoform; however, we found 549 cases (16% of 3,364) in which the C-terminal region is longer in the 471 \nalternative versus the reference isoform (Figure D). 472 \nWe also observed across all frameshift-driven events a depletion of isoform pairs in which a large 473 \nportion of the reference isoform (e.g., 2,000 AA or longer) is truncated due to a frameshift in the 474 \nalternative isoform (Figure D), a trend not observed in an experimentally predicted proteome (see section 475 \nbelow and Supplementary Figure S9B ), likely representing gene annotation decisions, as dramatically 476 \ntruncating frameshift events would lead to predicted NMD and filtered out or reassigned an NMD biotype 477 \n(Harrow et al. 2012). 478 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n23 \n479 \nFigure 5: Analysis of alternative C-terminal protein sequences in GENCODE v42. (A) Frequency of480 \nalternative C-terminal categories based on splicing or frameshifts being the primary driving factor. (B) 481 \n Distribution of the frequencies for various splice-driven patterns. (C) Scatterplot of the length of splice -482 \ndriven C-terminal variation in the reference versus alternative. An example of this category is observed in483 \nthe TDRD12 gene. TDRD12  undergoes splice- driven alteration causing an alternative terminal exon in484 \nTDRD12-202 to harbor the stop codon (D) Scatterplot of the length of frameshift-driven C- terminal485 \nvariation in the reference versus alternative. Biosurfer plot examples of frameshift- driven category486 \nillustrated in GIPC43 and FANCM genes. 487 \n 488 \n 489 \nCharacterization of altered protein regions across a long-read predicted proteome 490 \nThe process of defining the reference proteome heavily draws from sources of experimental491 \nevidence such as deeply sequenced cell and tissue types, and the protein isoform sequences represent a n492 \naggregate model of the human proteome. Therefore, to characterize potential isoform related protein493 \n \nof \n-\n in \n in \nal \nry \ntal \nn \nin \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n24 \nvariability in a specific biological condition, we employed a “long read proteogenomics” pipeline 494 \n(Kreitzer et al. 2013), generating a proteome predicted from long-read RNA-seq transcript sequences 495 \ncollected from a karyotypically normal human stem cell line (WTC-11). Using Biosurfer, we 496 \ncharacterized the landscape of protein isoforms in the WTC-11 proteome and found 44,962 protein 497 \nisoform pairs across 10,144 genes ( Supplementary Table S10 , p-block and c-block outputs in 498 \nSupplementary Table S11 and S12 ). Assigning the highest expressed transcript as the “reference”, we 499 \ndefined 53,915 altered protein regions ( Supplementary Figure S6 ). Compared to the GENCODE 500 \nanalysis, similar trends were observed for N-terminal ( Supplementary Figure S7 ) and internal region 501 \nvariation (Supplementary Figure S8 , snapback isoforms in Supplementary Table S13 ). Besides these 502 \nsimilar trends, the experimental proteome returned a higher number of C-terminal variations that involved 503 \nintron retention events (EXIT event type, see Figure 5B) as compared to GENCODE isoforms, matching 504 \nearlier findings from EST and cDNA data (Nakao et al. 2005)(Modrek et al. 2001)( Supplementary 505 \nFigure S9). 506 \n 507 \nDISCUSSION 508 \nTo study the functional impact of alternatively spliced protein isoforms, it is critical to track 509 \nprecise differences in protein isoform sequences and link such variations to the upstream explanatory 510 \nmechanisms. However, it is challenging to systematically characterize the full interplay between genomic 511 \nand proteomic variations, which hinders discoveries of novel biological variations represented in a long-512 \nread RNA-seq dataset. We developed Biosurfer, a computational approach, available as a Python 513 \npackage, that systematically extracts protein isoform sequence variations while maintaining the explicit 514 \nlinks to their underlying transcriptional and post-transcriptional mechanisms. 515 \nTo demonstrate the utility of Biosurfer, we characterized protein isoform differences across an 516 \nannotated (GENCODE) and long-read RNA-seq predicted proteome. Using Biosurfer’s interlinked 517 \ntranscript, codon, and protein data structures, we determined the upstream mechanisms explaining 518 \nisoform alterations, uncovering surprising complexity. First, we confirmed past observations of 519 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n25 \nalternative transcription underlying most N-terminal variations (Reyes and Huber 2018), and most 520 \ninternal protein sequence differences arising from single splicing events, such as exon skipping (Wang et 521 \nal. 2008). However, our study goes beyond these known trends. Upon detailed tracking of codons 522 \nflanking these altered protein regions, we found distinct split codon patterns that change the encoded 523 \namino acid residue identity and thus contribute to variation of AA residues. We also found an unusual 524 \nframeshift pattern that involves successive reading frame shifts that leads to a change in protein regions 525 \nthat is entirely internal to the protein, referred to here as \"snapback\" frameshifts. And last, C-terminal 526 \ndifferences are primarily splice-driven or frameshift-driven, and highly truncated alternative isoforms 527 \nfrom frameshifts are underrepresented in GENCODE annotations but not in an experimentally proteome 528 \npredicted from long-read RNA-seq data. 529 \nBiosurfer’s focus is distinct among the landscape of isoform tools, but the panoply of tools, 530 \nsteadily growing, may cause confusion as to the precise aspect of isoform biology being characterized by 531 \neach tool. Today, ma ny published tools process short read or long-read RNA seque ncing data for the 532 \npurpose of assembling full-length transcripts (Trinity)(Haas et al. 2013), discover novel transcript 533 \nvariations (Li et al. 2018), or quantifying splice events or entire isoforms (MISO, rMATs, RSEM, 534 \nKallisto)(Katz et al. 2010; Li and Dewey 2011; Shen et al. 2014; Bray et al. 2016). Other tools classify 535 \ntranscript exonic structures for novel transcripts derived from long-read RNA-seq analysis (e.g., 536 \nSQANTI, FLAIR, Bambu)(Tardaguila et al. 2018; Tang et al. 2020; Chen et al. 2023). Tools like SUPPA 537 \ndeconstruct full-length transcriptomes into individual splicing events (Alamancos et al. 2015). 538 \nCollectively, these tools analyze properties of transcripts, but with less focus on the protein 539 \neffects(Altschul et al. 1990). On the other hand, tools for comparison of proteins incorporate protein 540 \nsequence alignment (e.g., ClustalW, BLAST) (Chenna et al. 2003), but such alignments are disconnected 541 \nfrom information about the underlying genome. Indeed, there are several protein-to-genome alignment 542 \nalgorithms (e.g., Exonerate, miniprot)(Slater a nd Birney 2005; Li 2022), as well as, more recently, 543 \nmethods to align proteins with some knowledge of the underlying exonic and genomic locations of 544 \nresidues, such as the Mirage tool (Nord et al. 2018; Hanimann et al. 2022; Nord and Wheeler 2023). 545 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n26 \nRelated to the comparison of proteins, we had to designate a “reference” isoform, although the 546 \nbiological role of most isoforms are unknown (Yang et al. 2016; Reixachs /i1Solé and Eyras 2022), and 547 \nthus representative isoforms are chosen depending on the assumptions and goals of the research 548 \ncommunity (The UniProt Consortium 2017; Pozo et al. 2022). 549 \nBiosurfer analyses rely on user-defined protein isoforms. Only canonical start and stop codons are 550 \nassumed, unless non-canonical sites are annotated in a reference proteome (Mudge et al. 2022) or the 551 \nuser. Determination of the biologically relevant ORF remains an ongoing challenge.  Many ORF callers 552 \nlike transdecoder, CPAT, GMST and others predict ORFs, relying on heuristics and common features of 553 \ntranslation, which may not be the rule in every case. Currently, the prediction of proteins from deep 554 \ncoverage long-read RNA-seq datasets rely on heuris tics, such as prioritizing ORF from an alternative 555 \nisoform that shares the same start AUG codon with the reference, or selection of the most 5’ proximal 556 \nAUG (Tang et al. 2020; Miller et al. 2022), whereas ot hers have developed computationally efficient 557 \nscoring strategy that ranks more highly the ORFs with highest protein similarity to the reference 558 \n(Varabyou et al. 2022, 2023). To provide more reliable ORF annotations, experimental approaches like 559 \nRibo-Seq demarcate novel coding regions, including sites of non-canonical translation, which might be 560 \ninformation that could be incorporated in proteogenomic workflows (Mudge et al. 2022; Leblanc et al. 561 \n2024).  562 \nOur first version of Biosurfer proposes a new framework for detailed comparison of protein 563 \nisoforms, a first step towards inferring function. Further versions of this tool could map functional 564 \nelements, such as structural domains, active sites, post-translationally modified sites, or protein 565 \ninteractions, onto the altered protein regions, similar to the functionality of tappAS, isoTV, and DIGGER, 566 \nas well as other tools (de la Fuente et al. 2020; Annaldasula et al. 2021; Louadi et al. 2021). Furthermore, 567 \ngiven the link between genomic coordinates and effects on protein isoforms, Biosurfer could capture the 568 \nimpact of coding or splice-modifying genetic variants as carried through the lens of complex transcript 569 \nand protein variations, which might increase the accuracy of predicted genetic effects in ancestry- or 570 \npatient-specific populations (Rivas et al. 2015)(Cummings et al. 2017)(Yamaguchi et al. 2022)(Glinos et 571 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n27 \nal. 2022). With increasing appreciation for population scale diversity and the pan-genome, which, carried 572 \nforward, gives rise to a corresponding “pan-proteome”,  bioinformatic pipelines should be designed to 573 \nproduce automated results of all possible variations arising from newly sequenced sample. Ultimately 574 \ngenotype could be associated to the full repertoire of proteoform diversity (Aebersold et al. 2018), 575 \nespecially as proteomics approaches continues to capturing a greater swath of the protein isoform space 576 \n(Sinitcyn et al. 2023).  577 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n28 \nDATA ACCESS 578 \nAll raw data has been previously published and is described in the Methods. 579 \n 580 \nCOMPETING INTEREST STATEMENT 581 \nNo competing interests. 582 \n 583 \nACKNOWLEDGMENTS  584 \nThis work was supported by the National Library of Medicine (R01-LM014017) to G.M.S. and D.K. 585 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n29 \nREFERENCES 586 \n  587 \n  588 \nAbood A, Mesner LD, Jeffery ED, Murali M, Lehe M, Saquing J, Farber CR, Sheynkman GM. 2023. 589 \nLong-read proteogenomics to connect disease-associated sQTLs to the protein isoform effectors of 590 \ndisease. BioRxiv. 591 \nAebersold R, Agar JN, Amster IJ, Baker MS, Bertozzi CR, Boja ES, Costello CE, Cravatt BF, Fenselau 592 \nC, Garcia BA, et al. 2018. How many human proteoforms are there? Nat Chem Biol 14: 206–214. 593 \nAken BL, Ayling S, Barrell D, Clarke L, Curwen V, Fairley S, Fernandez Banet J, Billis K, García Girón 594 \nC, Hourlier T, et al. 2016. The Ensembl gene annotation system. Database (Oxford) 2016. 595 \nAlamancos GP, Pagès A, Trincado JL, Bellora N, Eyras E. 2015. Leveraging transcript quantification for 596 \nfast computation of alternative splicing profiles. RNA 21: 1521–1531. 597 \nAltschul SF, Gish W, Miller W, Myers EW, Lipman DJ. 1990. Basic local alignment search tool. J Mol 598 \nBiol 215: 403–410. 599 \nAnnaldasula S, Gajos M, Mayer A. 2021. IsoTV: processing and visualizing functional features of 600 \ntranslated transcript isoforms. Bioinformatics 37: 3070–3072. 601 \nAnvar SY, Allard G, Tseng E, Sheynkman GM, de Klerk E, Vermaat M, Yin RH, Johansson HE, 602 \nAriyurek Y, den Dunnen JT, et al. 2018. Full-length mRNA sequencing uncovers a widespread 603 \ncoupling between transcription initiation and mRNA processing. Genome Biol 19: 46. 604 \nBray NL, Pimentel H, Melsted P, Pachter L. 2016. Near-optimal probabilistic RNA-seq quantification. 605 \nNat Biotechnol 34: 525–527. 606 \nCarvill GL, Mefford HC. 2020. Poison exons in neurodevelopment and disease. Curr Opin Genet Dev 65: 607 \n98–102. 608 \nChen Y, Sim A, Wan YK, Yeo K, Lee JJX, Ling MH, Love MI, Göke J. 2023. Context-aware transcript 609 \nquantification from long-read RNA-seq data with Bambu. Nat Methods 20: 1187–1195. 610 \nChenna R, Sugawara H, Koike T, Lopez R, Gibson TJ, Higgins DG, Thompson JD. 2003. Multiple 611 \nsequence alignment with the Clustal series of programs. Nucleic Acids Res 31: 3497–3500. 612 \nClarke J, Wu H-C, Jayasinghe L, Patel A, Reid S, Bayley H. 2009. Continuous base identification for 613 \nsi\nngle-molecule nanopore DNA sequencing. Nat Nanotechnol 4: 265–270. 614 \nCooper TA, Wan L, Dreyfuss G. 2009. RNA and disease. Cell 136: 777–793. 615 \nCummings BB, Marshall JL, Tukiainen T, Lek M, Donkervoort S, Foley AR, Bolduc V, Waddell LB, 616 \nSandaradura SA, O’Grady GL, et al. 2017. Improving genetic diagnosis in Mendelian disease with 617 \ntranscriptome sequencing. Sci Transl Med 9. 618 \nde la Fuente L, Arzalluz-Luque Á, Tardáguila M, Del Risco H, Martí C, Tarazona S, Salguero P, Scott R, 619 \nLerma A, Alastrue-Agudo A, et al. 2020. tappAS: a comprehensive computational framework for the 620 \nanalysis of the functional impact of differential splicing. Genome Biol 21: 119. 621 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n30 \nDe Paoli-Iseppi R, Gleeson J, Clark MB. 2021. Isoform Age - Splice Isoform Profiling Using Long-Read 622 \nTechnologies. Front Mol Biosci 8: 711733. 623 \nde Souza VBC, Jordan BT, Tseng E, Nelson EA, Hirschi KK, Sheynkman G, Robinson MD. 2022. 624 \nTransformation of alignment files improves performance of variant callers for long-read RNA 625 \nsequencing data. BioRxiv. 626 \nEid J, Fehr A, Gray J, Luong K, Lyle J, Otto G, Peluso P, Rank D, Baybayan P, Bettman B, et al. 2009. 627 \nReal-time DNA sequencing from single polymerase molecules. Science 323: 133–138. 628 \nFiszbein A, McGurk M, Calvo-Roitberg E, Kim G, Burge CB, Pai AA. 2022. Widespread occurrence of 629 \nhybrid internal-terminal exons in human transcriptomes. Sci Adv 8: eabk1752. 630 \nFrankish A, Carbonell-Sala S, Diekhans M, Jungreis I, Loveland JE, Mudge JM, Sisu C, Wright JC, 631 \nArnan C, Barnes I, et al. 2023. GENCODE: reference annotation for the human and mouse genomes in 632 \n2023. Nucleic Acids Res 51: D942–D949. 633 \nFrankish A, Diekhans M, Jungreis I, Lagarde J, Loveland JE, Mudge JM, Sisu C, Wright JC, Armstrong 634 \nJ, Barnes I, et al. 2021. GENCODE 2021. Nucleic Acids Res 49: D916–D923. 635 \nGlinos DA, Garborcauskas G, Hoffman P, Ehsan N, Jiang L, Gokden A, Dai X, Aguet F, Brown KL, 636 \nGarimella K, et al. 2022. Transcriptome variation in human tissues revealed by long-read sequencing. 637 \nNature 608: 353–359. 638 \nGohr A, Irimia M. 2019. Matt: Unix tools for alternative splicing analysis. Bioinformatics 35: 130–132. 639 \nHaas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, 640 \nLieber M, et al. 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity 641 \nplatform for reference generation and analysis. Nat Protoc 8: 1494–1512. 642 \nHanimann J, Moch H, Zoche M, Kahraman A. 2022. IsoAligner: dynamic mapping of amino acid 643 \npositions across protein isoforms. F1000Res 11: 382. 644 \nHarrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, Aken BL, Barrell D, 645 \nZadissa A, Searle S, et al. 2012. GENCODE: the reference human genome annotation for The 646 \nENCODE Project. Genome Res 22: 1760–1774. 647 \nJoglekar A, Hu W, Zhang B, Narykov O, Diekhans M, Balacco J, Ndhlovu LC, Milner TA, Fedrigo O, 648 \nJarvis ED, et al. 2023. Single-cell long-read mRNA isoform regulation is pervasive across mammalian 649 \nbrain regions, cell types, and development. Bi oRxiv. 650 \nKatz Y, Wang ET, Airoldi EM, Burge CB. 2010. Analysis and design of RNA sequencing experiments 651 \nfor identifying isoform regulation. Nat Methods 7: 1009–1015. 652 \nKelemen O, Convertini P, Zhang Z, Wen Y, Shen M, Falaleeva M, Stamm S. 2013. Function of 653 \nalternative splicing. Gene 514: 1–30. 654 \nKozak M. 1978. How do eucaryotic ribosomes select initiation regions in messenger RNA? Cell 15: 655 \n1109–1123. 656 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n31 \nKreitzer FR, Salomonis N, Sheehan A, Huang M, Park JS, Spindler MJ, Lizarraga P, Weiss WA, So P-L, 657 \nConklin BR. 2013. A robust method to derive functional neural crest cells from human pluripotent stem 658 \ncells. Am J Stem Cells 2: 119–131. 659 \nLareau LF, Inada M, Green RE, Wengrod JC, Brenner SE. 2007. Unproductive splicing of SR genes 660 \nassociated with highly conserved and ultraconserved DNA elements. Nature 446: 926–929. 661 \nLeblanc S, Yala F, Provencher N, Lucier J-F, Levesque M, Lapointe X, Jacques J-F, Fournier I, Salzet M, 662 \nOuangraoua A, et al. 2024. OpenProt 2.0 builds a path to the functional characterization of alternative 663 \nproteins. Nucleic Acids Res 52: D522–D528. 664 \nLi B, Dewey CN. 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a 665 \nreference genome. BMC Bioinformatics 12: 323. 666 \nLi H. 2022. Protein-to-genome alignment with miniprot. arXiv [q-bioGN]. 667 \nLi YI, Knowles DA, Humphrey J, Barbeira AN, Dickinson SP, Im HK, Pritchard JK. 2018. Annotation-668 \nfree quantification of RNA splicing using LeafCutter. Nat Genet 50: 151–158. 669 \nLouadi Z, Yuan K, Gress A, Tsoy O, Kalinina OV, Baumbach J, Kacprowski T, List M. 2021. DIGGER: 670 \nexploring the functional role of alternative splicing in protein interactions. Nucleic Acids Res 49: 671 \nD309–D318. 672 \nMartelli PL, D’Antonio M, Bonizzoni P, Castrignanò T, D’Erchia AM, D’Onorio De Meo P, Fariselli P, 673 \nFinelli M, Licciulli F, Mangiulli M, et al. 2011. ASPicDB: a database of annotated transcript and 674 \nprotein variants generated by alternative splicing. Nucleic Acids Res 39: D80-5. 675 \nMehlferber MM, Jeffery ED, Saquing J, Jordan BT, Sheynkman L, Murali M, Genet G, Acharya BR, 676 \nHirschi KK, Sheynkman GM. 2022. Characterization of protein isoform diversity in human umbilical 677 \nvein endothelial cells via long-read proteogenomics. RNA Biol 19: 1228–1243. 678 \nMiller RM, Jordan BT, Mehlferber MM, Jeffery ED, Chatzipantsiou C, Kaur S, Millikin RJ, Dai Y, 679 \nTiberi S, Castaldi PJ, et al. 2022. Enhanced protein isoform characterization through long-read 680 \nproteogenomics. Genome Biol 23: 69. 681 \nModrek B, Resch A, Grasso C, Lee C. 2001. Genome-wide detection of alternative splicing in expressed 682 \nsequences of human genes. Nucleic Acids Res 29: 2850–2859. 683 \nMu\ndge JM, Ruiz-Orera J, Prensner JR, Brunet MA, Calvet F, Jungreis I, Gonzalez JM, Magrane M, 684 \nMartinez TF, Schulz JF, et al. 2022. Standardized annotation of translated open reading frames. Nat 685 \nBiotechnol 40: 994–999. 686 \nNakao M, Barrero RA, Mukai Y, Motono C, Suwa M, Nakai K. 2005. Large-scale analysis of human 687 \nalternative protein isoforms: pattern classification and correlation with subcellular localization signals. 688 \nNucleic Acids Res 33: 2355–2363. 689 \nNord A, Carey K, Hornbeck P, Wheeler T. 2018. Splice-Aware Multiple Sequence Alignment of Protein 690 \nIsoforms. ACM BCB 2018: 200–210. 691 \nNord AJ, Wheeler TJ. 2023. Mirage2’s high-quality spliced protein-to-genome mappings produce 692 \naccurate multiple-sequence alignments of isoforms. PLoS ONE 18: e0285225. 693 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n32 \nO’Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, Rajput B, Robbertse B, Smith-694 \nWhite B, Ako-Adjei D, et al. 2016. Reference sequence (RefSeq) database at NCBI: current status, 695 \ntaxonomic expansion, and functional annotation. Nucleic Acids Res 44: D733-45. 696 \nPardo-Palacios F, Reese F, Carbonell-Sala S, Diekhans M, Liang C, Wang D, Williams B, Adams M, 697 \nBehera A, Lagarde J, et al. 2021. Systematic assessment of long-read RNA-seq methods for transcript 698 \nidentification and quantification. Res Sq. 699 \nPozo F, Rodriguez JM, Martínez Gómez L, Vázquez J, Tress ML. 2022. APPRIS principal isoforms and 700 \nMANE Select transcripts define reference splice variants. Bioinformatics 38: ii89–ii94. 701 \nReese F, Williams B, Balderrama-Gutierrez G, Wyman D, Çelik MH, Rebboah E, Rezaie N, Trout D, 702 \nRazavi-Mohseni M, Jiang Y, et al. 2023. The ENCODE4 long-read RNA-seq collection reveals distinct 703 \nclasses of transcript structure diversity. BioRxiv. 704 \nReixachs-Solé M, Eyras E. 2022. Uncovering the impacts of alternative splicing on the proteome with 705 \ncurrent omics techniques. Wiley Interdiscip Rev RNA 13: e1707. 706 \nReyes A, Huber W. 2018. Alternative start and termination sites of transcription drive most transcript 707 \nisoform differences across human tissues. Nucleic Acids Res 46: 582–592. 708 \nRivas MA, Pirinen M, Conrad DF, Lek M, Tsang EK, Karczewski KJ, Maller JB, Kukurba KR, DeLuca 709 \nDS, Fromer M, et al. 2015. Human genomics. Effect of predicted protein-truncating genetic variants on 710 \nthe human transcriptome. Science 348: 666–669. 711 \nRodriguez JM, Maietta P, Ezkurdia I, Pietrelli A, Wesselink J-J, Lopez G, Valencia A, Tress ML. 2013. 712 \nAPPRIS: annotation of principal and alternative splice isoforms. Nucleic Acids Res 41: D110-7. 713 \nSharon D, Tilgner H, Grubert F, Snyder M. 2013. A single-molecule long-read survey of the human 714 \ntranscriptome. Nat Biotechnol 31: 1009–1014. 715 \nShen S, Park JW, Lu Z, Lin L, Henry MD, Wu YN, Zhou Q, Xing Y. 2014. rMATS: robust and flexible 716 \ndetection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci USA 717 \n111: E5593-601. 718 \nSinitcyn P, Richards AL, Weatheritt RJ, Brademan DR, Marx H, Shishkova E, Meyer JG, Hebert AS, 719 \nWestphall MS, Blencowe BJ, et al. 2023. Global detection of human variants and isoforms by deep 720 \nproteome sequencing. Nat Biotechnol 41: 1776–1786. 721 \nS\nlater GSC, Birney E. 2005. Automated generation of heuristics for biological sequence comparison. 722 \nBMC Bioinformatics 6: 31. 723 \nTanabe M, Sasai N, Nagata K, Liu XD, Liu PC, Thiele DJ, Nakai A. 1999. The mammalian HSF4 gene 724 \ngenerates both an activator and a repressor of heat shock genes by alternative splicing. J Biol Chem 725 \n274: 27845–27856. 726 \nTang AD, Soulette CM, van Baren MJ, Hart K, Hrabeta-Robinson E, Wu CJ, Brooks AN. 2020. Full-727 \nlength transcript characterization of SF3B1 mutation in chronic lymphocytic leukemia reveals 728 \ndownregulation of retained introns. Nat Commun 11: 1438. 729 \nTapial J, Ha KCH, Sterne-Weiler T, Gohr A, Braunschweig U, Hermoso-Pulido A, Quesnel-Vallières M, 730 \nPermanyer J, Sodaei R, Marquez Y, et al. 2017. An atlas of alternative splicing profiles and functional 731 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint \n\n  \n \n33 \nassociations reveals new regulatory programs and genes that simultaneously express multiple major 732 \nisoforms. Genome Res 27: 1759–1768. 733 \nTardaguila M, de la Fuente L, Marti C, Pereira C, Pardo-Palacios FJ, Del Risco H, Ferrell M, Mellado M, 734 \nMacchietto M, Verheggen K, et al. 2018. SQANTI: extensive characterization of long-read transcript 735 \nsequences for quality control in full-length transcriptome identification and quantification. Genome Res 736 \n28: 396–411. 737 \nTian L, Jabbari JS, Thijssen R, Gouil Q, Amarasinghe SL, Voogd O, Kariyawasam H, Du MRM, 738 \nSchuster J, Wang C, et al. 2021. Comprehensive characterization of single-cell full-length isoforms in 739 \nhuman and mouse with long-read sequencing. Genome Biol 22: 310. 740 \nTranchevent L-C, Aubé F, Dulaurier L, Benoit-Pilven C, Rey A, Poret A, Chautard E, Mortada H, 741 \nDesmet F-O, Chakrama FZ, et al. 2017. Identification of protein features encoded by alternative exons 742 \nusing Exon Ontology. Genome Res 27: 1087–1097. 743 \nThe UniProt Consortium. 2017. UniProt: the universal protein knowledgebase. Nucleic Acids Res 45: 744 \nD158–D169. 745 \nVarabyou A, Erdogdu B, Salzberg SL, Pertea M. 2023. Investigating Open Reading Frames in Known 746 \nand Novel Transcripts using ORFanage. BioRxiv. 747 \nVarabyou A, Sommer MJ, Erdogdu B, Shinder I, Minkin I, Chao K-H, Park S, Heinz J, Pockrandt C, 748 \nShumate A, et al. 2022. CHESS 3: an improved, comprehensive catalog of human genes and transcripts 749 \nbased on large-scale expression data, phylogenetic analysis, and protein structure. BioRxiv. 750 \nVeiga DFT, Nesta A, Zhao Y, Deslattes Mays A, Huynh R, Rossi R, Wu T-C, Palucka K, Anczukow O, 751 \nBeck CR, et al. 2022. A comprehensive long-read isoform analysis platform and sequencing resource 752 \nfor breast cancer. Sci Adv 8: eabg6711. 753 \nWang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, Kingsmore SF, Schroth GP, Burge CB. 754 \n2008. Alternative isoform regulation in human tissue transcriptomes. Nature 456: 470–476. 755 \nWang L, Park HJ, Dasari S, Wang S, Kocher J-P, Li W. 2013. CPAT: Coding-Potential Assessment Tool 756 \nusing an alignment-free logistic regression model. Nucleic Acids Res 41: e74. 757 \nWorkman RE, Tang AD, Tang PS, Jain M, Tyson JR, Razaghi R, Zuzarte PC, Gilpatrick T, Payne A, 758 \nQuick J, et al. 2019. Nanopore native RNA sequencing of a human poly(A) transcriptome. Nat 759 \nMethods 16: 1297–1305. 760 \nY\namaguchi K, Ishigaki K, Suzuki A, Tsuchida Y, Tsuchiya H, Sumitomo S, Nagafuchi Y, Miya F, 761 \nTsunoda T, Shoda H, et al. 2022. Splicing QTL analysis focusing on coding sequences reveals 762 \nmechanisms for disease susceptibility loci. Nat Commun 13: 4659. 763 \nYang X, Coulombe-Huntington J, Kang S, Sheynkman GM, Hao T, Richardson A, Sun S, Yang F, Shen 764 \nYA, Murray RR, et al. 2016. Widespread expansion of protein interaction capabilities by alternative 765 \nsplicing. Cell 164: 805–817.  766 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}