JUMPg: an Integrative Proteogenomics Pipeline Identifying Unannotated Proteins in Human Brain and Cancer Cells

Module 1: Customized Database Construction

JUMPg provides three options of constructing customized databases (Figure 2), including mutations, splice junctions, and non-canonical translation from all six frames.

Mutation peptides are directly input by genomic variant files (e.g. VCF files) or indirectly derived from RNA-seq raw data (FASTAQ file). RNA reads are aligned to human reference genome (hg19) by STAR⁴¹ with UCSC gene annotations.⁴² Mutations are detected by GATK following the “GATK Best Practices” protocol^{43, 44}, stored in VCF format, and classified by AnnoVar (version 23, August 2013)⁴⁵ into different functional categories (e.g. non-silent, intronic, intergenic, etc.). Amino acid (AA) sequences that flank (±30 AAs) non-silent mutations are extracted according to UCSC gene annotations into FASTA format.

Splice junction database is generated with a similar strategy described by Sheynkman et al. 2013.²⁶ Splice junction peptides are input by junction files from the STAR program,⁴¹ and filtered with following criteria: i) at least 2 uniquely mapped reads, and ii) containing canonical splice sites (GT/AG, GC/AG, or AT/AC). The two genomic sequence fragments that flank a pair of splice donor and acceptor sites are retrieved from reference genome, and translated in three open reading frames (ORFs). If the ORF contains a stop codon, the sequence is trimmed. If the stop codon locates upstream of the splice site, the sequence is discarded because no novel amino acid sequence is produced.

Non-canonical translated peptides are made from RNA-seq reads through all 6 frames. RNAseq reads are first assembled into transcripts using trinity (version 2.1.1) ⁴⁶ with default parameters (Supporting Information Notes). Each transcript is then translated in 6 frames, split by stop codons and only ORFs with at least 66 AAs (i.e., ~200 bp) are retained.

Module 2: Database Search

Our tag-based database search engine JUMP ⁴⁰ was modified to adapt multistage strategy and maintain regular search function (version 13.1.0). In the regular search, JUMP performs preprocessing, tag generation, MS/MS pattern matching and scoring as previously reported,⁴⁰ and outputs two temporary files: a DTAS file containing mass-corrected and intensity-normalized MS/MS peak information, and a TAGS file containing de novo tags for each MS/MS scan. We evaluated MS/MS quality by the two files and kept only unmatched high quality spectra for multistage analysis (see Module 5). In addition, we also used these files to bypass two time-consuming steps (preprocessing and tag generation) during the multistage analysis, which is termed tag recycling. Common parameters include ±6 ppm for precursor ion mass tolerance, dynamic mass shift for oxidized Met (+15.9949 Da), and the consideration of a, b, and y ions. Fixed modifications on the N termini and Lys (+229.1629 Da) are also included for the TMT dataset.

Module 3: Peptide-spectrum Match Filtering

The peptide-spectrum matches (PSMs) are filtered as previously reported.⁴⁷ The target-decoy strategy is used to evaluate false discovery rate.^{48, 49} PSMs are first filtered by user-specified parameters (e.g. minimal peptide length and minimum search score), then by precursor ion mass accuracy. The resulting PSMs are further grouped by precursor ion charge state, tryptic ends, and then filtered by matching scores (Jscore and ΔJscore) to achieve the user-specified level of false discovery rate.

Module 4: Peptide Annotation and Visualization

Accepted peptides from mutation or junction search space are examined respectively for the specific amino acid(s) mutated or that span the splicing junction. The survived peptides are condensed into events (e.g., by collapsing different modification forms), with genomic positions stored in BED format for visualization in the UCSC genome browser.⁵⁰

For peptides identified from RNA 6FT database, peptides from known search space are first removed and the corresponding RNA transcripts are aligned to reference genome (hg19) by BLAT (v36, default parameters).⁵¹ Top-scored alignment result for each transcript is retained from which the genomic position of each peptide is derived. Peptides are further classified according to UCSC gene annotations. In case one peptide can be assigned to multiple isoforms, the following priority will be applied: CDS > UTR > intron > non-coding genes > intergenic regions. Peptides that cannot be mapped to human genome are currently not considered.

Module 5: Spectrum Quality Scoring

The quality of tandem MS spectra is scored by the method of linear discriminant analysis (LDA)³⁸ with modification, using the 1^st stage database search result as a training set. The method integrates features reflecting tag information (i.e. the top tag length) and spectrum data (i.e. square root of MS2 peak number and logarithm transformed MS2 maximum peak intensity) by the equation:

Where f_i and c_i denote value and coefficient of each feature, respectively. The coefficients are obtained by training the LDA model. The threshold of spectrum quality score is adjusted to allow discarding less than 1% of accepted PSMs.

To evaluate whether the LDA method introduced overfitting during the analysis, we performed 10-fold cross validation analysis in the following steps: i) the whole dataset was randomly partitioned into 10 equal sized subsamples; ii) of these 10 subsamples, 9 subsamples were used for training the LDA model and deriving a training receiver operating characteristic (ROC) curve, while the remaining subsample was used as the testing data for a testing ROC curve; iii) we repeated the 2^nd step for 10 times to evaluate the variation of the analysis. If overfitting occurred, the area under ROC curve (AUC) of the training dataset would be larger than that of the independent testing dataset. Our results showed that the training and testing AUCs were nearly identical (Table S4 in Supporting Information), indicating no overfitting of LDA modeling in our analysis.

Data Analysis of Two Large Dataset (AD and MM cells)

The first dataset was collected from human Alzheimer’s disease postmortem brain, including label-free high resolution MS raw data (1.7 million MS/MS scans) and RNA-seq data (50 million reads; Table S1 in Supporting Information).⁵² The second dataset was obtained from a multiple myeloma cell line (ANBL6),⁵³ including tandem mass tag (TMT, Thermo Fisher Scientific) labeled raw data (6.2 million MS/MS scans) and RNA-seq data (242 million reads).

Proteogenomics analysis was performed by JUMPg in three stages. i) MS data were searched against a database combining mutations and splice junction peptides as well as UniProt human proteins (downloaded in February 2015) with fully tryptic restriction. The resulting PSMs in the 1^st stage were filtered to 1% unique protein FDR. ii) Partially tryptic search was carried out with the unmatched high quality spectra against proteins accepted during the 1^st stage. iii) The remaining spectra were matched to RNA-seq 6FT database, and filtered to 0% FDR. However, we recognized that the FDR estimation may not be accurate and further calculated the standard deviation of the FDR by previously reported method.⁴⁹ The actual FDR falls between 0% and 0.2% within a 99% confidence interval. Moreover, we compared peptides identified from novel search space with those from reference database, and found that all peptides showed highly similar distributions of matching score, precursor mass error and tag length (Figure S2 in Supporting Information), indicating that our stringent filtering yielded high quality PSMs in novel search space.

Improved Peptide Identification by Multistage Customized Databases and Partially Tryptic Search

The JUMPg program implements three independent methods for constructing customized peptide databases including mutations, splice junctions and RNA-seq 6FT translated peptides (Figure 1). While mutations and splice junctions focus on variant peptides in annotated CDS regions, there is emerging evidence for functional proteins encoded by unknown open reading frames.³ Therefore JUMPg constructs another database by RNA-seq de novo assembly and six frame translation²⁵ (Figure 2). This method shows the best balance between database comprehensiveness (the whole transcriptome) and specificity (compact search space), compared to whole genome translation (Supporting Information Note, Table S2 and Figure S3). As to search space, the mutation, splice junction and 6FT databases are equivalent to 1.1%, 3.0% and 71.4% of fully tryptic UniProt protein database (default reference database), respectively (Table 1). In addition, whole proteome data often contain a significant amount of partially tryptic peptides possibly due to endogenous protease activities.⁴⁸ The partially tryptic database is ~10 times larger than the fully tryptic counterpart.

To increase the search efficiency against these databases (Figure S4 in Supporting Information), a multistage strategy is developed to perform the analysis sequentially (Figure 1A). First, the MS/MS spectra are matched to a fully tryptic peptide database pooled from UniProt, mutations and splice junctions. The matching results are filtered to ~1% false discovery rate at protein level. Then only the unmatched high quality spectra are subjected to the 2^nd stage of analysis against the large partially tryptic database. Moreover, this partially tryptic database is shrunk by considering proteins only identified in the 1^st stage (Table S3 in Supporting Information). Finally, the remaining high quality spectra are searched against the 6FT database, as the majority of the peptides inside are false owing to six frame translation.

JUMPg exports the identified peptides in a text table, showing the peptide sequences, PSM counts, tags and scores of the best PSM, and database entries (Figure 1B). JUMPg also provides a visualization function by converting peptides to genomic locations via UCSC known gene annotations, which can be co-displayed with other genomic information on the UCSC genome browser (Figure 1C).⁵⁰ For example, two novel peptides (in magenta) were identified to indicate an event of intron retention, whereas another peptide (in green) was also identified to support the canonical protein isoform, annotated with known transcript (in navy) of the CASKIN1 gene in chromosome 16. Taken together, JUMPg maximizes novel peptide identification by three different customized database methods in an effective manner, and displays results for integrative visualization.

Acceleration of Multistage Database Search by Spectrum Quality Control and Amino Acid Tag Recycling

To reduce multistage database search time, we discard low quality spectra before the 2^nd stage analysis, based on MS/MS spectrum quality. To calculate a quality score for each spectrum, we use three features, spectrum intensity, total peak number, and the best tag length. All three features are highly correlated with PSM identification successful rate defined by accepted PSMs in the 1^st stage analysis (Figure 3A), consistent with previous reports.^{39, 54} Of these, the tag length (in black) shows the best performance by receiver operating characteristic (ROC) curves (Figure 3B). To maximize the discriminant capacity, we use the LDA method³⁸ to combine these features into one single quality score (in red), which outperforms any individual feature (Figure 3B). Next, a threshold of the quality score is selected to balance the removal of MS/MS spectra and the recovery of high quality spectra, depending on the score distribution and the accumulative curve of accepted PSMs in the 1^st stage analysis (Figure 3C). When preserving 99% of accepted PSMs, the threshold is able to eliminate ~35% of the MS/MS spectra. As database search time is linearly correlated with the spectrum number,⁵⁵ this filtering step can save ~35% of computational time (Figure S4C in Supporting Information).

To further improve database search speed, JUMPg implements a novel function by recycling amino acid tags of each MS/MS spectrum. The tag information enables high sensitivity and accuracy of peptide identification,^{40, 56, 57} but the process to derive tags is time-consuming. The tags generated at the 1^st stage analysis are recorded and re-used in the 2^nd and later stage analysis, which further accelerate database search by ~25%. Taken together, the spectrum quality control and tag recycling reduce search time by more than 50% for multistage analysis.

We also considered one potential caveat of multistage analysis, in which PSMs accepted in the 1^st stage could have higher scores in the 2^nd or 3^rd stage. To address this, more than 156,000 MS/MS spectra were searched against fully tryptic database, resulting in 52,114 accepted PSMs (1% FDR). These spectra were re-searched against partially tryptic database, leading to higher scores of only 97 (0.2% out of 52,114) PSMs. This small percentage is lower than the specified FDR level (1%). Thus, the multistage strategy saves database search time by omitting PSMs accepted in the 1^st stage, which does not introduce extra false positives of PSMs.

Application to Label Free Data: Discovery of a Novel Protein Coding Gene in Human Brain

We previously analyzed Alzheimer’s disease brain proteome by long gradient LC/LC-MS/MS and the transcriptome by RNA sequencing, identifying ~10,000 reference proteins.⁵² However, the sample-specific events such as DNA polymorphisms and brain specific alternative splicing events have not been addressed. Here we used the JUMPg pipeline to analyze the same dataset (dataset 1 with ~1.7 million label free MS/MS spectra), discovering a total of 496 novel peptides in three stage analyses (Table 1). These 496 novel peptides are classified into 14 types of events (Table 2), including amino acid substitution, deletion, frame shift, exon skipping and other alternative splicing events, as well as translation from antisense strand, non-coding genes, untranslated regions (UTR) and intergenic regions. Four examples are illustrated in detail (Figure 4): (i) a single amino acid substitution in ACO2 supported by three peptides, (ii) an in-frame deletion of PRUNE2, (iii) an exon skipping in SLC12A5, and (iv) 5′ UTR translation / alternative start codon demonstrated by the original and Met oxidized forms of a fully tryptic peptide.

As a high percentage of mutated peptides identified in proteogenomics studies are also found in public single nucleotide polymorphism database (dbSNP),^{24, 31, 37} we analyzed in this sample the 255 novel single-amino acid substitution events, and found that 239 (94%) are collected in the dbSNP (version 138). This consistent result provides additional support to the confidence of JUMPg-identified novel peptides.

Interestingly, we uncovered a novel protein coding gene specifically expressed in the human brain. During the 3^rd stage analysis with the 6FT database, we found nine novel peptides that locate within a 1.5 kb region of chromosome X (Figure 5A, 5B), a genomic region with no protein coding gene models recorded in databases (including refSeq, UCSC, GENCODE v19, Ensembl v75 and AceView). Closer examination of these nine peptides (in magenta) confirmed that they share the same open reading frame, suggesting that these peptides are translated from one novel protein coding gene. To obtain the full protein sequence of this new gene, we aligned the assembled RNA transcripts to the reference genome, translated this genomic region according to the reading frame defined by JUMPg-identified peptides, and derived 578 amino acid protein sequence (in green, Figure S4 in Supporting Information).

To functionally annotate this novel protein, the sequence was aligned to NCBI non-redundant protein database.⁵⁸ The top hit was a computationally predicted protein termed “putative paraneoplastic antigen-like protein 6B-like protein (PNMA6BL)” in a primate (Cercocebus atys) with 86% identity. Moreover, one “PNMA” domain was identified by Pfam⁵⁹ within the novel protein sequence, suggesting that this protein belongs to the paraneoplastic antigens Ma (PNMA) gene family. In human, this family includes 7 known members (PNMA1, 2, 3, 5, 6A and PNMAL1, 2). Phylogenetic analysis indicated that this novel protein was clustered within the PNMA6 sub-branch (Figure 5C), consistent with its homologue in the primate. Thus, we named this novel protein as “PNMA6BL”. In addition, pan-tissue RNA-seq analysis indicates that the PNMA6BL gene is almost exclusively expressed in human brain (Figure 5D), which is typical for the PNMA gene family.⁶⁰ This example demonstrates that our JUMPg pipeline allows the identification of previously unannotated genes by combining deep proteomics and transcriptomics data.

Application to Isotopically labeled Data: Identification of Immunoglobulin Light Chain Variable Region Peptide in Multiple Myeloma Cells

Isobaric labeling techniques (e.g. iTRAQ/TMT) have been widely used in quantitative proteomics,^{61, 62} but peptide identification is more challenging due to noise peaks induced by labeling reagents.⁶³ To further test the performance of JUMPg, we applied the same procedure (Figure 1A) to a large TMT dataset (dataset 2 with 6.2 million MS/MS spectra) from multiple myeloma (MM) cell lines. Intriguingly, we observed a high percentage of partially tryptic peptides in this dataset (Table 1), with 67,368 partially tryptic peptides (33% of total peptides) from 182,988 PSMs (20% of total PSMs). As MM cells are known to produce a high level of immunoglobulin, cell survival relies on excessive proteasome activity,⁶⁴ which might explain the presence of additional partially tryptic peptides. Similarly, high frequency of partial trypticity was also detected in human serum,⁶⁵ suggesting that in vivo proteases may be linked to this elevated partial cleavage.

In the three stage analyses, we accepted 198,103 peptides from 909,223 PSMs (Table 1). As expected, the PSM identification successful rate of the TMT dataset (18%) is lower than that of the label free AD dataset (36%). Nonetheless, we determined a total of 991 novel peptides, which were summarized into 531 mutation, 158 novel splicing and 48 non-canonical translation events (Table 2). Notably, the most abundant peptide matched by 211 PSMs displays three substitutions in the immunoglobulin light chain variable region (Figure 6). The result exemplifies the advantage of the 6FT method, in which bona fide RNA information is retained in the 6FT database. Because the immunoglobulin light chain variable region is unique to individual B cell clones, the capability of identifying these unique peptides highlights JUMPg as a potential tool towards personalized proteome.⁶⁶