Biosurfer for systematic tracking of regulatory mechanisms leading to protein isoform diversity

25 Long-read RNA sequencing has shed light on transcriptomic complexity, but questions remain about the 26 functionality of downstream protein products. We introduce Biosurfer, a computational approach for 27 comparing protein isoforms, while systematically tracking the transcriptional, splicing, and translational 28 variations that underlie differences in the sequences of the protein products. Using Biosurfer, we analyzed 29 the differences in 32,799 pairs of GENCODE annotated protein isoforms, finding a majority (70%) of 30 variable N-termini are due to the alternative transcription start sites, while only 9% arise from 5’ UTR 31 alternative splicing. Biosurfer’s detailed tracking of nucleotide-to-residue relationships helped reveal an 32 uncommonly tracked source of single amino acid residue changes arising from the codon splits at 33 junctions. For 17% of internal sequence changes, such split codon patterns lead to single residue 34 differences, termed “ragged codons”. Of variable C-termini, 72% involve splice- or intron retention-35 induced reading frameshifts. We found an unusual pattern of reading frame changes, in which the first 36 frameshift is closely followed by a distinct second frameshift that restores the original frame, which we 37 term a “snapback” frameshift. We analyzed long read RNA-seq-predicted proteome of a human cell line 38 and found similar trends as compared to our GENCODE analysis, with the exception of a higher 39 proportion of isoforms predicted to undergo nonsense-mediated decay. Biosurfer’s comprehensive 40 characterization of long-read RNA-seq datasets should accelerate insights of the functional role of protein 41 isoforms, providing mechanistic explanation of the origins of the proteomic diversity driven by the 42 alternative splicing. Biosurfer is available as a Python package at https://github.com/sheynkman-43 lab/biosurfer. 44 45

LRS Special Issue, Alternative splicing, long-read sequencing, protein isoforms, GENCODE, 46 protein sequence, alternative transcriptional start site (altTSS), reading frame shift, intron retention, 47 nonsense mediated decay (NMD), open reading frame (ORF), transcriptional start site (TSS), 48 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 3 transcriptional termination site (TSS), poison exon, Transcription Initiation Site (TIS), sequence 49 alignment 50 51

52 Through the isoform diversifying mechanisms of alternative transcription, splicing, and 53 polyadenylation, nearly every human gene can produce multiple protein products, with ~20K genes 54 giving rise to at least 180K annotated isoforms (Frankish et al. 2023). The pathway from gene to protein 55 is marked by several regulatory mechanisms that are highly tuned across development and cell states, 56 with disruption of this regulation producing aberrant isoforms that lead to pathophysiological states such 57 as cancer and cardiovascular disease (Cooper et al. 2009). Hence, approaches are needed to systematically 58 characterize the upstream regulatory causes and functional impacts of such protein isoform sequence 59 changes. 60 Transcript and, by extension, protein isoform diversity may now be globally characterized at 61 great depth for individual samples (Glinos et al. 2022; Reese et al. 2023). Transcript diversity can be 62 readily characterized by long read RNA-Seq, which employs single molecule sequencing of individual 63 cDNA or RNA molecules to determine the sequence across the entire length of spliced transcripts (Sharon 64 et al. 2013; Workman et al. 2019; Pardo/i1Palacios et al. 2021; Tian et al. 2021; Joglekar et al. 2023), with 65 platforms from Oxford Nanopore and PacBio being most commonly used (Clarke et al. 2009; Eid et al. 66 2009). Long read RNA sequencing captures long range connectivity between multiple exons of a 67 transcript and can reveal complex splice patterns unattainable by short read sequencing (De Paoli /i1Iseppi 68 et al. 2021), including dependencies across distal splicing events (Anvar et al. 2018) and alternative 5’ 69 and 3’ transcript usage. Given this readily characteri zed complexity afforded by long read sequencing, a 70 natural question is the extent to which such variations of the transcriptome lead to functional effects of the 71 proteome. Towards this goal, a necessary step is defining the potential proteome. Both our group and 72 others have reported methods in which long-read-derived transcript sequences serve as templates for 73 predicting full-length protein isoform sequences, thereby providing a global snapshot of potential protein 74 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 4 isoforms expressed in a particular biological cond ition (Miller et al. 2022; Veiga et al. 2022; Abood et al. 75 2023). 76 The complexity of alternative splicing (AS)—which for practical purposes in this manuscript we 77 define here as all transcriptional variations, including alternative transcription start sites (TSS) and 78 transcription termination sites (TTS)—can be observed and characterized at different levels: RNA 79 transcript, open reading frames (ORFs), and finally protein sequences (Reixachs /i1Solé and Eyras 2022). 80 The complex interplay between the changes occurring to AS variants and their ORF and protein products 81 are hard to characterize and quantify. Changes in mRNA sequence may lead to non-linear or traditionally 82 untracked variations. For example, a subtle splicing event could lead to a reading frame shift and thus to 83 more dramatic changes to the C-terminus of the protein than the originating small change at the RNA 84 level would suggest. Or, AS could occur at codon boundaries, leading to altered amino acid (AA) 85 identities of codons that technically overlap in genome-space but are differentially “split” across exon-86 exon junctions. And complex interplay may also be observed between transcriptional variations and ORF 87 choice, as alternative 5’ transcription or AS could lead to differentially availability of initiator codons, 88 delimiting start codon choice co-translationally. 89 We argue that rather than being an esoteric exercise, the ability to characterize all potential 90 interplay of RNA-protein variation is critical for fully elucidating the transcript and proteomic diversity 91 encoded within long-read RNA-seq datasets. As long-read RNA-seq approaches are increasingly adopted 92 in large-scale studies of hundreds of samples (Glinos et al. 2022; Reese et al. 2023) and are maturing into 93 stable tools being adopted by the community (Pardo /i1Palacios et al. 2021), extracting all sources of 94 biological molecular diversity is critical. Such interplay cannot be characterized by comparison of 95 isoforms using just one modality, such as transcript-focused annotation tools like SQANTI or Matt 96 (Tardaguila et al. 2018; Gohr and Irimia 2019), or conventional protein sequence alignment tools like 97 ClustalW (Chenna et al. 2003) or BLAST (Altschul et al. 1990). R ecently, multi-modal comparisons have 98 been reported. For example, ORFanage is an approach for large-scale annotation of ORFs across 99 predicted transcripts in the CHESS database the main focus being optimal selection of ORFs based on 100 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 5 protein-alignments (Varabyou et al. 2023). Other tools geared towards the phylogenetics community have 101 developed frame-aware alignment in which the alignment scoring system is includes penalties for frame-102 shift-inducing gaps (Evans and Loose 2015; Ranwez et al. 2011; Jammali et al. 2022) . However, these 103 tools do not comprehensively elucidate the interplay between transcript and protein variation. 104 To characterize how AS impacts protein sequence, several bioinformatic tools and databases have 105 been developed, such as VastDB, ASPicDB, ExonOntology, and DIGGER (Martelli et al. 2011; Tapial et 106 al. 2017; Tranchevent et al. 2017; Louadi et al. 2021). Tools such as tappAS and IsoTV annotate how 107 protein isoform sequence and potential functional differences (de la Fuente et al. 2020; Annaldasula et al. 108 2021). For example, tappAS is a Java application for quantifying differential isoform usage but also to 109 functionally annotate such isoforms, using the module IsoAnnot (de la Fuente et al. 2020). IsoAnnot maps 110 protein features (e.g., Pfam domains) across isoforms of a gene, and determines how splicing leads to 111 partial or full removal of protein features, indicating potential changes to molecular functions. Despite 112 existence of these tools, of need is the ability to systematically capture all possible effects to protein 113 isoforms, with the accompanying information of the underlying complex RNA-protein relationships. 114 We developed Biosurfer, a computational pipeline that tracks simultaneously the changes at all 115 three levels, to understand the impact of alternative splicing on transcriptome, ORFeome, and proteome 116 diversity. Biosurfer computes details not immediately apparent from genome annotation files or manual 117 inspection in genome browsers, such as how between isoforms of the same gene the frame of translation 118 and codon topology influences amino acid sequence identity changes. In order to accomplish this multi-119 layered comparison, a genome is used as a “scaffold ” to exactly position all nucleotides, codons, AA 120 elements, and associate with each element local and context-dependent attributes. The resulting data 121 structure of three distinct yet interlinked layers inform on the upstream biological mechanisms leading to 122 protein sequence changes. To demonstrate the utility of Biosurfer, we globally characterize variation 123 across protein isoforms in the reference human annotation (GENCODE) and protein isoforms predicted 124 from long-read RNA-seq of a human stem cell. This characterization includes comprehensive elaboration 125 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 6 of all sources of N-terminal, internal, and C-terminal protein variations observed, highlighting 126 mechanisms of RNA-protein interplay. 127 128

129 Biosurfer package for isoform analysis and visualization 130 Biosurfer is a computational pipeline that performs a multi-layered comparison between a pair of 131 isoforms, in which differences at three different levels: RNA (nucleotides, nt, ORF (codon), and protein 132 (AA residues, AA) are simultaneously tracked. The developed data structure enables not only the 133 comparison of AAs, but tracking frameshifts, patterns of codon splitting at junctions, and the attendant 134 upstream nucleotide and codon differences that explain AA changes. Such tracking aids in the systematic 135 annotation of explanatory mechanism(s) underlying AA residue changes, such as whether a substitution 136 of a stretch of AA residues is due to alternative splicing or a frameshift (Supplementary Figure S1). 137 The Biosurfer pipeline is organized into the following three stages ( Figure 1). First, an SQLite 138 database is populated with detailed information on each isoform at the transcript, ORF, and protein 139 product levels. For each isoform, the required inputs are (i) a transcript FASTA, (ii) a protein FASTA, 140 and (iii) a matching GTF with both Exon and CDS features. The inputs can either be extracted from the 141

annotations (e.g., GENCODE (Harrow et al. 2012)), or they can be user-defined ( e.g., predicted 142 protein isoform sequences from long-read RNA-seq data (Miller et al. 2022)by using ORF callers such as 143 CPAT (Wang et al. 2013) or Transdecoder (Haas et al. 2013)). Second, mu ltilayered isoform-level 144 alignments are generated, and each alignment is represented as three key data structures: t-blocks, c-145 blocks, and p-blocks corresponding to the view of the aligned isoforms at the transcript, codon sequence, 146 and protein sequence levels, respectively ( Supplementary Figure S2 ). Third, all information is 147 summarized in tabular format, to facilitate analysis of the transcription-level and codon-level mechanisms 148 driving the proteomic diversity. In addition, a visual representation of such mechanisms integrated with 149 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 7 the protein-level view of the isoforms can be output as png files.150 151 Figure 1: Biosurfer for analysis and visualization of protein isoform sequence differences. 152 Biosurfer analyzes protein isoforms from reference annotations (e.g., GENCODE) or proteins predicted 153 from long-read RNA-seq data. Analysis initiates with the creation of an SQLite database is populated 154 with isoform-relevant elements. Biosurfer performs a multi-layered comparison of transcript-, codon-, and 155 protein-level differences between pairs of protein isoforms. Variable regions as well as their associated 156 annotations are output in tabular format and visualization files, which includes protein-relevant details 157 such as the frame of translation. Note that the terms “Match”, “Deletion”, etc. represent very different 158 comparisons depending on the biological layer. GTF=Gene Transfer Format; iPSCs=induced pluripotent 159 stem cells. 160 161 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 8 Data structures for presenting transcript-, ORF-, and protein-level isoform alignments in Biosurfer 162 Biosurfer compares transcript-, ORF-, and protein-level sequences for each gene. One isoform is 163 selected as reference and the other isoforms are denoted as alternative. When aligning each alternative 164 isoform to the reference, the matched regions are referred to as matched blocks. The remaining regions on 165 each sequence that are not matched blocks are referred to as unmatched blocks. See Supplementary File 166 1 for step-by-step schema of the comparison process, which is summarized below. 167 Transcript blocks (t-blocks). T-blocks represent subsegments of the transcript sequence that are shared or 168 unique to the reference or alternative isoforms. Transcript-level differences are determined by analyzing 169 the transcript-to-genome coordinates alignment (GTF file) provided as an input by the user. Specifically, 170 Biosurfer defines the aligned exonic regions that are shared or unique to each isoform. The resulting 171 ranges, called t-blocks, are categorized as Match, Deletion, or Insertion t-blocks. Deletion or Insertion t-172 blocks are further annotated with the associated biological mechanism leading to the transcript nucleotide 173 change, e.g., alternative transcriptional start site, splicing event, or polyadenylation. The alternative 174 splicing events are then further classified into four basic types: retained intron, alternative donor, 175 alternative exon, or alternative acceptor. 176 Codon blocks (c-blocks). C-blocks represent a codon-centric layer defined through the comparison of the 177 protein coding regions of transcripts, i.e., open reading frames (ORF), between two isoforms. The c-block 178 data structure is the most complex, but critical, layer in Biosurfer that connects information between the 179 transcript and protein layers. 180 For ultimate granularity and precision, ORFs are compared based on the alignment of codons that 181 overlap in the genome space, in which one codon in the reference isoform is compared with another 182 codon in the alternative isoform (see Supplemental Methods for additional details). First, codons across 183 the two isoforms are “ paired” based on their mutual positions and base overlap: in a basic scenario, the 184 two codons are identical and their positions match ( Table1, Supplemental Figure S3A ); in a more 185 complex scenario, the codons are split and only partially overlapped (1 or 2 bases, see Table1, 186 Supplementary Fig S3B and C). In all other cases, the codons will be designated as “unpaired” ( Table1, 187 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 9 Supplementary Figure S3B). To keep track of unpaired codons within the data structure, they are linked188 to a “placeholder” codon, which serves to maintain consistency of the c-block structures. Once aligned,189 each single codon pair is categorized with respect to multiple attributes that explain their relationship.190 Overall, the paired and unpaired codons are classified into 9 categories based on their translation status,191 frameshift status, and codon topology (Table 1). 192 193 Table 1: Categories of overlapping codon pairs that are the basis for codon blocks (c-blocks). 194 195 Protein blocks (p-blocks). P-blocks represent a protein-centric layer defined through t- and c- block196 guided comparisons of two protein isoforms, thus focusing solely on the relationships between the AA197 residues of the respective protein isoforms. This decision was motivated by the fact that in most cases,198 functional differences exhibited by an alternative protein isoform are attributable to differences in protein199 rather than nucleotide sequences between the alternative and reference isoforms , and these differences200 could be conceptually decoupled from the knowledge of the upstream (ORF- or transcript-level) residue -201 altering mechanisms. 202 ed d, ip. us, ck A es, in es - .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 10 The p-blocks for a pair of isoforms are determined by translating c-blocks defined at the previous 203 level. The resulting data structure represents subsegments of contiguous stretches of AA residues from 204 each protein sequence. Each subsegment could be as short as a small translated portion of exon (even a 205 single residue, in the case of NAGNAG splicing), but it can also span multiple exons of a gene. Because 206 the actual matching happened at the two previous levels (t-blocks and c-blocks), no alignment is required 207 at the p-block level. The resulting protein subsegments can be either (1) fully matched (100% sequence 208 identity, across all residues in the subsegment), or (2) mismatched, where a subsegment in one protein 209 sequence will be aligned against a gapped region in the other sequence. Thus, a p-block represents either 210 a matched pair of protein subsequences or a subsequence matched against a gapped region. P-blocks are 211 then classified as Match, Insertion, Deletion, or Substitution changes. Substitution p-blocks must arise 212 from a combination of insertion/deletion/frameshift events found at the c-block level. 213 Overall, these p-block changes are agnostic to the upstream mechanisms, e.g., at the transcript or 214 ORF levels; however, at the same time, the corresponding upstream mechanisms can be retrieved and 215 analyzed within Biosurfer, unlike using traditional protein aligners. 216 217 Analysis of GENCODE isoforms using Biosurfer 218 The principal and alternative isoforms were defined using APPRIS annotation of genes in GENCODE 219 (Rodriguez et al. 2013). First, the set of APPRIS isoforms is identified for each gene from the input 220 genome data (GENCODE v42, basic annotation (Frankish et al. 2021)) by extracting the key transcript 221 features, such as 'transcript_id', 'transcript_name', and the associated APPRIS ‘tag’ (Rodriguez et al. 222 2013). Second, the transcript's APPRIS status is determined as 'principal', 'alternative', or ‘none’, based on 223 first rank-ordering of transcripts based on APPRIS tag and setting the transcript with the highest APPRIS 224 value as ‘principle’ and all other transcripts as ‘alternative’. We did not consider further genes with only 225 one annotated isoform and lacking an APPRIS tag (‘none’). Subsequently, we compiled a structured 226 dataset for each transcript, encompassing identifiers, gene associations, strand orientation, and APPRIS 227 status. During the annotation, the strand orientation is accounted for when necessary. 228 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 11 229 Long read RNA-seq analysis of WTC-11 cell line 230 Total RNA from WTC-11 cells was extracted using the RNeasy Kit (Qiagen) and analyzed on an 231 Agilent Bioanalyzer. We observed a RNA concentration of 30 ng/uL with the RNA Integrity Score (RIN) 232 of 9.9. As described previously, (Mehlferber et al. 2022) cDNA was synthesized from the extracted RNA 233 and the Iso-Seq Express Kit SMRT Bell Express Template prep kit 2.0 (Pacific Biosciences) was used on 234 a Sequel II system to obtain long-read sequence information and output Circular Consensus (CCS) reads. 235 Data is available at the Sequence Read Archive: SRR18130587 and previously published (de Souza et al. 236 2022). 237 We analyzed the WTC-11 data with a proteogenomics Nextflow pipeline we previously 238 developed (Miller et al. 2022) (Mehlferber et al. 2022). The output CCS reads from long-read sequencing 239 were processed with Iso-Seq3 and SQANTI3 (version 1.3) for transcript isoform classification and quality 240 assessment. The CPAT(Wang et al. 2013) algorithm was used to predict Open Reading Frames (ORFs), 241 which were grouped into protein isoforms. 242 243 Biosurfer is implemented as a Python package freely available at GitHub repository: 244 https://github.com/sheynkman-lab/biosurfer. The analysis code is available at 245 https://github.com/sheynkman-lab/biosurfer_analysis. All necessary input files and intermediate and final 246 output files from Biosurfer analysis are uploaded to Zenodo at: https://zenodo.org/records/10822882. 247 248

249 Characterization of altered protein regions in the human annotation (GENCODE) 250 Here, we demonstrate the utility of Biosurfer through a genome-wide analysis of protein isoforms 251 in the GENCODE annotation (basic annotation, version 42 (Frankish et al. 2021)). We analyzed 35,083 252 reference-alternative protein isoform pairs across 11,815 genes ( Figure 2A ). Each pair consists of the 253

protein isoform for a gene—the APPRIS principal isoform (Rodriguez et al. 2013)—and an 254 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 12 alternative protein isoform. The number of isoform pairs correspond to the number of “alternative” 255 isoforms for a gene ( Supplementary Table S1). Genes lacking an APPRIS annotation (8,166 of 19,981) 256 were excluded from the analysis (Supplementary Table S2). 257 Globally, we found a total of 44,326 altered protein regions with an average of 1.3 altered regions 258 per isoform. Note that we are using the term altered protein region interchangeably with p-block (see 259 Methods). Altered protein regions are contiguous regions of altered protein sequence relative to the 260 “reference” (i.e., APPRIS principal) protein, which includes p-block Insertions, Deletions, and 261 Substitutions. A majority (79%, 27,672 of 35,083) of isoform pairs contain a single altered region (Figure 262 2B). Still, 7,411 alternative isoforms contain two or more discontinuous altered regions. Notably, some 263 isoform pairs exhibited an extreme number of regions. For example, 14 regions are found for proteins of 264 DNAH14 (Reference: DNAH14-220, Alternative: DNAH14-211), which is explained by its extremely 265 large number of exons (86 exons for DNAH14-220). 266 Among the 44,326 altered protein regions ( Figure 2C ), the median number of affected amino 267 acids (lost or gained due to a protein insertion, deletion, or substitution) is 49 AA, with the first and third 268 quartiles containing 21 AA and 128 AA, respectively. Among the altered regions, 14% (6,189) are 269 insertions, 47% (20,780) are deletions, and 39% (17,357) are substitutions (Figure 2D ). Full annotations 270 for these altered regions at the protein and codon-level, representing the Biosurfer output for protein-271 blocks (p-blocks) and codon-blocks (c-blocks) can be found in Supplementary Table S3, S4. 272 The lengths of p-block insertions tend to be shorter than the length of deletions (p < 2.2e-16, 273 Mann-Whitney U test) or substitutions (p < 2.2e-16, Mann-Whitney U test) ( Figure 2E-H ). Since the 274 deletion or insertion status of a protein region is dependent on which isoform is denoted as the reference, 275 this trend may reflect the tendency that longer isofor ms are more likely to be defined as the “reference” 276 isoform, which may be biologically driven or influenced by genome annotation guidelines (Rodriguez et 277 al. 2013; O’Leary et al. 2016; The UniProt Consortium 2017; Frankish et al. 2021; Varabyou et al. 2023). 278 We found a similar trend for p-block substitutions; the lengths for substituted regions tended to be longer 279 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 13 for the sequence affected in the reference ( Figure 2G ) versus the alternative ( Figure 2H ) isoform;280 although, for 3,263 (18% of cases), the affected regions in alternative isoforms can be longer (Figure 2I). 281 282 Figure 2: Characterization of altered protein regions (Biosurfer p-blocks) across the GENCODE283 annotated human proteome. ( A) Schematic of genes with one or two alternative protein isoforms and284 distribution of the number of alternative protein isoforms per gene. (B) Schematic of the altered protein285 regions (here, highlighted in pink and b lue), displayed relative to the underlying transcript structures. The286 proteome-wide distribution of the number of affected protein regions observed per alternative isoform.287 (C) Distribution of the length of altered protein regions across the annotated proteome. Differences 288 greater than 600 AAs are not included (2,347 cases, 5.3% of the data). These altered protein regions289 include cases of 1) deleted regions, 2) inserted regions, and 3) the region (in the reference isoform) in290 which one polypeptide subsegme nt is substituted for another. In other words, this distribution (C)291 represents an aggregation of the distributions shown in panels E- G. (D) Fraction of altered protein regions292 m; E nd in he m. ns in C) ns .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 14 affected by insertions, deletions, and substitutions. (E-H) Distribution of the lengths of altered protein 293 regions for (E) deletions, (F) insertions, (G) substituted region in the reference isoform, and (H) 294 substituted region in the alternative isoform. (I) Comparison of the lengths of altered regions in the 295

versus alternative isoforms for substitutions. 296 297 298 Analysis of N-terminal protein variations 299 Variations in the N-termini are found in 28% (12,504 of 44,326) of all possible reference-300 alternative isoform pairs, corresponding to 5,872 genes (Supplementary Table S5). We examined for the 301 N-terminal variations the explanatory mechanism, including alternative transcription start sites (TSSs), 302 alternative splicing, and alternative translation initiation sites (TISs). 303 All cases of variable N-termini involve two initiation codons (AUGs), one upstream and one 304 downstream, relative to the genome. A major category we first observed are those in which the N-305 terminus is different due to start codons that are mutually exclusively present across the two transcript 306 isoforms (Figure 3A ). Specifically, the start codon present in the reference isoform is absent from the 307 alternative isoform, and vice versa. These “mutually exclusive start codons” or MXS were observed for 308 3,123 reference-alternative isoform pairs ( Figure 3A). MXS may arise either from an alternative TSS or 309 from alternative splicing in the 5’ UTR. Strikingly, we found that nearly all (99%, 3,097 of 3,123) cases 310 are caused by alternative TSS usage ( Figure 3A, hatched region of the bar), with only a small, but non-311 zero, fraction (1%, 26 of 3,123) of MXS cases arise from splicing of the 5’ UTR, in which splicing 312 regulation is influencing N-terminal usage. An example of TSS-driven MXS for PRKACA is shown in 313 Figure 3B (pair, PRKACA-201 and PRKACA-202). 314 The second category is when the upstream start codon is transcribed in only one of the two 315 isoforms, but the downstream start codon is present in both transcripts. We refer to this scenario as shared 316 downstream starts (SDS), of which there were 6,878 cases ( Figure 3A). Like with cases of MXS, SDS 317 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 15 arises primarily from alternative TSS usage (80% of cases) versus 5’ UTR alternative splicing (20% of 318 cases). 319 MXS and SDS are common patterns underlying alterations of the N-termini or protein, driven by 320 differential availability of initiator codons in the mature transcript. We asked if there may be differences 321 in length of such N-terminal alterations for SDS versus MXS events. Measuring the differential length of 322 the affected N-terminal regions between reference-alternative isoform pairs, we found that, on average, 323 SDS tends to affect a greater proportion of protein length, as compared to MXS ( Figure 3C , 324 Supplementary Figure S4; p = 2.8e-148, Mann-Whitney U test). The larger differences in length driven 325 by SDS could be explained by cases in which transcription is initiated from internal sites of the gene, 326 giving rise to an ORF that corresponds to a subsequence of the ORF in the other isoform, theoretically 327 producing a truncated C-terminal-containing subsequence of the full-length protein. 328 Recently, a mechanism related to SDS was described in which internal exons (not the 5’ most 329 exon, first exon, of a transcript) in one transcript can be immediately downstream of a DNA element of 330 novel promoter activity and thus serve as the first transcribed exon in other isoforms (Fiszbein et al. 331 2022). Such so-called hybrid exons thus can operate as both sites of transcription initiation and alternative 332 splicing, in effect, swapping their roles depending on the regulatory context. Of the cases of SDS, we 333 observed that ~22% correspond to these hybrid exon swaps ( Figure 3D, Supplementary Table S6). The 334 functional consequences of hybrid exon usage are not well understood; however, one potential function 335 could be the production of protein with a truncated N-terminus, which could remove signal peptide 336 sequences or binding domains (Kelemen et al. 2013). 337 In addition to MXS and SDS, wherein start codon availability is controlled through differential 338 transcription, we also observed many cases in which both upstream and downstream start codons co-339 occur in one or both transcript isoforms of a pair. In these cases, the choice of start codon may be 340 influenced by co-translational regulation, e.g., ribosome initiates translation at alternative initiation sites 341 (altTIS). 342 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 16 For 2,054 cases, we found altTIS, which we classified as instances of a mutually shared start 343 (MSS) codon. We also found 449 cases in which the upstream start codon is present in both isoforms, but 344 the downstream start codon is only present in one isoform. 345 In such cases of shared upstream start (SUS) ( Figure 3A ), the explanatory co-translational 346 mechanism is not as clear, based on the ribosomal scanning model of translation initiation (Kozak 1978). 347 The upstream start codon present in both isoforms would need to be bypassed by the ribosome only in the 348 alternative isoform. Therefore, some of the SUS annotations may need to be validated or may be 349 erroneous ORF calls, as early ORF prediction workflows attribute higher scores to longer ORFs (Wang et 350 al. 2013; Varabyou et al. 2023), under-annotating ORFs that utilize the upstream (annotated) start codon 351 but is much shorter than the reference due to a reading frame shift (Wang et al. 2013). 352 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 17 353 354 Figure 3: Analysis of mechanisms underlying variable N- terminal proteins across the GENCODE355 annotated human proteome. (A) Distribution of the types of alternative N- terminal regions, classified356 based on presence and translational status of the start codon. Hatches denote the fraction of alternative N -357 terminal regions associated with alternative transcription start sites, as opposed to 5’ UTR alternative358 splicing. (B) Biosurfer output of altered N-terminal regions for PRKACA gene that undergo MXS (light359 green arrows) and SDS (dark green arrows). In the example of MXS, the yellow bars above and below 360 E ed - ve ht .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 18 PRKACA-202 (alternate) transcript indicate the N-terminal ranges that differ between the reference and 361 alternative isoform. The yellow bar above PRKACA-202 shows the N-terminal protein sequence that is 362 specific to the reference (PRKACA-201) and the bar below the transcript indicates the N-terminal region 363 specific to the alternative isoform. The Biosurfer bars span the intronic lines between exons, but intronic 364 regions do not contribute to the protein sequence differences. In the example of SDS, the pink Biosurfer 365 bar above the transcript of PRKACA-203 represents the range of transcript sequence that is translated in 366 the reference, but not translated in the alternative isoform. (C) Scatterplot of the length of affected N-367 terminal variation in the reference versus alternative, faceted by mutually exclusive starts (MXS) or 368 shared downstream start (SDS) status. An interesting case of an SDS leading to unique N-terminal 369 sequence in the reference is caused by usage of a different frame at the initiation of translation, with an 370 example shown for isoforms of the gene FHOD3. (D) Fraction of shared downstream starts (SDS) caused 371 by hybrid exon swaps. 372 373 374 Analysis of internal protein variations 375 The internal regions of the protein isoforms account for 43% (19,263 of 44,326) of possible 376 reference-alternative isoform pairs, corresponding to 6,673 genes ( Supplementary Table S7 ). A large 377 majority of these regions (80%, 15,444 of 19,263) are caused by single simple splicing events: exon 378 skipping, alternative acceptor, alternative donor, or an in-frame retained intron. As expected, exon 379 skipping events are most numerous making up 9,505 of cases. In terms of the general effect on protein 380 sequence, most altered regions (70%, 13,453 of 19,263) lead to a deletion or removal of AA residues 381 (Figure 4A), and, again, in such cases, exon skipping is most common (51%, 6,908 of 13,453). 382 Going beyond simple splicing events, we observed that 19% of variable internal regions (3,701 of 383 19,263) were associated with multiple events, which we refer to as compound, or linked, splicing events. 384 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 19 We found that 56% (2,099 of 3,701) compound events involve multi-exon skipping, the rest being 385 combinations of alternative donor/acceptor sites with exon skipping/inclusion (Figure 4B). 386 Complex nucleotide to amino acid relationships that affect internal protein sequence 387 Previous studies of the impact splicing on proteins have typically focused on cases in which 388 differential splicing of transcript regions directly corresponds to changes in protein sequence 389 (Reixachs/i1Solé and Eyras 2022). However, in many instances, there is not a simple one-to-one 390 relationship between nucleotides in a transcript and the corresponding amino acid identities in the 391 encoded protein. Such complexity alters the protein sequence in non-intuitive ways. Using the detailed 392 codon tracking afforded by Biosurfer, we systematically characterized the protein-level impact of 393 variations not commonly described: codons that span junctions and unusual reading frame shifts. 394 To characterize differentially split codons, we examined all paired codons (see c-block section in 395 Methods) that are split across junctions and determined the identity of the associated AAs 396 (Supplemental Figure S3 ). Across all internal altered protein regions, we found 17% (3,213 of 19,263) 397 of regions that are flanked by one or more split codon pairs that encode different amino acid residues 398 (Figure 4C , also see Table 1 ). These so-called ragged codons affect a single residue and are always 399 adjacent to an altered protein region. Interestingly, while split codons are not particularly enriched by 400 splice event type (Chi-square test: p-value = 5.54e-99), we found that ragged codons are more frequently 401 found in protein insertions compared to deletions or substitutions (Chi-square test: p-value = 4.11e-106). 402 For a majority of splice-driven frame shifts, the shifted frame is maintained to the end of the 403 protein, leading to a protein variation affecting the C-terminus. Such frameshifted proteins could lead to 404 truncated proteins or destabilization of the transcript via mechanisms such as nonsense mediated decay 405 (NMD). Surprisingly, we found an un commonly characterized pattern of successive readin g frame shifts 406 that exclusively affects the internal residues of a protein. In these cases, the alternative isoform’s reading 407 frame is shifted due to one splicing event, but then shifts back into register of the reference frame due to a 408 second, independent splicing event. We refer to these events as “snapback” frameshifts, as there is a 409 return back to the original reading frame. Snapback frameshifts have been previously observed, such as in 410 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 20 the gene HSF4 , which produces internally frameshifted isoforms that have been demonstrated to exert411 different regulatory effects (Tanabe et al. 1999), but the snapback phenomenon generally speaking has not412 been systematically described. Within GENCODE, we found 118 examples of suc h snapback isoforms413 across 95 genes, including FAS , and PLEKHJ1 (Figure 4D, Supplementary Table S8). What is notable414 about these cases is that the same underlying genomic sequence encodes different amino acid residues,415 and genetic mutations could lead to two different residue changes depending on the isoform. 416 417 Figure 4: Analysis of internal protein altered regions across proteins in GENCODE. (A) Frequencies418 of the categories of the splicing mechanism underlyi ng internal protein sequence changes, split by their419 protein-coding impact (Deletion, Insertion, Substitution). All sequence regions involve a reference -420 alternative isoform pair. (B) Frequency of compound splicing events across the altered internal regions.421 (C) Proportion of each internal protein region type for which there exists a split codon pattern near its 422 boundaries that would cause a single amino acid difference, or “ragged” codon. (D) Examples of423 successive frameshifting (snapback” frameshift) that leads to an affected protein region that is wh olly424 internal to the protein, for genes FAS and PLEKHJ1. 425 426 Analysis of C-terminal variations 427 ert ot ms le es, ies eir - ns. of lly .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 21 C-terminal changes make up 28% (12,241 of 44,326) of reference-alternative pairs, 428 corresponding to 6,138 genes (Supplementary Table S9). 429 To break down the sources of C-termini variability, we found it useful to distinguish the most 430 direct preceding cause of altered C-termini. In principle, all C-terminal changes must arise from an 431 upstream splicing event that influences the termination codon used (notwithstanding post-translational 432 cleavage events). However, such changes could be further classified. The altered C-terminus could arise 433 from alternative terminal (i.e., last) exons, each harboring a different stop codon, so that the splicing event 434 more or less directly influences the stop codon availability (i.e., direct splice-driven events). In other 435 instances, C-terminal changes could arise from a somewhat indirect relationship to the splice event, such 436 as when a splicing event causes a translational frameshift, in effect, “revealing” in the other frame a new 437 stop codon that is now decoded by the ribosome (i.e., frameshift-driven events) ( Supplementary Table 438 S9). 439 In direct splice-driven events, the stop codon availability is dictated by the actively transcribed 440 regions that contain the stop codon. We find this scenario for 72% (8,877 of 12,241) of all C-terminal 441 variations (Figure 5A). These variations can be further classified based on the pattern of splicing at the C 442 terminus. The first pattern involves an exon extension into the intron region, introducing a premature stop 443 codon. These exon extensions introduce termination, or “EXIT”, make up 3,498 (39.4% of 8,877) cases 444 (Figure 5B). The second pattern involves usage of alternative terminal coding exons or “ATE”, making 445 up 3,301 (37.1% of 8,877) of cases ( Figure 5B ). Overall, EXIT and ATE changes lead to a shorter C-446 termini in the alternative isoform (distribution shown in Figure 5C). 447 EXIT versus ATE reflect how different “modes” of spliceosome regulation could lead to distinct 448 C-terminal consequences. In EXIT, the reference-containing donor splice site fails to be spliced in the 449 alternative isoform, leading to partial or full intron retention. An example of EXIT is shown in the right 450 panel of Figure 5C for the pair, TDRD12-206 and TDRD102-202. In ATE, on the other hand, the 451 spliceosome catalyzes splicing at one of two splice site acceptor sites, influencing terminal exon identity 452 and thus stop codon used. A well-known pattern of ATE are poison exons, a mechanism by which 453 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 22 inclusion of an exon leads to a premature termination codon that either elicits nonsense mediated decay or 454 generates a truncated protein product (Carvill and Mefford 2020). Poison exons are evolutionarily 455 conserved and likely play a role in downregulating gene expression (Lareau et al. 2007). We found 550 456 (6.2% of total 8,877) cases of potential poison exons. Other ATE patterns include one in which the 457 alternative last exon of the alternative isoform resides in the UTR region of the reference isoform, 458 suggesting that such sequences in the 3’ UTR could dually code both transcript and protein functional 459 elements. We found 819 (9% of 8,877) cases of such alternative last exon in UTR (ALE in UTR) ( Figure 460 5B). A third ATE variation, referred to as a cut-out splice terminal exon (COSTE), the 5’ end of the last 461 exon is shared, but the alternative isoform utilizes a splice site that skips over the remaining portion of the 462 last exon in the reference, thereby creating a different last exon not found in the original reference. We 463 identified 273 cases of this pattern (3% of 8,877) (Figure 5B and Supplementary Figure S5). 464 Frameshift-driven events influence stop codon usage somewhat indirectly through shifts in the 465 translational reading frame. In such cases, a splice-induced reading frame shift causes all downstream 466 codons to be read in a different frame and stop codons are “revealed”, or decoded, by the ribosome. We 467 found 3,364 (28% of 12,241) cases of frameshift-induced C-terminal changes ( Figure 5A, Examples are 468 shown in Figure 5D , pairs, GIPC3-201 and GIPC3-202, along with FANCM-201 and FANCM-231). 469 Generally speaking, frameshifts lead to a dramatic shortening of the C-terminal region in the alternative 470 isoform; however, we found 549 cases (16% of 3,364) in which the C-terminal region is longer in the 471 alternative versus the reference isoform (Figure D). 472 We also observed across all frameshift-driven events a depletion of isoform pairs in which a large 473 portion of the reference isoform (e.g., 2,000 AA or longer) is truncated due to a frameshift in the 474 alternative isoform (Figure D), a trend not observed in an experimentally predicted proteome (see section 475 below and Supplementary Figure S9B ), likely representing gene annotation decisions, as dramatically 476 truncating frameshift events would lead to predicted NMD and filtered out or reassigned an NMD biotype 477 (Harrow et al. 2012). 478 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 23 479 Figure 5: Analysis of alternative C-terminal protein sequences in GENCODE v42. (A) Frequency of480 alternative C-terminal categories based on splicing or frameshifts being the primary driving factor. (B) 481 Distribution of the frequencies for various splice-driven patterns. (C) Scatterplot of the length of splice -482 driven C-terminal variation in the reference versus alternative. An example of this category is observed in483 the TDRD12 gene. TDRD12 undergoes splice- driven alteration causing an alternative terminal exon in484 TDRD12-202 to harbor the stop codon (D) Scatterplot of the length of frameshift-driven C- terminal485 variation in the reference versus alternative. Biosurfer plot examples of frameshift- driven category486 illustrated in GIPC43 and FANCM genes. 487 488 489 Characterization of altered protein regions across a long-read predicted proteome 490 The process of defining the reference proteome heavily draws from sources of experimental491 evidence such as deeply sequenced cell and tissue types, and the protein isoform sequences represent a n492 aggregate model of the human proteome. Therefore, to characterize potential isoform related protein493 of - in in al ry tal n in .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 24 variability in a specific biological condition, we employed a “long read proteogenomics” pipeline 494 (Kreitzer et al. 2013), generating a proteome predicted from long-read RNA-seq transcript sequences 495 collected from a karyotypically normal human stem cell line (WTC-11). Using Biosurfer, we 496 characterized the landscape of protein isoforms in the WTC-11 proteome and found 44,962 protein 497 isoform pairs across 10,144 genes ( Supplementary Table S10 , p-block and c-block outputs in 498 Supplementary Table S11 and S12 ). Assigning the highest expressed transcript as the “reference”, we 499 defined 53,915 altered protein regions ( Supplementary Figure S6 ). Compared to the GENCODE 500 analysis, similar trends were observed for N-terminal ( Supplementary Figure S7 ) and internal region 501 variation (Supplementary Figure S8 , snapback isoforms in Supplementary Table S13 ). Besides these 502 similar trends, the experimental proteome returned a higher number of C-terminal variations that involved 503 intron retention events (EXIT event type, see Figure 5B) as compared to GENCODE isoforms, matching 504 earlier findings from EST and cDNA data (Nakao et al. 2005)(Modrek et al. 2001)( Supplementary 505 Figure S9). 506 507

508 To study the functional impact of alternatively spliced protein isoforms, it is critical to track 509 precise differences in protein isoform sequences and link such variations to the upstream explanatory 510 mechanisms. However, it is challenging to systematically characterize the full interplay between genomic 511 and proteomic variations, which hinders discoveries of novel biological variations represented in a long-512 read RNA-seq dataset. We developed Biosurfer, a computational approach, available as a Python 513 package, that systematically extracts protein isoform sequence variations while maintaining the explicit 514 links to their underlying transcriptional and post-transcriptional mechanisms. 515 To demonstrate the utility of Biosurfer, we characterized protein isoform differences across an 516 annotated (GENCODE) and long-read RNA-seq predicted proteome. Using Biosurfer’s interlinked 517 transcript, codon, and protein data structures, we determined the upstream mechanisms explaining 518 isoform alterations, uncovering surprising complexity. First, we confirmed past observations of 519 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 25 alternative transcription underlying most N-terminal variations (Reyes and Huber 2018), and most 520 internal protein sequence differences arising from single splicing events, such as exon skipping (Wang et 521 al. 2008). However, our study goes beyond these known trends. Upon detailed tracking of codons 522 flanking these altered protein regions, we found distinct split codon patterns that change the encoded 523 amino acid residue identity and thus contribute to variation of AA residues. We also found an unusual 524 frameshift pattern that involves successive reading frame shifts that leads to a change in protein regions 525 that is entirely internal to the protein, referred to here as "snapback" frameshifts. And last, C-terminal 526 differences are primarily splice-driven or frameshift-driven, and highly truncated alternative isoforms 527 from frameshifts are underrepresented in GENCODE annotations but not in an experimentally proteome 528 predicted from long-read RNA-seq data. 529 Biosurfer’s focus is distinct among the landscape of isoform tools, but the panoply of tools, 530 steadily growing, may cause confusion as to the precise aspect of isoform biology being characterized by 531 each tool. Today, ma ny published tools process short read or long-read RNA seque ncing data for the 532 purpose of assembling full-length transcripts (Trinity)(Haas et al. 2013), discover novel transcript 533 variations (Li et al. 2018), or quantifying splice events or entire isoforms (MISO, rMATs, RSEM, 534 Kallisto)(Katz et al. 2010; Li and Dewey 2011; Shen et al. 2014; Bray et al. 2016). Other tools classify 535 transcript exonic structures for novel transcripts derived from long-read RNA-seq analysis (e.g., 536 SQANTI, FLAIR, Bambu)(Tardaguila et al. 2018; Tang et al. 2020; Chen et al. 2023). Tools like SUPPA 537 deconstruct full-length transcriptomes into individual splicing events (Alamancos et al. 2015). 538 Collectively, these tools analyze properties of transcripts, but with less focus on the protein 539 effects(Altschul et al. 1990). On the other hand, tools for comparison of proteins incorporate protein 540 sequence alignment (e.g., ClustalW, BLAST) (Chenna et al. 2003), but such alignments are disconnected 541 from information about the underlying genome. Indeed, there are several protein-to-genome alignment 542 algorithms (e.g., Exonerate, miniprot)(Slater a nd Birney 2005; Li 2022), as well as, more recently, 543

to align proteins with some knowledge of the underlying exonic and genomic locations of 544 residues, such as the Mirage tool (Nord et al. 2018; Hanimann et al. 2022; Nord and Wheeler 2023). 545 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 26 Related to the comparison of proteins, we had to designate a “reference” isoform, although the 546 biological role of most isoforms are unknown (Yang et al. 2016; Reixachs /i1Solé and Eyras 2022), and 547 thus representative isoforms are chosen depending on the assumptions and goals of the research 548 community (The UniProt Consortium 2017; Pozo et al. 2022). 549 Biosurfer analyses rely on user-defined protein isoforms. Only canonical start and stop codons are 550 assumed, unless non-canonical sites are annotated in a reference proteome (Mudge et al. 2022) or the 551 user. Determination of the biologically relevant ORF remains an ongoing challenge. Many ORF callers 552 like transdecoder, CPAT, GMST and others predict ORFs, relying on heuristics and common features of 553 translation, which may not be the rule in every case. Currently, the prediction of proteins from deep 554 coverage long-read RNA-seq datasets rely on heuris tics, such as prioritizing ORF from an alternative 555 isoform that shares the same start AUG codon with the reference, or selection of the most 5’ proximal 556 AUG (Tang et al. 2020; Miller et al. 2022), whereas ot hers have developed computationally efficient 557 scoring strategy that ranks more highly the ORFs with highest protein similarity to the reference 558 (Varabyou et al. 2022, 2023). To provide more reliable ORF annotations, experimental approaches like 559 Ribo-Seq demarcate novel coding regions, including sites of non-canonical translation, which might be 560 information that could be incorporated in proteogenomic workflows (Mudge et al. 2022; Leblanc et al. 561 2024). 562 Our first version of Biosurfer proposes a new framework for detailed comparison of protein 563 isoforms, a first step towards inferring function. Further versions of this tool could map functional 564 elements, such as structural domains, active sites, post-translationally modified sites, or protein 565 interactions, onto the altered protein regions, similar to the functionality of tappAS, isoTV, and DIGGER, 566 as well as other tools (de la Fuente et al. 2020; Annaldasula et al. 2021; Louadi et al. 2021). Furthermore, 567 given the link between genomic coordinates and effects on protein isoforms, Biosurfer could capture the 568 impact of coding or splice-modifying genetic variants as carried through the lens of complex transcript 569 and protein variations, which might increase the accuracy of predicted genetic effects in ancestry- or 570 patient-specific populations (Rivas et al. 2015)(Cummings et al. 2017)(Yamaguchi et al. 2022)(Glinos et 571 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 27 al. 2022). With increasing appreciation for population scale diversity and the pan-genome, which, carried 572 forward, gives rise to a corresponding “pan-proteome”, bioinformatic pipelines should be designed to 573 produce automated results of all possible variations arising from newly sequenced sample. Ultimately 574 genotype could be associated to the full repertoire of proteoform diversity (Aebersold et al. 2018), 575 especially as proteomics approaches continues to capturing a greater swath of the protein isoform space 576 (Sinitcyn et al. 2023). 577 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 28 DATA ACCESS 578 All raw data has been previously published and is described in the Methods. 579 580 COMPETING INTEREST STATEMENT 581 No competing interests. 582 583 ACKNOWLEDGMENTS 584 This work was supported by the National Library of Medicine (R01-LM014017) to G.M.S. and D.K. 585 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 29

586 587 588 Abood A, Mesner LD, Jeffery ED, Murali M, Lehe M, Saquing J, Farber CR, Sheynkman GM. 2023. 589 Long-read proteogenomics to connect disease-associated sQTLs to the protein isoform effectors of 590 disease. BioRxiv. 591 Aebersold R, Agar JN, Amster IJ, Baker MS, Bertozzi CR, Boja ES, Costello CE, Cravatt BF, Fenselau 592 C, Garcia BA, et al. 2018. How many human proteoforms are there? Nat Chem Biol 14: 206–214. 593 Aken BL, Ayling S, Barrell D, Clarke L, Curwen V, Fairley S, Fernandez Banet J, Billis K, García Girón 594 C, Hourlier T, et al. 2016. The Ensembl gene annotation system. Database (Oxford) 2016. 595 Alamancos GP, Pagès A, Trincado JL, Bellora N, Eyras E. 2015. Leveraging transcript quantification for 596 fast computation of alternative splicing profiles. RNA 21: 1521–1531. 597 Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. 1990. Basic local alignment search tool. J Mol 598 Biol 215: 403–410. 599 Annaldasula S, Gajos M, Mayer A. 2021. IsoTV: processing and visualizing functional features of 600 translated transcript isoforms. Bioinformatics 37: 3070–3072. 601 Anvar SY, Allard G, Tseng E, Sheynkman GM, de Klerk E, Vermaat M, Yin RH, Johansson HE, 602 Ariyurek Y, den Dunnen JT, et al. 2018. Full-length mRNA sequencing uncovers a widespread 603 coupling between transcription initiation and mRNA processing. Genome Biol 19: 46. 604 Bray NL, Pimentel H, Melsted P, Pachter L. 2016. Near-optimal probabilistic RNA-seq quantification. 605 Nat Biotechnol 34: 525–527. 606 Carvill GL, Mefford HC. 2020. Poison exons in neurodevelopment and disease. Curr Opin Genet Dev 65: 607 98–102. 608 Chen Y, Sim A, Wan YK, Yeo K, Lee JJX, Ling MH, Love MI, Göke J. 2023. Context-aware transcript 609 quantification from long-read RNA-seq data with Bambu. Nat Methods 20: 1187–1195. 610 Chenna R, Sugawara H, Koike T, Lopez R, Gibson TJ, Higgins DG, Thompson JD. 2003. Multiple 611 sequence alignment with the Clustal series of programs. Nucleic Acids Res 31: 3497–3500. 612 Clarke J, Wu H-C, Jayasinghe L, Patel A, Reid S, Bayley H. 2009. Continuous base identification for 613 si ngle-molecule nanopore DNA sequencing. Nat Nanotechnol 4: 265–270. 614 Cooper TA, Wan L, Dreyfuss G. 2009. RNA and disease. Cell 136: 777–793. 615 Cummings BB, Marshall JL, Tukiainen T, Lek M, Donkervoort S, Foley AR, Bolduc V, Waddell LB, 616 Sandaradura SA, O’Grady GL, et al. 2017. Improving genetic diagnosis in Mendelian disease with 617 transcriptome sequencing. Sci Transl Med 9. 618 de la Fuente L, Arzalluz-Luque Á, Tardáguila M, Del Risco H, Martí C, Tarazona S, Salguero P, Scott R, 619 Lerma A, Alastrue-Agudo A, et al. 2020. tappAS: a comprehensive computational framework for the 620 analysis of the functional impact of differential splicing. Genome Biol 21: 119. 621 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 30 De Paoli-Iseppi R, Gleeson J, Clark MB. 2021. Isoform Age - Splice Isoform Profiling Using Long-Read 622 Technologies. Front Mol Biosci 8: 711733. 623 de Souza VBC, Jordan BT, Tseng E, Nelson EA, Hirschi KK, Sheynkman G, Robinson MD. 2022. 624 Transformation of alignment files improves performance of variant callers for long-read RNA 625 sequencing data. BioRxiv. 626 Eid J, Fehr A, Gray J, Luong K, Lyle J, Otto G, Peluso P, Rank D, Baybayan P, Bettman B, et al. 2009. 627 Real-time DNA sequencing from single polymerase molecules. Science 323: 133–138. 628 Fiszbein A, McGurk M, Calvo-Roitberg E, Kim G, Burge CB, Pai AA. 2022. Widespread occurrence of 629 hybrid internal-terminal exons in human transcriptomes. Sci Adv 8: eabk1752. 630 Frankish A, Carbonell-Sala S, Diekhans M, Jungreis I, Loveland JE, Mudge JM, Sisu C, Wright JC, 631 Arnan C, Barnes I, et al. 2023. GENCODE: reference annotation for the human and mouse genomes in 632 2023. Nucleic Acids Res 51: D942–D949. 633 Frankish A, Diekhans M, Jungreis I, Lagarde J, Loveland JE, Mudge JM, Sisu C, Wright JC, Armstrong 634 J, Barnes I, et al. 2021. GENCODE 2021. Nucleic Acids Res 49: D916–D923. 635 Glinos DA, Garborcauskas G, Hoffman P, Ehsan N, Jiang L, Gokden A, Dai X, Aguet F, Brown KL, 636 Garimella K, et al. 2022. Transcriptome variation in human tissues revealed by long-read sequencing. 637 Nature 608: 353–359. 638 Gohr A, Irimia M. 2019. Matt: Unix tools for alternative splicing analysis. Bioinformatics 35: 130–132. 639 Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, 640 Lieber M, et al. 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity 641 platform for reference generation and analysis. Nat Protoc 8: 1494–1512. 642 Hanimann J, Moch H, Zoche M, Kahraman A. 2022. IsoAligner: dynamic mapping of amino acid 643 positions across protein isoforms. F1000Res 11: 382. 644 Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, Aken BL, Barrell D, 645 Zadissa A, Searle S, et al. 2012. GENCODE: the reference human genome annotation for The 646 ENCODE Project. Genome Res 22: 1760–1774. 647 Joglekar A, Hu W, Zhang B, Narykov O, Diekhans M, Balacco J, Ndhlovu LC, Milner TA, Fedrigo O, 648 Jarvis ED, et al. 2023. Single-cell long-read mRNA isoform regulation is pervasive across mammalian 649 brain regions, cell types, and development. Bi oRxiv. 650 Katz Y, Wang ET, Airoldi EM, Burge CB. 2010. Analysis and design of RNA sequencing experiments 651 for identifying isoform regulation. Nat Methods 7: 1009–1015. 652 Kelemen O, Convertini P, Zhang Z, Wen Y, Shen M, Falaleeva M, Stamm S. 2013. Function of 653 alternative splicing. Gene 514: 1–30. 654 Kozak M. 1978. How do eucaryotic ribosomes select initiation regions in messenger RNA? Cell 15: 655 1109–1123. 656 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 31 Kreitzer FR, Salomonis N, Sheehan A, Huang M, Park JS, Spindler MJ, Lizarraga P, Weiss WA, So P-L, 657 Conklin BR. 2013. A robust method to derive functional neural crest cells from human pluripotent stem 658 cells. Am J Stem Cells 2: 119–131. 659 Lareau LF, Inada M, Green RE, Wengrod JC, Brenner SE. 2007. Unproductive splicing of SR genes 660 associated with highly conserved and ultraconserved DNA elements. Nature 446: 926–929. 661 Leblanc S, Yala F, Provencher N, Lucier J-F, Levesque M, Lapointe X, Jacques J-F, Fournier I, Salzet M, 662 Ouangraoua A, et al. 2024. OpenProt 2.0 builds a path to the functional characterization of alternative 663 proteins. Nucleic Acids Res 52: D522–D528. 664 Li B, Dewey CN. 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a 665

genome. BMC Bioinformatics 12: 323. 666 Li H. 2022. Protein-to-genome alignment with miniprot. arXiv [q-bioGN]. 667 Li YI, Knowles DA, Humphrey J, Barbeira AN, Dickinson SP, Im HK, Pritchard JK. 2018. Annotation-668 free quantification of RNA splicing using LeafCutter. Nat Genet 50: 151–158. 669 Louadi Z, Yuan K, Gress A, Tsoy O, Kalinina OV, Baumbach J, Kacprowski T, List M. 2021. DIGGER: 670 exploring the functional role of alternative splicing in protein interactions. Nucleic Acids Res 49: 671 D309–D318. 672 Martelli PL, D’Antonio M, Bonizzoni P, Castrignanò T, D’Erchia AM, D’Onorio De Meo P, Fariselli P, 673 Finelli M, Licciulli F, Mangiulli M, et al. 2011. ASPicDB: a database of annotated transcript and 674 protein variants generated by alternative splicing. Nucleic Acids Res 39: D80-5. 675 Mehlferber MM, Jeffery ED, Saquing J, Jordan BT, Sheynkman L, Murali M, Genet G, Acharya BR, 676 Hirschi KK, Sheynkman GM. 2022. Characterization of protein isoform diversity in human umbilical 677 vein endothelial cells via long-read proteogenomics. RNA Biol 19: 1228–1243. 678 Miller RM, Jordan BT, Mehlferber MM, Jeffery ED, Chatzipantsiou C, Kaur S, Millikin RJ, Dai Y, 679 Tiberi S, Castaldi PJ, et al. 2022. Enhanced protein isoform characterization through long-read 680 proteogenomics. Genome Biol 23: 69. 681 Modrek B, Resch A, Grasso C, Lee C. 2001. Genome-wide detection of alternative splicing in expressed 682 sequences of human genes. Nucleic Acids Res 29: 2850–2859. 683 Mu dge JM, Ruiz-Orera J, Prensner JR, Brunet MA, Calvet F, Jungreis I, Gonzalez JM, Magrane M, 684 Martinez TF, Schulz JF, et al. 2022. Standardized annotation of translated open reading frames. Nat 685 Biotechnol 40: 994–999. 686 Nakao M, Barrero RA, Mukai Y, Motono C, Suwa M, Nakai K. 2005. Large-scale analysis of human 687 alternative protein isoforms: pattern classification and correlation with subcellular localization signals. 688 Nucleic Acids Res 33: 2355–2363. 689 Nord A, Carey K, Hornbeck P, Wheeler T. 2018. Splice-Aware Multiple Sequence Alignment of Protein 690 Isoforms. ACM BCB 2018: 200–210. 691 Nord AJ, Wheeler TJ. 2023. Mirage2’s high-quality spliced protein-to-genome mappings produce 692 accurate multiple-sequence alignments of isoforms. PLoS ONE 18: e0285225. 693 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 32 O’Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, Rajput B, Robbertse B, Smith-694 White B, Ako-Adjei D, et al. 2016. Reference sequence (RefSeq) database at NCBI: current status, 695 taxonomic expansion, and functional annotation. Nucleic Acids Res 44: D733-45. 696 Pardo-Palacios F, Reese F, Carbonell-Sala S, Diekhans M, Liang C, Wang D, Williams B, Adams M, 697 Behera A, Lagarde J, et al. 2021. Systematic assessment of long-read RNA-seq methods for transcript 698 identification and quantification. Res Sq. 699 Pozo F, Rodriguez JM, Martínez Gómez L, Vázquez J, Tress ML. 2022. APPRIS principal isoforms and 700 MANE Select transcripts define reference splice variants. Bioinformatics 38: ii89–ii94. 701 Reese F, Williams B, Balderrama-Gutierrez G, Wyman D, Çelik MH, Rebboah E, Rezaie N, Trout D, 702 Razavi-Mohseni M, Jiang Y, et al. 2023. The ENCODE4 long-read RNA-seq collection reveals distinct 703 classes of transcript structure diversity. BioRxiv. 704 Reixachs-Solé M, Eyras E. 2022. Uncovering the impacts of alternative splicing on the proteome with 705 current omics techniques. Wiley Interdiscip Rev RNA 13: e1707. 706 Reyes A, Huber W. 2018. Alternative start and termination sites of transcription drive most transcript 707 isoform differences across human tissues. Nucleic Acids Res 46: 582–592. 708 Rivas MA, Pirinen M, Conrad DF, Lek M, Tsang EK, Karczewski KJ, Maller JB, Kukurba KR, DeLuca 709 DS, Fromer M, et al. 2015. Human genomics. Effect of predicted protein-truncating genetic variants on 710 the human transcriptome. Science 348: 666–669. 711 Rodriguez JM, Maietta P, Ezkurdia I, Pietrelli A, Wesselink J-J, Lopez G, Valencia A, Tress ML. 2013. 712 APPRIS: annotation of principal and alternative splice isoforms. Nucleic Acids Res 41: D110-7. 713 Sharon D, Tilgner H, Grubert F, Snyder M. 2013. A single-molecule long-read survey of the human 714 transcriptome. Nat Biotechnol 31: 1009–1014. 715 Shen S, Park JW, Lu Z, Lin L, Henry MD, Wu YN, Zhou Q, Xing Y. 2014. rMATS: robust and flexible 716 detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci USA 717 111: E5593-601. 718 Sinitcyn P, Richards AL, Weatheritt RJ, Brademan DR, Marx H, Shishkova E, Meyer JG, Hebert AS, 719 Westphall MS, Blencowe BJ, et al. 2023. Global detection of human variants and isoforms by deep 720 proteome sequencing. Nat Biotechnol 41: 1776–1786. 721 S later GSC, Birney E. 2005. Automated generation of heuristics for biological sequence comparison. 722 BMC Bioinformatics 6: 31. 723 Tanabe M, Sasai N, Nagata K, Liu XD, Liu PC, Thiele DJ, Nakai A. 1999. The mammalian HSF4 gene 724 generates both an activator and a repressor of heat shock genes by alternative splicing. J Biol Chem 725 274: 27845–27856. 726 Tang AD, Soulette CM, van Baren MJ, Hart K, Hrabeta-Robinson E, Wu CJ, Brooks AN. 2020. Full-727 length transcript characterization of SF3B1 mutation in chronic lymphocytic leukemia reveals 728 downregulation of retained introns. Nat Commun 11: 1438. 729 Tapial J, Ha KCH, Sterne-Weiler T, Gohr A, Braunschweig U, Hermoso-Pulido A, Quesnel-Vallières M, 730 Permanyer J, Sodaei R, Marquez Y, et al. 2017. An atlas of alternative splicing profiles and functional 731 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint 33 associations reveals new regulatory programs and genes that simultaneously express multiple major 732 isoforms. Genome Res 27: 1759–1768. 733 Tardaguila M, de la Fuente L, Marti C, Pereira C, Pardo-Palacios FJ, Del Risco H, Ferrell M, Mellado M, 734 Macchietto M, Verheggen K, et al. 2018. SQANTI: extensive characterization of long-read transcript 735 sequences for quality control in full-length transcriptome identification and quantification. Genome Res 736 28: 396–411. 737 Tian L, Jabbari JS, Thijssen R, Gouil Q, Amarasinghe SL, Voogd O, Kariyawasam H, Du MRM, 738 Schuster J, Wang C, et al. 2021. Comprehensive characterization of single-cell full-length isoforms in 739 human and mouse with long-read sequencing. Genome Biol 22: 310. 740 Tranchevent L-C, Aubé F, Dulaurier L, Benoit-Pilven C, Rey A, Poret A, Chautard E, Mortada H, 741 Desmet F-O, Chakrama FZ, et al. 2017. Identification of protein features encoded by alternative exons 742 using Exon Ontology. Genome Res 27: 1087–1097. 743 The UniProt Consortium. 2017. UniProt: the universal protein knowledgebase. Nucleic Acids Res 45: 744 D158–D169. 745 Varabyou A, Erdogdu B, Salzberg SL, Pertea M. 2023. Investigating Open Reading Frames in Known 746 and Novel Transcripts using ORFanage. BioRxiv. 747 Varabyou A, Sommer MJ, Erdogdu B, Shinder I, Minkin I, Chao K-H, Park S, Heinz J, Pockrandt C, 748 Shumate A, et al. 2022. CHESS 3: an improved, comprehensive catalog of human genes and transcripts 749 based on large-scale expression data, phylogenetic analysis, and protein structure. BioRxiv. 750 Veiga DFT, Nesta A, Zhao Y, Deslattes Mays A, Huynh R, Rossi R, Wu T-C, Palucka K, Anczukow O, 751 Beck CR, et al. 2022. A comprehensive long-read isoform analysis platform and sequencing resource 752 for breast cancer. Sci Adv 8: eabg6711. 753 Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, Kingsmore SF, Schroth GP, Burge CB. 754 2008. Alternative isoform regulation in human tissue transcriptomes. Nature 456: 470–476. 755 Wang L, Park HJ, Dasari S, Wang S, Kocher J-P, Li W. 2013. CPAT: Coding-Potential Assessment Tool 756 using an alignment-free logistic regression model. Nucleic Acids Res 41: e74. 757 Workman RE, Tang AD, Tang PS, Jain M, Tyson JR, Razaghi R, Zuzarte PC, Gilpatrick T, Payne A, 758 Quick J, et al. 2019. Nanopore native RNA sequencing of a human poly(A) transcriptome. Nat 759

16: 1297–1305. 760 Y amaguchi K, Ishigaki K, Suzuki A, Tsuchida Y, Tsuchiya H, Sumitomo S, Nagafuchi Y, Miya F, 761 Tsunoda T, Shoda H, et al. 2022. Splicing QTL analysis focusing on coding sequences reveals 762 mechanisms for disease susceptibility loci. Nat Commun 13: 4659. 763 Yang X, Coulombe-Huntington J, Kang S, Sheynkman GM, Hao T, Richardson A, Sun S, Yang F, Shen 764 YA, Murray RR, et al. 2016. Widespread expansion of protein interaction capabilities by alternative 765 splicing. Cell 164: 805–817. 766 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted March 18, 2024. ; https://doi.org/10.1101/2024.03.15.585320doi: bioRxiv preprint