Loading Gene Clusters From FASTA Files
+ + + Source:vignettes/LoadFastaFiles.Rmd
+ LoadFastaFiles.RmdIntro +
+This tutorial demonstrates creating gene cluster visualizations from
+FASTA files using geneviewer. For guidance on importing
+clusters from GenBank files see: link.
+geneviewer parses FASTA headers to extract key-value pairs,
+denoted by =. For instance, in the below FASTA header it identifies
+‘protein’ as ‘putative phosphoenolpyruvate synthase’ and similarly
+processes other pairs. The ‘location’ tag is crucial because it allows
+the extraction of ‘start’, ‘end’, and ‘strand’ information for each
+gene.
>lcl|HQ3864.1_prot_ADX475.1_1 protein=putative phosphoenolpyruvate synthase protein_id=ADX66475.1 location=complement(<1..2247) gbkey=CDS
+MRATGLVRGRAAKRFGRAGGAGDAIGQCRATGHDCLGGSAMIEQYVWDLHEVDETQVAVVGGKGAHLGGL
+SRIEGIRVPAGFCVTTDAFRRIMAEAPSIDDGLDQLSRLNPDDREAIRTLSAQIRRTIEGIAIPGDLAAA
+ITRALARLGEHAACAVRSSATAEDLPTASFAGQQDTYLNVVGPTAILQHVSRCWASLFTERAVTYRQRNGMaterials +
+The FASTA files uses for the tutorial were originally downloaded from +the NCBI website but +can be downloaded directly from the geneviewer-tutorials +repository.
+Loading FASTA files +
+The example FASTA files used for this tutorial are each from a
+different organism and holds the genes of a cluster that produces a
+polyene antifungal agent. Download the files and place them under the
+same folder. We can read the files into R using the
+read_fasta function from geneviewer. The
+function uses the Biostrings library so we will need to
+load this as well. The function will parse all the FASTA headers and
+create a data.frame from the extracted information. We have set sequence
+to FALSE since we do not need the sequences to create the
+cluster visualization.
+library(geneviewer)
+# if (!require("BiocManager", quietly = TRUE))
+# install.packages("BiocManager")
+#
+# BiocManager::install("Biostrings")
+library(Biostrings)
+
+# change the path to the folder where you have saved the file
+file_path <- "~/path/to/folder/"
+fasta_df <- read_fasta(file_path, sequence = FALSE)
+View(fasta_df) # Inspect the data.frame in RstudioThe data.frame holds the protein IDs, start, end and strand +information for each gene. In addition there is the cluster column that +tells us from which file each gene is coming.
+| protein | +protein_id | +location | +gbkey | +gene | +start | +end | +strand | +cluster | +locus_tag | +
|---|---|---|---|---|---|---|---|---|---|
| putative | +ADX66475.1 | +complement(<1..2247) | +CDS | +NA | +1 | +2247 | +complement | +BGC0000108 | +NA | +
| putative | +ADX66476.1 | +complement(3301..3870) | +CDS | +sigD | +3301 | +3870 | +complement | +BGC0000108 | +NA | +
| conserved | +ADX66477.1 | +complement(4356..5282) | +CDS | +calU-c | +4356 | +5282 | +complement | +BGC0000108 | +NA | +
| ScnRII | +ADX66474.1 | +complement(5576..6154) | +CDS | +scnRII | +5576 | +6154 | +complement | +BGC0000108 | +NA | +
| ScnRI | +ADX66458.1 | +6555..10085 | +CDS | +scnRI | +6555 | +10085 | +forward | +BGC0000108 | +NA | +
| putative | +ADD85130.2 | +10250..11626 | +CDS | +scnK | +10250 | +11626 | +forward | +BGC0000108 | +NA | +
Cluster Visualization +
+Using the data.frame we created from the FASTA files we can now +easily make the cluster visualization.
+
+# Make the cluster chart
+chart <- GC_chart(
+ fasta_df,
+ start = "start",
+ end = "end",
+ strand = "strand",
+ cluster = "cluster",
+ height = "400px"
+) %>%
+ GC_clusterLabel(unique(fasta_df$cluster))
+chartUniprot Keywords +
+We can get additional information for each gene using the
+Uniprot.ws package. We make the connection to the database
+using the UniProt.ws function. Using
+keytypes(up) we can see all the keys we can use to query
+the database and with columns(up) we can see all the
+different fields we can get in return. Our gene IDs are
+EMBL-GenBank-DDBJ_CDS identifiers. For this tutorial we want to return
+the keywords associated with each gene such that we can use it to group
+the genes based on their function.
+library(UniProt.ws)
+
+up <- UniProt.ws()
+
+keytypes(up)
+columns(up)
+
+keywords <- UniProt.ws::select(
+ up,
+ column = c("keyword"),
+ keys = fasta_df$protein_id,
+ keytype = "EMBL-GenBank-DDBJ_CDS"
+ )This will return a data.frame with our original identifiers, +uniprotIDs and the keywords associated to each gene.
+| From | +Entry | +Keywords | +
|---|---|---|
| ADX66475.1 | +F1CLA2 | +Pyruvate | +
| ADX66476.1 | +F1CLA3 | +Transcription;Transcription regulation | +
| ADX66477.1 | +F1CLA4 | +NA | +
| ADX66474.1 | +B9V2R3 | +DNA-binding | +
| ADX66458.1 | +F1CLA6 | +DNA-binding;Two-component regulatory system | +
| ADD85130.2 | +D4P8Q9 | +Glycosyltransferase;Transferase | +
We now do a bit of data wrangling with the help of the
+dplyr package performing the following steps:
-
+
- We split the keywords for each gene on “;” +
- Count the occurrence of each keyword +
- For each gene keep only the highest occurring keyword +
- Replace any keyword that has less then 10 occurrences to +“Other” +
- Bind the keywords to our fasta_df by protein ID +
- Replace any genes without keyword to Other +
- Select only the columns we need for our visualization +
+library(dplyr)
+
+keywords_count <- keywords %>%
+ # Separate the Keywords into different rows
+ tidyr::separate_rows(Keywords, sep = ";") %>%
+ # Replace NA values with 'Other'
+ mutate(Keywords = replace(Keywords, is.na(Keywords), "Other")) %>%
+ # Add a Count column that counts each Keyword
+ add_count(Keywords)
+
+# Keep top ranking keywords for each protein ID
+keywords_count <- keywords_count %>%
+ group_by(From) %>%
+ filter(rank(-n, ties.method = "first") == 1) %>%
+ ungroup()
+
+# Replace any keyword that has a count less then 10 with "Other"
+keywords_count <- keywords_count %>%
+ mutate(Keywords = ifelse(n < 10, "Other", Keywords))
+
+# Bind Keywords back to fasta_df
+fasta_df_with_keywords <-
+ left_join(fasta_df, keywords_count, by = c("protein_id" = "From")) %>%
+ mutate(Keywords = replace(Keywords, is.na(Keywords), "Other")) %>%
+ select(cluster, start, end, strand, protein_id, protein, Keywords, n)| cluster | +start | +end | +strand | +protein_id | +protein | +Keywords | +n | +
|---|---|---|---|---|---|---|---|
| BGC0000108 | +1 | +2247 | +complement | +ADX66475.1 | +putative | +Other | +1 | +
| BGC0000108 | +3301 | +3870 | +complement | +ADX66476.1 | +putative | +Other | +3 | +
| BGC0000108 | +4356 | +5282 | +complement | +ADX66477.1 | +conserved | +Other | +13 | +
| BGC0000108 | +5576 | +6154 | +complement | +ADX66474.1 | +ScnRII | +DNA-binding | +23 | +
| BGC0000108 | +6555 | +10085 | +forward | +ADX66458.1 | +ScnRI | +DNA-binding | +23 | +
| BGC0000108 | +10250 | +11626 | +forward | +ADD85130.2 | +putative | +Transferase | +44 | +
We can now make our visual again but use the keywords to group the +genes. We alter the order of the legend such that it is in alphabetical +order and the Other group is placed last and finally we customize the +tooltip to also show the protein Id on hoover.
+
+
+key_legend <- c(setdiff(sort(unique(fasta_df_with_keywords$Keywords)), "Other"), "Other")
+
+# Make the cluster chart
+chart <- GC_chart(
+ fasta_df_with_keywords,
+ start = "start",
+ end = "end",
+ strand = "strand",
+ group = "Keywords",
+ cluster = "cluster",
+ height = "400px"
+) %>%
+ GC_clusterLabel(unique(fasta_df$cluster)) %>%
+ GC_legend(order = key_legend) %>%
+ GC_tooltip(
+ formatter = "
+ <b>ProteinID:</b> {protein_id}<br>
+ <b>Region:</b> {start} - {end}
+ "
+ )
+chartFor further customization options like setting a title, altering the +colors adding gene links and much more see Get Started and Examples.
+