Introduction

Proteogenomics involves the use of mass spectrometry (MS) based proteomics data against genomics and transcriptomics data to identify peptides and to understand protein-level evidence of gene expression. In the first section of the tutorial, we will create a protein database (FASTA) using RNA-sequencing files (FASTQ) and then perform sequence database searching using the resulting FASTA file with the MS data to identify peptides corresponding to novel proteoforms. Then, we will assign the genomic coordinates and annotations for these identified peptides and visualize the data for its spectral quality and genomic localization

Proteogenomics integrates RNA-Seq data for generating customized protein sequence databases with mass spectrometry-based proteomics data, which are matched to these databases to identify protein sequence variants. (Chambers et al. 2017).

Agenda

In this tutorial, we will cover:

1. Introduction
2. Data upload
3. Change the chromosome names in the Ensembl GTF to match a UCSC genome reference
4. Aligning FASTQ files to the mouse genome
5. Variant Analysis
   1. Variant Calling
   2. Variant annotation and Genome Mapping
6. Transcript Assembly
   1. RNA-seq to transcripts
   2. Evaluate the assembly with annotated transcripts
   3. Translate transcripts
7. Creating FASTA Databases
   1. Genomic mapping database
   2. Variant Annotations database
   3. List of Reference Protein Accession names
8. Conclusion

Overview

This tutorial focuses on creating a FASTA database generated from RNA-seq data. There are two outputs from this workflow: (1) a sequence database consisting of variants and known reference sequences and (2) mapping files containing genomic and variant mapping data.

The first part of the tutorial deals with creating the FASTA file containing sequences with single amino acid variants (SAVs), insertions and deletions (indels). The second part of the workflow helps in creating a FASTA file with transcript assemblies (splicing variants).

Data upload

In this tutorial, protein and the total RNA sample was obtained from the early development of B-cells from mice. It was obtained at two developmental stages of B-cells: Ebf1 -/- pre-pro-B and Rag2 -/- pro-B. Please refer to the original study Heydarian, M. et al.for details.

Hands-on: data upload and organization

1. Create a new history and name it something meaningful (e.g. Proteogenomics DB creation)

   Tip: Creating a new history

   Click the new-history icon at the top of the history panel.

   If the new-history is missing:

   1. Click on the galaxy-gear icon (History options) on the top of the history panel
   2. Select the option Create New from the menu

2. Import the Uniprot FASTA, FASTQ file and the GTF file from Zenodo

   https://zenodo.org/api/files/bf8c34bc-ed55-4b1c-9b69-edfe0926ea84/Trimmed_ref_5000_uniprot_cRAP.fasta
   https://zenodo.org/record/1489208/files/FASTQ_ProB_22LIST.fastqsanger
   https://zenodo.org/record/1489208/files/Mus_musculus.GRCm38.86.gtf
   

   Tip: Importing via links

   • Copy the link location

   • Open the Galaxy Upload Manager (galaxy-upload on the top-right of the tool panel)

   • Select Paste/Fetch Data

   • Paste the link into the text field

   • Press Start

   • Close the window

3. Rename the datasets with more descriptive names (strip off the url prefixes)

   Tip: Renaming a dataset

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, change the Name field
   • Click the Save button

4. Make sure that the datatype of file FASTQ_ProB_22List.fastqsanger is set to fastqsanger

   Tip: Changing the datatype

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, click on the galaxy-chart-select-data Datatypes tab on the top
   • Select fastqsanger
     • tip: you can start typing the datatype into the field to filter the dropdown menu
   • Click the Save button

5. Make sure the Database/Build (dbkey) is set to Mouse.Dec 2011 (GRCm38/mm10)(mm10)

   Tip: Changing database/build (dbkey)

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, change the Database/Build field
   • Select your desired database key from the dropdown list: Mouse.Dec 2011 (GRCm38/mm10)(mm10)
   • Click the Save button

6. Note: If you are running the workflow taken from the GTN, then make sure the tool “Tabular-to-FASTA” (tool number 26) has Title Column labeled as “1” and Sequence Column as “2”.

Change the chromosome names in the Ensembl GTF to match a UCSC genome reference

UCSC prefaces chromosome names with chr while Ensembl does not. To perform this action, we use the Replace Text in a specific column - a regular expression tool.

Hands-on: Replace Text in a specific column

1. Replace Text in a specific column tool with the following parameters:
2. param-file “File to process”: Mus_musculus.GRCm38.86.gtf
3. Change numbered chromosome names
   • param-select “in column”: 1
   • param-select “Find pattern”: ^([1-9][0-9]*)$
   • param-select “Replace with: chr\\1
4. param-repeat Insert Replacement to add another replacement option to change XY chromosomes
   • param-select “in column”: 1
   • param-select “Find pattern”: ^([XY])$
   • param-select “Replace with: chr\\1
5. param-repeat Insert Replacement to add another replacement option to change mitochondrial chromosome
   • param-select “in column”: 1
   • param-select “Find pattern”: ^MT$
   • param-select “Replace with: chrM

1. Rename the output to Mus_musculus.GRCm38.86.fixed.gtf

Comment: Note on chromosome name changes

This step is necessary to help identify and correct errors due to inputs having non-identical chromosome identifiers and/or different chromosome sequence content. Galaxy Support

Aligning FASTQ files to the mouse genome

The first tool in the workflow is the HISAT2 alignment tool. It maps next-generation sequence (NGS) reads to the reference genome. This tool requires an RNA-seq file (.FASTQ) and a reference genome file in Gene transfer format (GTF). This .gtf file is obtained from the Ensembl database. When successful, this tool outputs a .bam file (binary version of a SAM: Sequence Alignment/Map).

Hands-on: Alignment with HISAT2

1. HISAT2 tool with the following parameters:
   • param-select “Source for the reference genome”: Use a built-in genome
   • param-select “Select a reference genome”: mm10
   • param-select “Single-end or paired-end reads”: Single end
   • param-file “Input FASTQ files”: FASTQ_ProB_22LIST.fastqsanger
   • param-select “Specify strand information”: Unstranded
2. Inspect galaxy-eye the resulting files
3. Rename the output to HISAT_Output.BAM

Comment: Note on strandedness

Note that if your reads are from a stranded library, you need to choose the appropriate setting for “Specify strand information” mentioned above. For single-end reads, use F or R. F means that a read corresponds to a forward transcript. R means that a read corresponds to the reverse complemented counterpart of a transcript.

For paired-end reads, use either FR or RF. With this option being used, every read alignment will have an XS attribute tag: + means a read belongs to a transcript on the positive + strand of the genome. - means a read belongs to a transcript on the negative - strand of the genome. (TopHat has a similar option, --library--type option, where fr - first strand corresponds to R and RF; fr - second strand corresponds to F and FR.)

Variant Analysis

Variant Calling

FreeBayes is a Bayesian genetic variant detector designed to find small polymorphisms, specifically SNPs (single-nucleotide polymorphisms), indels (insertions and deletions), MNPs (multi-nucleotide polymorphisms), and complex events (composite insertion and substitution events) smaller than the length of a short-read sequencing alignment.

Provided with some BAM dataset(s) and a reference sequence, FreeBayes will generate a VCF dataset describing SNPs, indels, and complex variants in samples in the input alignments. By default, FreeBayes will consider variants supported by at least two observations in a single sample (-C) and also by at least 20% of the reads from a single sample (-F). These settings are suitable for low to high depth sequencing in haploid and diploid samples, but users working with polyploid or pooled samples may adjust them depending on the characteristics of their sequencing data.

FreeBayes is capable of calling variant haplotypes shorter than a read length where multiple polymorphisms segregate on the same read. The maximum distance between polymorphisms phased in this way is determined by the --max-complex-gap, which defaults to 3bp. In practice, this can comfortably be set to half the read length. Ploidy may be set to any level (-p), but by default all samples are assumed to be diploid. FreeBayes can model per-sample and per-region variation in copy-number (-A) using a copy-number variation map.

FreeBayes can act as a frequency-based pooled caller and can describe variants and haplotypes in terms of observation frequency rather than called genotypes. To do so, use --pooled-continuous and set input filters to a suitable level. Allele observation counts will be described by AO and RO fields in the VCF output.

Hands-on: Variant Calling with FreeBayes

1. FreeBayes tool with the following parameters:
   • param-select “Choose the source for the reference genome”: Locally cached file
     • param-check “Run in batch mode?”: Run Individually
     • param-file “BAM dataset:: HISAT_Output.BAM
     • param-select “Using reference genome”: mm10
   • param-select “Limit variant calling to a set of regions?”: Do not Limit
   • param-select “Choose parameter selection level”: Simple diploid calling
2. Inspect galaxy-eye the resulting files.

Note: FreeBayes may take some time…be patient!!!!

Comment: FreeBayes options

Galaxy allows five levels of control over FreeBayes options, provided by the Choose parameter selection level menu option. These are:

1. Simple diploid calling: The simplest possible FreeBayes application. Equivalent to using FreeBayes with only a BAM input and no other parameter options.

2. Simple diploid calling with filtering and coverage: Same as #1 plus two additional options: (1) -0 (standard filters: --min-mapping-quality 30 --min-base-quality 20 --min-supporting-allele -qsum 0 --genotype-variant-threshold 0) and (2) --min-coverage.

3. Frequency-based pooled calling: This is equivalent to using FreeBayes with the following options: --haplotype-length 0 --min-alternate-count 1 --min-alternate-fraction 0 --pooled -continuous --report- monomorphic. This is the best choice for calling variants in mixtures such as viral, bacterial, or organellar genomes.

4. Frequency-based pooled calling with filtering and coverage: Same as #3 but adds -0 and --min-coverage like in #2.

Complete list of all options: Gives you full control by exposing all FreeBayes options as Galaxy parameters.

Variant annotation and Genome Mapping

CustomProDB (Wang and Zhang 2013) generates custom protein FASTA sequences files from exosome or transcriptome data. Once Freebayes creates the .vcf file, CustomProDB uses this file to generate a custom protein FASTA file from the transcriptome data. For this tool, we use Ensembl 89 mmusculus (GRm38.p5) (dbsnp142) as the genome annotation. We create three FASTA files from CustomProDB: (1) a variant FASTA file for short indels, (2) a Single Amino acid Variant (SAV) FASTA file, an SQLite database file (genome mapping and variant mapping) for mapping proteins to a genome and (3) an RData file for variant protein coding sequences. The reference protein set can be filtered by transcript expression level (RPKM calculated from a BAM file), and variant protein forms can be predicted based on variant calls (SNPs and indels reported in a VCF file).

Annotations CustomProDB depends on a set of annotation files (in RData format) to create reference and variant protein sequences. Galaxy administrators can use the CustomProDB data manager to create these annotations to make them available for users.

Hands-on: Generate protein FASTAs from exosome or transcriptome data

1. CustomProDB tool with the following parameters:
   • param-select “Will you select a genome annotation from your history or use a built-in annotation?”: Use built in genome annotation
     • param-select “Using reference genome”: Ensemble 89 mmusculus (GRm38.p5) (dbsnp142)
     • param-file “BAM file”: HISAT_Output.BAM
     • param-file “VCF file”: Freebayes.vcf
     • param-check “Annotate SNPs with rsid from dbSNP”: No
     • param-check “Annotate somatic SNPs from COSMIC (human only)”: No
   • param-text “Transcript Expression Cutoff (RPKM)”: 1
   • param-check “Create a variant FASTA for short insertions and deletions”: Yes
   • param-check “Create SQLite files for mapping proteins to genome and summarizing variant proteins”: Yes
   • param-check “Create RData file of variant protein coding sequences”: Yes
2. Inspect galaxy-eye the resulting files

There are 6 files are generated through the CustomProDB tool:

• rpkm.fasta: reference proteins filtered by transcript expression level (RPKM calculated from a BAM file)
• snv.fasta: Single Nucleotide Variants detected by FreeBayes (.vcf file)
• indel.fasta: Insertions/Deletions detected by FreeBayes (.vcf file)
• Genomic_mapping.sqlite: mapping proteins to genome
• Variant_annotation.sqlite: used to visualize the variants in the Multi-omics Visualization Platform
• RData: set of annotation files (in RData format) to create reference and variant protein sequences

For mapping proteins to genome and an RData file for variant protein coding sequences. Similar to the genomic mapping, a variant mapping file is also generated from CustomProDB. This SQLite file is also converted to tabular format and made SearchGUI-compatible. This variant annotation file will be used to visualize the variants in the Multi-omics Visualization Platform (in-house visualization platform developed by Galaxy-P senior developers).

Transcript Assembly

RNA-seq to transcripts

StringTie is a fast and highly efficient assembler of RNA-Seq alignments into potential transcripts. It uses a network flow algorithm as well as an optional de novo assembly step to assemble and quantitate full-length transcripts representing multiple splice variants for each gene locus.

Its input can include not only the alignments of raw reads used by other transcript assemblers, but also alignments of longer sequences that have been assembled from those reads. To identify differentially expressed genes between experiments, StringTie’s output can be processed by specialized software like Ballgown, Cuffdiff or other programs (DESeq2, edgeR, etc.).

Hands-on: Transcript assembly with StringTie

1. StringTie tool with the following parameters:
   • param-file “Input mapped reads”: HISAT_Output.BAM
   • param-select “Specify strand information”: Unstranded
   • param-select “Use a reference file to guide assembly?”: Use Reference GTF/GFF3
     • param-select “Reference file”: Use file from History
     • param-file “GTF/GFF3 dataset to guide assembly”: Mus_musculus.GRCm38.86.fixed.gtf
     • param-select “Use Reference transcripts only?”: No
     • param-select “Output files for differential expression?”: No additional output
     • param-select “Output coverage file?”: No
2. Inspect galaxy-eye the resulting files.

3. Rename the output to Stringtie_output.gtf

   Tip: Renaming a dataset

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, change the Name field
   • Click the Save button

4. Make sure the datatype is gtf

   Tip: Changing the datatype

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, click on the galaxy-chart-select-data Datatypes tab on the top
   • Select your desired datatype
     • tip: you can start typing the datatype into the field to filter the dropdown menu
   • Click the Save button

StringTie accepts a BAM (or SAM) file of paired-end RNA-seq reads, which must be sorted by genomic location (coordinate position). This file contains spliced read alignments and can be produced directly by programs such as HISAT2. We recommend using HISAT2 as it is a fast and accurate alignment program. Every spliced read alignment (i.e. an alignment across at least one junction) in the input BAM file must contain the tag XS to indicate the genomic strand that produced the RNA from which the read was sequenced. Alignments produced by HISAT2 (when run with the --dta option) already include this tag, but if you use a different read mapper you should check that this XS tag is included for spliced alignments.

Details: Note about parameter settings

Be sure to run HISAT2 with the --dta option for alignment (under ‘Spliced alignment options’), for the best results.

Also note that if your reads are from a stranded library, you need to choose the appropriate setting for ‘Specify strand information’ above. As, if Forward (FR) is selected, StringTie will assume the reads are from a --fr library, while if Reverse (RF) is selected, StringTie will assume the reads are from a --rf library, otherwise it is assumed that the reads are from an unstranded library (the widely-used, although now deprecated, TopHat had a similar --library-type option, where fr-firststrand corresponded to RF and fr-secondstrand corresponded to FR). If you don’t know whether your reads are from a stranded library or not, you could use the tool ‘RSeQC Infer Experiment’ to try to determine.

As an option, a reference annotation file in GTF/GFF3 format can be provided to StringTie. In this case, StringTie will prefer to use these “known” genes from the annotation file, and for the ones that are expressed it will compute coverage, TPM and FPKM values. It will also produce additional transcripts to account for RNA-seq data that aren’t covered by (or explained by) the annotation. Note that if option -e is not used, the reference transcripts need to be fully covered by reads in order to be included in StringTie’s output. In that case, other transcripts assembled from the data by StringTie and not present in the reference file will be printed as well.

We highly recommend that you provide annotation if you are analyzing a genome that is well annotated, such as human, mouse, or other model organisms.

Evaluate the assembly with annotated transcripts

GffCompare compares and evaluates the accuracy of RNA-Seq transcript assemblers (Cufflinks, Stringtie). * collapse (merge) duplicate transcripts from multiple GTF/GFF3 files (e.g. resulted from assembly of different samples) * classify transcripts from one or multiple GTF/GFF3 files as they relate to reference transcripts provided in a annotation file (also in GTF/GFF3 format)

The original form of this program is also distributed as part of the Cufflinks suite, under the name CuffCompare. Most of the options and parameters of CuffCompare are supported by GffCompare, while new features will likely be added to GffCompare in the future.

Hands-on: compare assembled transcripts against a reference annotation

1. GffCompare tool with the following parameters:
   • param-file “GTF inputs for comparison”Stringtie_output.gtf
   • param-select “Use Reference Annotation”: yes
     • param-select “Choose the source for the reference annotation”: History
       • param-file “Reference Annotation”: Mus_musculus.GRCm38.86.gtf
     • param-select “Ignore reference transcripts that are not overlapped by any input transfrags”: No
     • param-select “Ignore input transcripts that are not overlapped by any reference transcripts”: No
   • param-select “Use Sequence Data”: No
   • param-select “Discard (ignore) single-exon transcripts”: No
   • param-text “Max. Distance for assessing exon accuracy”: 100
   • param-text “Max distance for transcript grouping”: 100
   • param-select “discard intron-redundant transfrags sharing 5’“: No

Details: GffCompare vs CuffCompare

A notable difference between GffCompare and CuffCompare is that when a single query GTF/GFF file is given as input along with a reference annotation (-r option), GFFCompare switches into “annotation mode” and it generates a .annotated.gtf file instead of the .merged.gtf produced by CuffCompare with the same parameters. This file has the same general format as CuffCompare’s .merged.gtf file (with “class codes” assigned to transcripts as per their relationship with the matching/overlapping reference, transcript) but the original transcript IDs are preserved, so GffCompare can thus be used as a simple way of annotating a set of transcripts.

Another important difference is that the input transcripts are no longer discarded when they are found to be “intron redundant”, i.e. contained within other, longer isoforms. CuffCompare had the -G option to prevent collapsing of such intron redundant isoforms into their longer “containers”, but GffCompare has made this the default mode of operation (hence the -G option is no longer needed and is simply ignored when given).

Translate transcripts

Hands-on: Translate transcripts to protein sequence

First, we must convert the GffCompare annotated GTF file to BED format.

1. Convert gffCompare annotated GTF to BED tool with the following parameters:
   • param-file “GTF annotated by gffCompare”: output from gff compare
   • param-select “filter GffCompare class_codes to convert”:
   • j : Potentially novel isoform (fragment): at least one splice junction is shared with a reference transcript
   • e : Single exon transfrag overlapping a reference exon and at least 10 bp of a reference intron, indicating a possible pre-mRNA fragment.
   • i : A transfrag falling entirely within a reference intron
   • p : Possible polymerase run-on fragment (within 2Kbases of a reference transcript)
   • u : Unknown, intergenic transcript

   Next, we translate transcripts from the input BED file into protein sequences.

2. Translate BED transcripts cDNA in 3frames or CDS tool with the following parameters:
   • param-file “A BED file with 12 columns”: Convert GffCompare-annotated GTF to BED
   • param-select “Source for Genomic Sequence Data”: Locally cached File
   • param-select “Select reference 2bit file”: mm10
   • param-select “Fasta ID Options”
   • param-select “fasta ID source, e.g. generic”: generic

   Finally, we convert a BED format file of the proteins from a proteomics search database into a tabular format for the Multiomics Visualization Platform (MVP).

3. bed to protein map tool with the following parameters:
   • param-file “A BED file with 12 columns, thickStart and thickEnd define protein coding region”: Translate cDNA_minus_CDS

Comment: Note on data formats

• The tabular output can be converted to a sqlite database using the Query_Tabular tool.
• The sqlite table should be named: feature_cds_map.
• The names for the columns should be: name, chrom, start, end, strand, cds_start, cds_end
• This SQL query will return the genomic location for a peptide sequence in a protein (multiply the amino acid position by 3 for the cds location)

Creating FASTA Databases

In this section we will perform the following tasks:

• The Protein Database Downloader tool is used to download the FASTA database from UniProt and cRAP (common Repository of Adventitious Proteins) database containing known/reference mouse proteins. For the purpose of this tutorial, we will use the Reference_5000_uniprot_cRAP.fasta file as this database.
• The FASTA Merge Files and Filter Unique Sequences tool is used to merge the databases obtained from CustomProDB and “Translate BED” tool along with the UniProt and cRAP databases.
• The Regex Text Manipulation tool is used to manipulate the FASTA file to make it SearchGUI-compatible.

Hands-on: Hands-on

1. FASTA Merge Files and Filter Unique Sequences tool
   • Click param-repeat “Insert Input FASTA File(s)” twice so that there are a total of three “Input FASTA File(s)” blocks.
   • Use the following parameters:
   • param-check “Run in batch mode?”: Merge individual FASTAs (output collection if input is collection)
   • param-file “Input FASTA File(s)” : Input Custom ProDB Fasta File output

     1.HISAT_Output.rpkm
     2.HISAT_Output.snv
     3.HISAT_Output.indel
     

   • param-select “How are sequences judged to be unique?”: Accession and Sequence
   • param-text “Accession Parsing Regular Expression”: ^>([^ |]+).*$

Comment: Tool parameters explained

This tool concatenates FASTA database files together.

• If the uniqueness criterion is “Accession and Sequence”, only the first appearence of each unique sequence will appear in the output. Otherwise, duplicate sequences are allowed, but only the first appearance of each accession will appear in the output.
• The default accession parser will treat everything in the header before the first space as the accession.

Comment: CustomProDB generated protein Accession names need to be changed for PeptideShaker

PeptideShaker does not allow greaterthan sign, comma, or spaces in the protein accession name. PeptideShaker also requires the FASTA ID to start with a source, e.g. >sp|Q8C4J7| for SwissProt entries. Nonstandard entries must be prefixed generic|. We will also need to similarily modify the accessions in the genomic mapping and and variant annotation outputs. For eg: ENSEMBL-PRO (output from CustomProDB) and STRG (from StringTie) are not standard format.

Hands-on: FASTA to Tabular

1. FASTA-to-Tabular tool: with the default parameters
   • param-file “Convert these sequences”: Merged and Filtered FASTA from (fasta)
2. Column Regex Find And Replace tool with the following parameters:
   • param-file “Select cells from”: FASTA-to-Tabular output from previous step
   • param-select “Using”: column 1
   • param-repeat Insert Check
     • param-text “Find Regex” : ^(ENS[^_]+_\d+:)([ACGTacgt]+)>([ACGTacgt]+)\s*
     • param-text “Replacement”: \1\2_\3
   • param-repeat Insert Check
     • param-text “Find Regex”: ,([A-Y]\d+[A-Y]?)\s*
     • param-text “Replacement”: .\1
   • param-repeat Insert Check
     • param-text “Find Regex”: ^(ENS[^ |]*)\s*
     • param-text “Replacement”: generic|\1
3. Tabular-to-FASTA tool with the following parameters:
   • “Title column”: 1
   • “Sequence Column”:2 - Rename it as CustomProDB Merged Fasta

Comment: FASTA-to-Tabular and Tabular-to-FASTA versions

Please make sure the tool version is 1.1.1 Also, while running the workflow, the Title column for Tabular-to-FASTA might come empty, so please make sure you select Title Column as “1”.

For visualization purposes we also use the concatenate tool to concatenate the genomic mapping with the protein mapping dataset. This output will be used for visualization in MVP to view the genomic coordinates of the variant peptide.

We also need to modify the CustomProDB protein accessions the same as was done for the FASTA database.

An SQLite database containing the genomic mapping SQLite, variant annotation and information from the protein mapping file is concatenated to form a single genomic mapping SQLite database later used as an input for the ‘Peptide Genomic Coordinate’ tool. For that we need to follow the steps below:

Genomic mapping database

Hands-on: Create Database for genomic mapping

1. SQLite to tabular for SQL query hands_on with the following parameters:
   • param-file “SQLite Database”: genomic_mapping.sqlite from CustomProDB
   • param-text “Query”:

       SELECT pro_name, chr_name, cds_chr_start - 1, cds_chr_end,strand,cds_start - 1, cds_end
       FROM genomic_mapping
       ORDER BY pro_name, cds_start, cds_end
     

   • param-file “Omit column headers from tabular output”: No

2. Rename output to genomic_mapping_sqlite

   Tip: Renaming a dataset

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, change the Name field
   • Click the Save button

   The output is further processed so that the results are compatible with the Multiomics Visualization Platform. We need to modify the CustomProDB protein accessions the same as was done for the FASTA database.

3. Column Regex Find And Replace tool with the following parameters:
   • param-file “Select cells from”: genomic_mapping_sqlite' (tabular)
   • param-select “Using”: column 1
   • param-repeat Insert Check
     • param-text “Find Regex” : ^(ENS[^_]+_\d+:)([ACGTacgt]+)>([ACGTacgt]+)\s*
     • param-text “Replacement”: \1\2_\3
   • param-repeat Insert Check
     • param-text “Find Regex”: ,([A-Y]\d+[A-Y]?)\s*
     • param-text “Replacement”: .\1
   • param-repeat Insert Check
     • param-text “Find Regex”: ^(ENS[^ |]*)\s*
     • param-text “Replacement”: \1

4. Rename the output to SAV_INDEL

   This tool goes line by line through the specified input file and if the text in the selected column matches a specified regular expression pattern, it replaces the text with the specified replacement.

   Next, we will concatenate the output from this tool with the “Bed to protein map” output.

5. Concatenate multiple datasets tool with the following parameters:
   • param-files “Concatenate Datasets”: Select the output from the previous tool named as SAV_INDEL and the output from the “bed to protein map” tool and concatenate the tool.

   Output will be the “Genomic_Protein_map” rename it as that

   Now, we load this tabular datasets into an SQLite database.

6. Query Tabular using SQLite tool with the following parameters:
   • param-repeat Insert Database Table
     • param-file “Tabular Dataset for Table”: Genomic_Protein_map
     • Section Table Options:
       • param-text “Specify Name for Table”: feature_cds_map
       • param-text “Specify Column Names (comma-separated list)”: name,chrom,start,end,strand,cds_start,cds_end
       • “Only load the columns you have named into database”: No
       • param-repeat Insert Table Index:
         • param-check “This is a unique index”: No
         • param-text “Index on columns”: name,cds_start,cds_end
       • param-select “Save the sqlite database in your history”:Yes

7. Rename the output as “genomic_mapping_sqlite”

   

Variant Annotations database

We will repeat the process for the variant annotations (this file will be used as a input to the new updates that will be made on the Multiomics Visualization platform). We need to modify the CustomProDB protein accessions the same as was done for the FASTA database.

Hands-on: Create database for variant annotations

1. SQLite to tabular tool with the following parameters:
   • param-file “SQLite Database”: variant_annotations.sqlite from CustomProDB
   • param-text “Query”:

     SELECT var_pro_name,pro_name,cigar,annotation
     FROM variant_annotation
     

   We will subject the output to text manipulation so that the results are compatible with the Multiomics Visualization Platform.

2. Column Regex Find And Replace tool with the following parameters:
   • param-file “Select cells from”: variant_annotations_sqlite' (tabular)
   • param-select “Using:” column 1
   • param-repeat Insert Check
     • param-text “Find Regex”: ^(ENS[^_]+_\d+:)([ACGTacgt]+)>([ACGTacgt]+)\s*
     • param-text “Replacement”: \1\2_\3
   • param-repeat Insert Check:
     • param-text “Find Regex”: ,([A-Y]\d+[A-Y]?)\s*
     • param-text “Replacement”: .\1
   • param-repeat Insert Check:
     • param-text “Find Regex”: ^(ENS[^ |]*)\s*
     • param-text “Replacement”: \1

   Next, we will load this tabular dataset into an SQLite database.

3. Query Tabular tool with the following parameters:
   • param-repeat Insert Database Table
     • param-file “Tabular Dataset for Table”: variant_annotation
     • Section Table Options:
       • param-text “Specify Name for Table”: variant_annotation
       • param-text “Specify Column Names (comma-separated list)”: name,reference,cigar,annotation
       • param-check “Only load the columns you have named into database”: No
     • param-repeat Insert Table Index
       • param-check “This is a unique index”: No
       • param-text “Index on columns”: name,cigar
     • param-select “Save the sqlite database in your history”:Yes
4. Rename the output as “Variant_annotation_sqlitedb”

These SQLite databases, which contain the genomic mapping SQLite and variant annotation information from the protein mapping file, will be utilized by MVP to visualize the genomic loci of any variant peptides.

Finally, we can create a database which can be used to search Mass spectrometry (raw/mgf) files to identify peptides that match to this database.

To do so:

Hands-on: Hands-on

1. FASTA Merge Files and Filter Unique Sequences tool with the following parameters:
   • param-check “Run in batch mode?”: Merge individual FASTAs (output collection if input is collection)
   • param-files “Input FASTA File(s)” : Input Custom ProDB Merged Fasta File output
   • param-files “Input FASTA File(s)” : Input Translate BED transcripts Fasta output
   • param-files “Input FASTA File(s)” : Input Reference_5000_uniprot_cRAP.fasta
   • param-select “How are sequences judged to be unique?”: Accession and Sequence
   • param-text “Accession Parsing Regular Expression”: ^>([^ |]+).*$

Comment: Tool parameters explained

This tool concatenates FASTA database files together.

• If the uniqueness criterion is “Accession and Sequence”, only the first appearence of each unique sequence will appear in the output. Otherwise, duplicate sequences are allowed, but only the first appearance of each accession will appear in the output.
• The default accession parser will treat everything in the header before the first space as the accession.

List of Reference Protein Accession names

Generate a list of Reference Proteins. Identify peptides that are contained in the reference proteins will be filtered out in the next tutorial.

Hands-on: FASTA to Tabular

1. FASTA-to-Tabular tool: with the default parameters
   • param-file “Convert these sequences”: HISAT_Output.rpkm (fasta)
2. Filter Tabular tool: with these filters:
   • param-file “Tabular Dataset to filter”: HISAT_Output.rpkm (tabular)
   • param-text “Filter by”: select columns
   • param-text “enter column numbers to keep”: 1 - param-text “Filter by”: regex replace value in column
   • param-text “enter column number to replace”: 1
   • param-text “regex pattern”: ^([^ |]+).*$
   • param-text “replacement expression”: \1
3. FASTA-to-Tabular tool: with the default parameters
   • param-file “Convert these sequences”: Trimmed_ref_5000_uniprot_cRAP.fasta (fasta)
4. Filter Tabular tool: with these filters:
   • param-file “Tabular Dataset to filter”: Trimmed_ref_5000_uniprot_cRAP.fasta (tabular)
   • param-text “Filter by”: select columns
   • param-text “enter column numbers to keep”: 1 - param-text “Filter by”: regex replace value in column
   • param-text “enter column number to replace”: 1
   • param-text “regex pattern”: ^[^|]+[|]([^| ]+).*$
   • param-text “replacement expression”: \1
5. Concatenate multiple datasets tool with the following parameters:
   • param-files “Concatenate Datasets”: Select the output from the previous 2 “Filter Tabular” outputs.

Comment: FASTA-to-Tabular and Tabular-to-FASTA versions

Please make sure the tool version is 1.1.1

Comment: Tool versions

• All the tools mentioned in this tutorial are subjected to change when the tool version is upgraded .

Conclusion

This completes the walkthrough of the proteogenomics database creation workflow. This tutorial is a guide to have a database and mapping files ready for Database searching and novel peptide analysis. Researchers can use this workflow with their data also, please note that the tool parameters, reference genomes and the workflow will be needed to be modified accordingly.

This workflow was developed by the Galaxy-P team at the University of Minnesota. For more information about Galaxy-P or our ongoing work, please visit us at galaxyp.org

Key points

• Generating variant protein database

• Generating genomic and variant mapping files for visualization

Frequently Asked Questions

Have questions about this tutorial? Check out the tutorial FAQ page or the FAQ page for the Proteomics topic to see if your question is listed there. If not, please ask your question on the GTN Gitter Channel or the Galaxy Help Forum

Useful literature

Further information, including links to documentation and original publications, regarding the tools, analysis techniques and the interpretation of results described in this tutorial can be found here.

References

1. Wang, X., and B. Zhang, 2013 customProDB: an R package to generate customized protein databases from RNA-Seq data for proteomics search. Bioinformatics 29: 3235–3237. 10.1093/bioinformatics/btt543
2. Chambers, M. C., P. D. Jagtap, J. E. Johnson, T. McGowan, P. Kumar et al., 2017 An Accessible Proteogenomics Informatics Resource for Cancer Researchers. Cancer Research 77: e43–e46. 10.1158/0008-5472.can-17-0331

Feedback

Did you use this material as an instructor? Feel free to give us feedback on how it went.
Did you use this material as a learner or student? Click the form below to leave feedback.

Citing this Tutorial

1. Subina Mehta, Timothy J. Griffin, Pratik Jagtap, Ray Sajulga, James Johnson, Praveen Kumar, 2022 Proteogenomics 1: Database Creation (Galaxy Training Materials). https://training.galaxyproject.org/training-material/topics/proteomics/tutorials/proteogenomics-dbcreation/tutorial.html Online; accessed TODAY
2. Batut et al., 2018 Community-Driven Data Analysis Training for Biology Cell Systems 10.1016/j.cels.2018.05.012

BibTeX

@misc{proteomics-proteogenomics-dbcreation,
author = "Subina Mehta and Timothy J. Griffin and Pratik Jagtap and Ray Sajulga and James Johnson and Praveen Kumar",
title = "Proteogenomics 1: Database Creation (Galaxy Training Materials)",
year = "2022",
month = "10",
day = "18"
url = "\url{https://training.galaxyproject.org/training-material/topics/proteomics/tutorials/proteogenomics-dbcreation/tutorial.html}",
note = "[Online; accessed TODAY]"
}
@article{Batut_2018,
    doi = {10.1016/j.cels.2018.05.012},
    url = {https://doi.org/10.1016%2Fj.cels.2018.05.012},
    year = 2018,
    month = {jun},
    publisher = {Elsevier {BV}},
    volume = {6},
    number = {6},
    pages = {752--758.e1},
    author = {B{\'{e}}r{\'{e}}nice Batut and Saskia Hiltemann and Andrea Bagnacani and Dannon Baker and Vivek Bhardwaj and Clemens Blank and Anthony Bretaudeau and Loraine Brillet-Gu{\'{e}}guen and Martin {\v{C}}ech and John Chilton and Dave Clements and Olivia Doppelt-Azeroual and Anika Erxleben and Mallory Ann Freeberg and Simon Gladman and Youri Hoogstrate and Hans-Rudolf Hotz and Torsten Houwaart and Pratik Jagtap and Delphine Larivi{\`{e}}re and Gildas Le Corguill{\'{e}} and Thomas Manke and Fabien Mareuil and Fidel Ram{\'{\i}}rez and Devon Ryan and Florian Christoph Sigloch and Nicola Soranzo and Joachim Wolff and Pavankumar Videm and Markus Wolfien and Aisanjiang Wubuli and Dilmurat Yusuf and James Taylor and Rolf Backofen and Anton Nekrutenko and Björn Grüning},
    title = {Community-Driven Data Analysis Training for Biology},
    journal = {Cell Systems}
}
                   

Congratulations on successfully completing this tutorial!

Do you want to extend your knowledge? Follow one of our recommended follow-up trainings:

• Proteomics
  • Proteogenomics 2: Database Search: tutorial hands-on
  • Proteogenomics 3: Novel peptide analysis: tutorial hands-on