Pathways Project: Annotation Walkthrough

We are now viewing the sequence for the 182 amino acids in the translated protein of Rheb-RA in D. melanogaster (Figure 14, right).

In Part 2.1, we retrieved the protein sequence for Rheb-PA in D. melanogaster. Now we need to perform a BLAST search (Figure 16) of Rheb-PA against the entire D. yakuba genome assembly.

Since we are looking for the orthologous region of Rheb in D. yakuba, we want to use BLAST to search the entire genome of D. yakuba to identify regions of local similarity with the protein sequence of Rheb in D. melanogaster we obtained in Part 2.1. In other words, we need to perform a BLAST search of our Rheb protein sequence from D. melanogaster against the entire genome of D. yakuba to narrow down the possible regions where the Rheb ortholog could be in D. yakuba.

For this BLAST search, we will use the tblastn program to search the translated nucleotide database of our target species, D. yakuba, using the protein sequence of Rheb in D. melanogaster as our query.

When performing a search, BLAST may return any number of matches (often referred to as “hits”) for regions of local similarity between our query sequence and database; however, each hit is not necessarily statistically significant. BLAST provides statistical scores to help us determine which alignments between the two sequences are statistically significant and which are spurious (i.e., likely occurred by chance and, therefore, are not evidence of real biological conservation). If BLAST returns multiple good hits (i.e., more than one match with a low E-value and a high sequence identity), we will need to investigate them all further to determine the most likely ortholog.

Our tblastn search found five regions within the translated D. yakuba genome that show similarities with the protein sequence of Rheb in D. melanogaster (Figure 21); however, only one of these is a good hit (Drosophila yakuba strain Tai18E2 chromosome 3R, Prin_Dyak_Tai18E2_2.1, whole genome shotgun sequence; sequence identity: 97.14% and E-value: 2e-78^[1]). The second hit has a much higher E-value (7e-37) and much lower percent identity (43.75%), and this pattern of increasingly lower quality matches continues through the other three matches. Therefore, we will continue our analysis based on the hypothesis that the putative ortholog of Rheb-PA in D. yakuba is in the scaffold of chromosome 3R.

Notice that each of the genome regions of our five hits has a unique accession number (Figure 20).

The accession number for the chromosome 3R scaffold in D. yakuba is NC_052530 (Figure 21, red arrow).

We should now see 7 ranges listed in ascending order of the subject start position (coordinate). Each range corresponds to a contiguous portion of the subject sequence (D. yakuba genome) that shows significant sequence similarity (i.e., matches) with a portion of the query sequence (D. melanogaster Rheb-PA). Here, the 7 ranges show matches between Rheb-PA in D. melanogaster (Query) and a range of coordinates from the NC_052530 scaffold in D. yakuba (Sbjct). For example, Range 1 shows a match between Rheb-PA in D. melanogaster (Query: 6 –​ 44) and coordinates 11,148,568 –​ 11,148,431 of NC_052530 scaffold in D. yakuba when translated in Frame -3. Range 1 has an E-value (Expect) of 0.002 and a sequence identity of 46% (Figure 24, green boxes).

Since the NC_052530 scaffold in D. yakuba is 30,730,773 base pairs (bp) long (Figure 24, purple circle), it is possible we will have some ranges (i.e., alignment matches or hits) that don’t correspond to our ortholog; therefore, we need to examine each match more closely.

Remember that we are looking for matches with low E-values and high sequence identities, and there are five matches (ranges 3 –​ 7) that fit these criteria (E-value of 2e-78 and sequence identities that range from 83% to 97%). These five alignment matches are also collinear and appear on the same strand of DNA (+ Frame) (Figure 25).

The best collinear set of alignments to Rheb-PA is located at 19,150,809 –​ 19,151,699 on the NC_052530 scaffold of the D. yakuba genome assembly and the five alignment matches cover all 182 amino acids of Rheb-PA (Figure 25, arrow). Therefore, we will continue our analysis based on the hypothesis that the putative ortholog of Rheb-PA is located at approximately 19,150,809 –​ 19,151,699 on the NC_052530 scaffold of the D. yakuba genome assembly.

In the Genome Browser image, we should now see the following tracks (Figure 29):

Looking at the NCBI RefSeq Genes track, we notice that the alignment for the coding regions of the D. yakuba RefSeq transcript XM_039375862 against the NC_052530 scaffold of D. yakuba line up with the Spaln alignment to the D. melanogaster proteins Rheb-PA and Rheb-PB. Furthermore, the coding portions of the five alignment blocks are in congruence with the placements of the five coding exons (CDS’s) predicted by GeMoMa, N-SCAN, and Augustus.

According to the RefSeq transcript XM_039375862 (Figure 29, red arrow), the putative (probable) ortholog of Rheb-PA is located at NC_052530:19,150,510-19,152,080 in D. yakuba and the region spanning NC_052530:19,150,809-19,151,702, within the putative ortholog, corresponds to the alignment to the coding region of the RefSeq transcript.

Notice the position of XM_039375862 is “NC_052530:19,150,510-19,152,080” (Figure 31, black arrow).

We are now viewing the computationally predicted protein sequence of the putative ortholog to Rheb-PA in D. yakuba.

When we run BLAST for this protein, we can enter the accession number “XP_039231796” and the program will use the translated protein sequence shown in the grey box in Figure 32.

Our blastp search found 45 proteins within the D. melanogaster reference proteins (refseq_protein) database that show similarities with the protein sequence of the putative ortholog of Rheb-PA in D. yakuba (XP_039231796); however, only one of these is a good hit (Ras homolog enriched in brain, isoform A [Drosophila melanogaster]; accession: NP_730950, sequence identity of 97.25%, and an E-value of 6e-131). The remaining hits had much higher E-values (5e-41 to 0.002) and much lower percent identities (43.75% to 23.48%) (Figure 36). Therefore, the blastp search result shows that the D. yakuba RefSeq prediction XP_039231796 is most similar to the A isoform of Rheb among all the annotated proteins in D. melanogaster. Hence, based on parsimony, the blastp search result supports the hypothesis that XP_039231796 is the putative ortholog of Rheb in D. yakuba.

Now that we’ve examined the genomic neighborhood of the putative ortholog (Figure 39), we need to identify the direction of transcription for the putative ortholog of Rheb-PA and the neighboring genes in D. yakuba and then use this information to draw a sketch of the genomic neighborhood. To do so, we need to repeat the process we followed in Part 1 (i.e., zoom into an intron (or exon) of each gene and draw arrows in the correct directions on our sketch).

If the genomic neighborhood looks similar between Rheb in D. melanogaster (Figure 11) and the putative ortholog in D. yakuba (Figure 40), we can be confident we have found the best candidate for the ortholog. However, if any of the information is inconsistent with this being a locally syntenic region (local synteny: conservation of genomic neighborhood), we should inspect our other hits in the tblastn search (Part 2) to see if a different genomic region is a better match overall.

Examination of the genomic regions surrounding the Rheb gene in D. melanogaster (Figure 41; top) and the putative Rheb ortholog in D. yakuba (Figure 41; bottom) shows that the relative gene order (i.e., CG12746, CG2931, Rheb, CRMP, and CG2926) and orientations (+, -, +, + , -) are the same in the two species. Hence the synteny analysis supports the assignment of the D. yakuba feature at NC_052530:19,150,510-19,152,080 as an ortholog of Rheb.

We need to ascertain whether the tblastn alignment for CDS-1 at 19,150,809-19,150,856 (Figure 57) is supported by RNA-Seq data and, if so, determine the location of the start codon.

The “Aggregated RNA-Seq Coverage” track show high RNA-Seq read depth within the tblastn alignment block (19,150,809-19,150,856), consistent with the hypothesis that this region is being transcribed in D. yakuba.

Notice that Frame +1 of this region has two stop codons (denoted in red with an “*”) while Frames +2 and +3 have an Open Reading Frame (i.e., no stop codons are shown within Frames +2 or +3 of CDS-1).

In Part 5 (Figure 54), we found our tblastn alignment for CDS-1 begins at approximately 19,150,809 (Figure 57, highlighted blue). Examination of this region using the Genome Browser shows us that Frame +3 has a start codon in this location and is supported by multiple evidence tracks—​the NCBI RefSeq Genes and Spaln alignment of D. melanogaster proteins, and gene predictors GeMoMa, N-SCAN, and Augustus (Figure 57).

The tblastn alignment for CDS-1 (CDS 1_9812_0) of Rheb in D. melanogaster against the D. yakuba scaffold encompasses all 16 amino acids of the CDS (Figure 54), and the alignment begins with a start codon at 19,150,809-19,150,811. Hence the most parsimonious gene model for Rheb-PA in D. yakuba would use the start codon at 19,150,809-19,150,811 in scaffold NC_052530.

In Part 5 (Figure 55), we found our tblastn alignment for CDS-5 ends at approximately 19,151,702 when translated in Frame +3. The alignment covers all 30 amino acids of CDS-5 and ends with a stop codon (Figure 58).

The stop codon at 19,151,700-19,151,702 is consistent with the NCBI RefSeq Genes and Spaln alignment of D. melanogaster proteins (stop codons don’t become part of the protein rather they signal when translation should terminate and the newly made protein should be released), and gene predictors GeMoMa, N-SCAN, and Augustus.

Based on the tblastn alignment for CDS-5 and the available evidence on the Genome Browser, the stop codon for the Rheb-PA ortholog is placed at 19,151,700-19,151,702, and the last codon (S; Serine), before the stop codon, ends at 19,151,699 (Figure 59).

Note that the “Aggregated RNA-Seq Coverage” track indicates that transcription extends beyond the stop codon. Based on the gene structure of Rheb-RA in D. melanogaster, the region with RNA-Seq read coverage that extends beyond the stop codon likely corresponds to the 3' untranslated region (UTR) of the last exon in Rheb.

Because the tblastn alignment for CDS-1 of Rheb terminates at 19,150,856 (Figure 54), we expect to find the splice donor site for CDS-1 at around position 19,150,856.

The GT splice donor site closest to 19,150,856 (Figure 61, green box) is located at 19,150,858-19,150,859 (Figure 61, red box) which is supported by multiple lines of evidence—​NCBI RefSeq Genes, Spaln alignment of D. melanogaster proteins, and the GeMoMa, N-SCAN, and Augustus gene predictions, and the aggregated RNA-Seq read coverage.

Using tblastn in Part 5, we determined the approximate location of CDS-1 was 19,150,809-19,150,856. Using the Genome Browser, in Part 6.1 we verified the start codon of CDS-1 was in Frame +3 at 19,150,809-19,150,811 (Figure 57). We visually inspected the region surrounding the approximate end of CDS-1 and placed the splice donor site at 19,150,858-19,150,859, after which we confirmed Frame +3 of CDS-1 has an ORF. Now we need to determine the phase of the splice donor site for CDS-1.

In Figure 62, we see the last complete codon (i.e., containing three nucleotides) of CDS-1 before the splice donor site codes for the amino acid Tryptophan in Frame +1, Arginine in Frame +2, and Valine in Frame +3.

The splice donor site at 19,150,858-19,150,859 is in phase 0 relative to Frame +1 because CDS-1 ends in a complete codon and thus has no extra nucleotides. In contrast, Frame +2 and Frame +3 don’t end in complete codons so they each have extra nucleotides. The splice donor site is in phase 2 relative to Frame +2 because there are 2 extra nucleotides (GG) after the last complete codon. The splice donor site is in phase 1 relative to Frame +3 since there is 1 extra nucleotide (G) after the last complete codon (Figure 62).

We previously determined CDS-1 is translated in Frame +3 (Figure 57); therefore, the last complete codon of CDS-1 (GTG which codes for Valine (V)) is located at 19,150,854-19,150,856 and there is one extra nucleotide (G at 19,150,857) between the last complete codon and the splice donor site. Hence, the CDS-1 splice donor site is in phase 1.

Our tblastn alignment in Part 5 placed CDS-2 at approximately 19,150,987 –​ 19,151,055 (Figure 55); therefore, we expect to find the splice acceptor site for CDS-2 around position 19,150,987.

There is only one potential canonical (“standard rule”) splice acceptor site (AG) in the 30 bp region surrounding the start of the tblastn alignment to CDS-2 (Figure 63, green box). The splice acceptor site is located at 19,150,983-19,150,984 (Figure 63, red box) and is supported by the NCBI RefSeq Genes and Spaln alignment of D. melanogaster proteins, the GeMoMa, N-SCAN, and Augustus gene predictions, and the aggregated RNA-Seq read coverage. Thus, CDS-2 starts at 19,150,985.

Now we need to determine the frame in which CDS-2 is translated.

Each frame can have extra nucleotides (0, 1, or 2) between the beginning of CDS-2 at 19,150,985 and the first complete codon of each frame (Figure 64).

Since the CDS-1 splice donor site at 19,150,858 –​ 19,150,859 is in phase 1 relative to Frame +3 (i.e., CDS-1 ends with one extra nucleotide), the CDS-2 splice acceptor site must be in phase 2 (i.e., CDS-2 must begin with two extra nucleotides) to maintain the ORF after Intron-1 has been removed. Since CDS-1 had one extra nucleotide (G) between the last complete codon (GTG = Valine; V) and the splice donor site, CDS-2 must start with two extra nucleotides to join the extra one from CDS-1 because we need 3 nucleotides to make a complete codon.

Since CDS-2 is in Frame +1 and the first complete codon (AAA codes for K) is located at 19,150,987- 19,150,989, there are two nucleotides (GC at 19,150,985- 19,150,986) between the potential splice acceptor site and the first complete codon. Hence, the CDS-2 splice acceptor site is in phase 2.

The extra nucleotides near the splice sites (i.e., G + GC) will form an additional amino acid (Glycine/G) after splicing (Figure 65). Collectively, our analysis suggests that CDS-1 ends at 19,150,857 with a phase 1 splice donor site and CDS-2 begins at 19,150,985 with a phase 2 splice acceptor site.

The same annotation strategy can be used to determine the phases for the remaining splice donor and splice acceptor sites between CDS-2 and CDS-3, CDS-3 and CDS-4, and CDS-4 and CDS-5 (Figure 66).

We reviewed the phases of splice donor and acceptor sites in Figure 60. Recall that to maintain the open reading frame (ORF) across adjacent CDS’s, the phases of the donor and acceptor sites of adjacent CDS’s must be compatible with each other (i.e., the sum of the donor and acceptor phases of adjacent CDS’s must either be 0 or 3).

Looking at the splice sites between CDS-1 and CDS-2, we see the splice donor phase of CDS-1 is one and the splice acceptor phase of CDS-2 is two. Thus, the sum of the donor and acceptor phases for Intron-1 is three (Figure 67).

There is only one splice junction predicted in this region.

The splice junction between the phase 1 donor site of CDS-1 and the phase 2 acceptor site of CDS-2 is supported by the NCBI RefSeq Genes and Spaln alignment of D. melanogaster proteins, gene predictors GeMoMa, N-SCAN, and Augustus, and the RNA-Seq data.

The score tells us how many spliced RNA-Seq reads there are that support a predicted splice junction. Since JUNC00113127 has a score of 6257, that splice junction prediction is supported by 6,257 spliced RNA-Seq reads. Splice junction predictions are color-coded based on the number of spliced RNA-Seq reads that support the junction (i.e., their scores). Based on the color-coded table in the “Description” section, JUNC00113127 will be red in the Genome Browser image since greater than 1,000 spliced RNA-Seq reads support the feature (Figure 70, bottom).

Based on multiple lines of evidence, we can conclude the splice junction prediction JUNC00113127 is consistent with our splice donor site for CDS-1 at 19,150,858-19,150,859 and our splice acceptor site for CDS-2 at 19,150,983-19,150,984 that we annotated in Part 6.3.

The same annotation strategy can be used to determine the coordinates for Intron-2 between CDS-2 and CDS-3.

There are two splice junctions predicted in this region, JUNC00113128 and JUNC00113129.

The tblastn alignment for CDS-2 ends at 19,151,055 (Figure 55). The potential splice donor site at 19,151,057-19,151,058 for CDS-2 is supported by the splice junction predictions JUNC00113128 and JUNC00113129, the NCBI RefSeq Genes and Spaln alignment of D. melanogaster proteins, gene predictors GeMoMa, N-SCAN, and Augustus, and the RNA-Seq data. There is one nucleotide (A at 19,151,056) between the last complete codon (AAC) and the potential splice donor site. Hence, this splice donor site is in phase 1 relative to Frame +1.

The tblastn alignment for CDS-3 spans from 19,151,156-19,151,359 in Frame +2 (Figure 55). The potential splice acceptor site at 19,151,152-19,151,153 for CDS-3 is supported by the splice junction prediction JUNC00113128, the NCBI RefSeq Genes and Spaln alignment of D. melanogaster proteins, gene predictors GeMoMa, N-SCAN, and Augustus, and the RNA-Seq data. There are two nucleotides (CC at 19,151,154-19,151,155) between the first complete codon (TTC) and the potential splice acceptor site. Hence, this splice acceptor site is in phase 2 relative to Frame +2. This phase 2 splice acceptor site is compatible with the phase 1 splice donor site for CDS-2.

In addition to the splice junction JUNC00113128, which supports the proposed splice acceptor site for CDS-3 at 19,151,152-19,151,153, there is another splice junction which suggests a different splice acceptor site (JUNC00113129) (Figure 71).

Thus, we need to investigate whether predicted junction JUNC00113129 is supported by other lines of evidence. CDS-3 in D. yakuba includes two methionine in Frame +2 (at 19,151,255-19,151,257 and 19,151,348-19,151,350). Hence, the splice junction JUNC00113129 could indicate the presence of a novel isoform of Rheb in D. yakuba. However, when we assess the number of spliced RNA-Seq reads that support the splice junction JUNC00113129, we find that this junction is weakly supported by only 8 spliced RNA-Seq reads; thus, there is little evidence to postulate a novel isoform of Rheb in D. yakuba based on this splice junction prediction.

Taking into account the splice site phases we found in Part 6.3 (Figure 66) and using Parts 6.4 and 6.5 as a guide, you will be able to find the start and end coordinates of the remaining CDSs for the gene model of Rheb-PA in D. yakuba (Figure 73).

Our analysis of the CDS-by-CDS tblastn alignments and the evidence tracks on the Genome Browser allowed us to precisely define the start and end positions of each of the five coding exons (CDS’s) of Rheb-PA. To verify that our proposed gene model satisfies the basic biological constraints (e.g., begins with a start codon, has compatible splice sites, and ends with a stop codon), we will check our gene model coordinates using the Gene Model Checker.

Once the analysis is complete, the right panel of the Gene Model Checker contains the results. The “Checklist” tab enumerates the list of criteria that have been checked by the Gene Model Checker (Figure 75). For example, the Gene Model Checker verifies that our proposed gene model begins with a start codon and ends with a stop codon. It also verifies that the splice junctions contain the canonical (“standard”) splice donor (GT) and acceptor (AG) sites. Some of the items on the checklist have been skipped because they do not apply to a complete gene (e.g., CDS-1 doesn’t have a splice acceptor site and CDS-5 doesn’t have a splice donor site).

In addition to verifying the basic gene structure, the Gene Model Checker also compares the proposed gene model against the putative D. melanogaster ortholog using a protein alignment and a dot plot (Figure 77).

The alternating color boxes in the dot plot correspond to the different CDS’s in the two sequences. Dots in the dot plot correspond to regions of similarity between the D. melanogaster protein and the submitted D. yakuba gene model. If the submitted sequence is identical to the D. melanogaster ortholog, then the dot plot will show a straight diagonal line with a slope of 1. Changes in the size of the submitted model compared to the D. melanogaster ortholog will alter the slope of this line.

In this case, the dot plot shows that the five CDS’s of Rheb-PA in D. melanogaster and D. yakuba have similar lengths (compare the length shown on the x-axis to the length shown on the y-axis for each CDS). However, within a small region of CDS-4, the dot plot did not detect sequence similarity between the submitted model for D. yakuba and the D. melanogaster ortholog (Figure 76).

To further investigate the dot plot, we will examine the protein alignment between the two sequences.

The protein alignment shows the comparison of the D. melanogaster protein against the conceptual translation for the submitted D. yakuba gene model. Like the dot plot, alternating colors correspond to the different CDS’s (Figure 78).

The protein alignment between the D. melanogaster ortholog and the D. yakuba gene model shows that the five CDS’s are 97.3% identical. The symbols in the match line denote the level of similarity (“*” indicates conserved amino acids, “:” denotes amino acids with highly similar chemical properties). Hence, the tiny gap in the dot plot within CDS-1 can be attributed to similar, but not identical, amino acids near its center.

The protein alignment between the D. melanogaster Rheb-PA and the submitted gene model in D. yakuba shows that that the gap within CDS-4 in the dot plot can be attributed to differences in three amino acid residues.

We can view the submitted gene model within the context of the other evidence tracks in the Genome Browser.

Now we see our entire submitted gene model for Rheb-PA in D. yakuba (Figure 80).

In addition to the Pathways Project: Annotation Form, you must prepare three additional data files to submit a project to the GEP — a General Feature Format File (GFF), a Transcript Sequence File (fasta), and a Peptide Sequence File (pep). The Gene Model Checker automatically creates these three files for a specific isoform (e.g., Rheb-PA) when you verify a gene model.

Recall that Rheb in D. yakuba has two isoforms; the files we just downloaded were for the Rheb-PA isoform. Since the coding exons (CDS’s) for both of our isoforms are identical, we don’t need to verify any additional models. However, if your project has multiple unique isoforms (i.e., don’t all have identical coding regions) then you should repeat Parts 7.1 and 7.2 for each unique isoform.

Now we need to combine the GFF, transcript sequence, and peptide sequence files for all unique isoforms into a single file prior to project submission (we’ll submit one merged file for each file type). We do NOT need to do this for Rheb in D. yakuba; however, we are including the steps for merging two separate files, so you’ll know how to do so. Again, if your project had only one unique isoform, you’d be done at this point. The following is merely to show you how to merge your files IF your project had multiple unique isoforms.

Let’s merge our two isoform GFF files first.

Now we need to merge our two isoform Transcript Sequence Files (.fasta).

Lastly, we need to merge our two isoform Peptide Sequence Files (.pep).