Semidán Robaina, September 2022.

In this Notebook, we will use Pynteny through its command-line interface to find candidate peptide sequences beloging to the sox operon within unannotated, fully sequenced, genomes of marine prokaryotic organisms. Note that we could have conducted the same search through Pynteny's python API. Find more info in the documentation pages. Let's start by importing some required modules.

Download PGAP profile HMM database¶

First, let's download PGAP's profile HMM database from the NCBI webpage. To this end, we will use pynteny subcommand download, which will unzip and store files in the specified output directory. The metadata file will be parsed and filtered to remove HMM entries that are not available in the downloaded database (this is to avoid possible downstream errors).

Build peptide sequence database¶

For this example, we are going to use the MAR reference database (currently version v7), a collection of 970 fully sequenced prokaryotic genomes from the marine environment. Specifically, we will use the assembly data file containing the assembled nucleotide sequences.

Our final goal is to build a peptide sequence database in a single FASTA file where each record corresponds to an inferred ORF, which will display the positional information (i.e. ORF number within the parent contig as well as the DNA strand). To this end, we will run pynteny's subcommand build, which will take care of:

• Predict and translate ORFs with prodigal
• Label each ORF with a unique identifier and add positional metadata (with respect to the parent contig)

To follow this example, you should have previously downloaded the assembly data file, assembly.fa, from MAR ref. Here is what the first lines of assembly.fa look like, each record corresponds to a single, assembled genome:

Let's run pynteny build to generate a peptide database labeled with positional information. The labels are organized following the structure:

<genome ID>__<contig ID>_<gene position>_<locus start>_<locus end>_<strand>

where gene position, locus start, and locus end are taken with respect to the contig.

Here are some position-labeled predicted peptides corresponding to the assembled genome displayed above (CP000435.1):

Search synteny structure in MAR ref¶

Finally, we are going to use pynteny's search subcommand to search for a specific syntenic block within the previously built peptide database. Specifically, we are interested in the sox operon:

>soxX 0 >soxY 0 >soxZ 0 >soxA 0 >soxB 0 >soxC

We this synteny structure, we are searching for peptide sequences matching the profile HMM corresponding to these gene symbols, which are also arranged in this particular order, all in the positive (sense) strand, as indicated by >, and which are located exactly next to each other in the same contig (no ORFs allowed between them, as indicated by a maximum number of in-between ORFs of 0 in all cases.)

We see that pynteny parse has found a number of profile HMMs matching the gene symbols in the provided synteny structure. Additionally, in two cases it has found two HMMs matching a single gene symbol, which are displayed within parentheses and separated by "|":

>TIGR04485.1 0 >TIGR04488.1 0 >TIGR04490.1 0 >(TIGR01372.1|TIGR04484.1) 0 >(TIGR01373.1|TIGR04486.1) 0 >TIGR04555.1

In these cases, pynteny search will match sequences by either or all of the HMMs in each group within parentheses.

Alright, now that we know that our HMM database contains models for all the gene symbols in our synteny structure, let's run pynteny search to find matches in our peptide sequence database.

Some notes:

• Since we are using gene symbols instead of HMM names, we need to add the flag --gene_ids

• We could have directly input the synteny string composed of HMM names. In that case, we wouldn't need to provide the path to the HMM metadata file (--hmm_meta) and we would remove the flag --gene_ids

• Providing a path to the HMM database directory (--hmm_dir) is optional. If not provided, then pynteny will download and store the PGAP HMM database (only once if not previously downloaded) and use it to run the search. A custom HMM database provided in --hmm_dirwill override pynteny's default database

• We can also manually download the PGAP HMM database with the subcommand pynteny download

Pynteny has generated a number of output files in the provided output directory. HMMER3 hit results are stored within the subdirectory hmmer_outputs. The main output file, synteny_matched.tsv contains the labels of the matched sequences grouped by synteny block and sorted by gene number within their parent contig. The remaining (FASTA) files contain the retrieved peptide sequences for each gene symbol / HMM name in the synteny structure.

Displaying the first synteny match¶

So far, we have identified sequences (putatively) belonging to the sox operon in the MAR ref database. However, we don't know to which organisms these sequences belong. Luckily, the MAR ref database provides a metadata file that contains GTDB taxonomical information for each genome. Let's extract that info and add it to Pynteny's output table.

Displayed above the first synteny match in our peptide database, we see that all peptides are located within the same parent contig and respect the positional restrictions of our input synteny structure. Furthermore, all sequences belong to Ruegeria pomeroyii, and alphaproteobacteria of the family Rhodobacteraceae.

Additional notes:

• The previous results are strand-specific (all ORFs must be located in the positive or sense strand). However, we could have made them strand-agnostic by omitting the strand symbols in the synteny structure (i.e., using soxX 0 soxY 0 soxZ 0 soxA 0 soxB 0 soxC)

• We could have made the search even more general dropping the constraint on the arrangement by adding the flag pynteny search --unordered. In which case, Pynteny would match any group of 6 ORFs corresponding to the provided HMM names, located in the same contig and adjacent to each other, but not necessarily arranged in the same order displayed by the synteny structure. In other words, --unordered enables searching for "true" synteny, as opposed to the, more restrictive, collinearity.

Alright, finally, let's get the peptide sequences in our original input database that correspond to the identified synteny block displayed above:

There you go! All six peptide sequences corresponding to the synteny block >soxX 0 >soxY 0 >soxZ 0 >soxA 0 >soxB 0 >soxC. We have come from unannotated nucleotide assembly data to identify and locate (part of) the sox operon of Ruegeria pomeroyi. Not bad...

Here is the KEGG genome visualization centered around the sox operon. Click on the image to open it in the browser, then change gene identifiers from KID to Symbol to display gene symbols. The sox operon is depicted in purple.