6.1 fix the gene names

The gene names are too long for NCBI to accept. It is time to swap them out for the locus tags. Pull out the gene ID from each element, make a sub-table of genes only, and number them sequentially.

gff.base$ID <- unlist(lapply(gff.base$V9, scSplit))

At this stage, I generate locus tags and use them also as gene-ID

6.2 adding custom locus tags

locus tags are added to the gene elements:

ID=scaffold_0-gene-0.2274;locus_tag=GMORB2_0001;

I already built a locus tag map, but the goal is to make an easily reproducible method here.

library(stringr)
library(tidyr)
gff.gene <- subset(gff.base, V3 == "gene")
t <- gff.gene
t$short <- "locus_tag=GMORB2_"
t$number <- 1
t$shorter <- "ID=GMORB2"
for (n in (1:nrow(t))) {
  t[n,12] <- str_pad(n, 4, pad = "0")
}
t <- t %>% unite(locus, short, number, sep = "")
t$number <- 1

for (n in (1:nrow(t))) {
  t[n,13] <- str_pad(n, 4, pad = "0")
}
t <- t %>% unite (shorter, shorter, number, sep = "_")
locus.map <- data.frame(t$ID, t$shorter)
colnames(locus.map) <- c("gene","new.name")
t <- t %>% unite (V9, shorter, locus, sep = ";")

t$ID <- str_replace(t$ID, "ID", "note")
t <- t %>% unite(V9, V9, ID, sep=";")
gff.gene <- t

the object “locus.map” can be used to fix CDS and mRNA frames

we make partial gff with only mRNA lines. It will be easier moving forward if we have three separate gffs, one each for gene, CDS, and mRNA. We can amend them as needed and then merge and order after the fact.

gff.mrna <- gff.gene
gff.mrna$V3 <- "mRNA"
scSplit <- function(x){
  strsplit(x,";")[[1]][1]
}
gff.mrna$ID <- unlist(lapply(gff.mrna$V9, scSplit))
gff.mrna$ID <- str_replace(gff.mrna$ID, "ID=", "ID=mRNA:")
scSplit <- function(x){
  strsplit(x,";")[[1]][2]
}
gff.mrna$parent <- unlist(lapply(gff.mrna$V9, scSplit))
gff.mrna$parent <- str_replace(gff.mrna$parent, "locus_tag", "Parent")
gff.mrna <- gff.mrna %>% unite(V9, ID, parent, sep = ";")

next we make sure the lineage of the CDS are correct

gff.cds <- subset(gff.base, V3 == "CDS")
colnames(gff.cds)[10] <- "gene"
gff.cds$gene <- str_replace(gff.cds$gene, "CDS:", "")
library(dplyr)
gff.cds <- left_join(gff.cds, locus.map)
gff.cds$new.name <- str_replace(gff.cds$new.name, "ID=", "ID=CDS:")
gff.cds$parent <- str_replace(gff.cds$new.name, "ID=CDS:", "Parent=mRNA:")
gff.cds <- gff.cds %>% unite(V9, new.name, parent, sep = ";")
gff.cds$gene <- str_replace(gff.cds$gene, "ID=","")
gff.cds <- gff.cds[c(1:8,10,9)]

All of the lineages are now corrected and names have been revised. I have three separate data frames, gene, mRNA, CDS.

that’s it.

7.1 cleaning up illegal characters

described in https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md, several characters must be escaped when they are used in descriptions

* ; semicolon (%3B)

* = equals (%3D)

* & ampersand (%26)

* , comma (%2C)

#many odd things here were detected in the GenBank discrepancy report upon initial submission
cleanfun <- function(x) {
  x <- str_replace_all(x, ",", "%2C")
  x <- str_replace_all(x, ";", "%3B")
  x <- str_replace_all(x, "=", "%3D")
  x <- str_replace_all(x, "_", " ")
  x <- str_replace_all(x, "Inherit from COG: ", "")
  x <- str_replace_all(x, "FOG: ", "")
  x <- str_replace_all(x, "\\. ", "\\, ")
  x <- str_replace_all(x, "a.k.a. ", "")
  x <- str_replace_all(x, "the structure of ", "")
  x <- str_replace_all(x, " homolog", "-like protein")
  x <- str_replace_all(x, " homologue", "-like protein")
  x <- str_replace_all(x, " activity", "")
  x <- str_replace_all(x, " gene", "")
  x <- str_replace_all(x, "Partial ", "")
  str_replace_all(x, "&", "%26")
}
name.kegg$description <- cleanfun(name.kegg$description)
name.cog$description <- cleanfun(name.cog$description)
map.en$description <- cleanfun(map.en$description)

All three files are now in place. Now I will make a map with a gff-friendly orthology for each, a single frame for each system.

Again, it should look like this:

product=metaxin;inference=protein motif:KEGG:K17776

*note that there is JUNK in the eggnog map, should be manually curated, look for non-canonical EggNOG orthology numbers and remove those rows

#merge the ortholgoy and name mapping files together and make O+description gff format
#kegg
map.kegg <- left_join(map.kegg, name.kegg)
df <- map.kegg
df$dummy <- "product"
df <- df %>% unite(product, dummy, description, sep = "=")
map.kegg <- df
df <- map.kegg
df$dummy <- "inference=protein motif:KEGG:"
df <- df %>% unite(inference, dummy, KO, sep="")
map.kegg <- df
df <- map.kegg
df <- df %>% unite(orthology, product, inference, sep=";")
map.kegg <- df

#COG
df <- left_join(map.cog, name.cog)
df$dummy <- "product"
df <- df %>% unite(product, dummy, description, sep = "=")
df$dummy <- "inference=protein motif:COG:"
df <- df %>% unite(inference, dummy, COG, sep="")
df <- df %>% unite(orthology, product, inference, sep=";")
map.cog <- df

#EggNOG
df <- map.en
df$dummy <- "product"
df <- df %>% unite(product, dummy, description, sep = "=")
df$dummy <- "inference=protein motif:EggNOG:"
df <- df %>% unite(inference, dummy, NOG, sep="")
df <- df %>% unite(orthology, product, inference, sep=";")
map.en <- df

we now have our clean maps with

• one row for each map

• genes with no map excluded

• mapped-to values are gff compatible

7.2 hierarchical annotation

Now is the tricky part. I want to annotate in preference of KEGG > COG > EggNOG.

I can handle this with an “is.in” type of list check

t <- gff.cds
#create a dummy column to receive annotation information
t$annot <- ""

spSplit <- function(x){
  strsplit(x,";")[[1]][1]
}

#add kegg, this is the easy part
for (n in (1:nrow(t))) {
  if (t[n,10] %in% map.kegg$gene) {
    t[n,11] <- map.kegg[map.kegg$gene == t[n,10],2]
  } 
}
t2<-t
#add cog, requires two tests
for (n in (1:nrow(t2))) {
  if (t2[n,10] %in% map.cog$gene & t2[n,11] == "") {
    t2[n,11] <- map.cog[map.cog$gene == t2[n,10],2]
  }
}
t3 <- t2
#add eggnog using same strategy
for (n in (1:nrow(t3))) {
  if (t3[n,10] %in% map.en$gene & t3[n,11] == "") {
    t3[n,11] <- map.en[map.en$gene == t3[n,10],2]
  }
}

gff.cds <- t3

t <- gff.cds %>% unite(V9, V9, annot, sep=";")
t <- t[-10]
gff.cds <- t

9.1 initial revisions

“output_file.stats” provides an overview of suspicious issues with the validation. I manually reformatted the file output_file.stats so that I could read into here easily:

stats <- read.csv("output_file.stats",sep = "\t",header = F)
datatable(stats)

Most of these are not something I am going to worry about. However there are 2 types of problems I should investigate (they are pseudogenes. First, import the full *.val file and filter down to the important problems. This is mostly just a list of entries, but I can parse out important information

val <- read.csv("output_file.val", sep = "\t",header = F, stringsAsFactors = F)
spSplit <- function(x){
  strsplit(x,"\\[")[[1]][2]
}
val$code <- unlist(lapply(val$V1, spSplit))
spSplit <- function(x){
  strsplit(x,"\\]")[[1]][1]
}
val$code <- unlist(lapply(val$code, spSplit))
val$dummy <- ":"
library(tidyr)
val <- val %>% unite(code, code, dummy, sep = "")

#make sub-stats table only the guys I will investigate
stats <- subset(stats, V3 == "investigate")

#then use the substats table to filter the val table
val <- subset(val, code %in% stats$V1)
write.csv(val, "val.curated.csv")

The validation items now drop from 7800 to 28. Look into parsing it a bit, writing, and consider what to do next.

the only thing I altered was conversion of 14 genes into pseudogenes (manually)

9.2 run curated gff trhough tbl2asn

Run again through tbl2asn_gff and check the validation;

/home/robert/Tools/table2asn_GFF.Linux-3.10.0-1062.9.1.el7.x86_64-x86_64 \
-M n -J -c efswx -euk -a r10k -X E -augustus-fix \
-t template.sbt \
-gaps-min 10 -l paired-ends -j "[organism=Geosmithia morbida] [isolate=1262]" -i \
refgenome/G.morbida.refgenome.fa -f gmorb2_curate.gff -o output_file.sqn -Z

Validation looks fine.

9.3 Consider previous submission errors

The discrepancy report from NCBI on the previous submission was domianted by issues that this pipeline fixes (bad kegg names, data structure problems). The only thing I address now is the overlapping genes;

Remove

Scaffold 37 was flagged as pure contaminat. Remove completely prior to final submission.

There were also errors listed in the email they sent, some overlapping genes that I should remove.

FIND_OVERLAPPED_GENES: 4 genes completely overlapped by other genes
ORIG/output_file.sqn:Gene   GMORB2_6716 lcl|scaffold_6:816662-817629    GMORB2_6716
ORIG/output_file.sqn:Gene   GMORB2_6718 lcl|scaffold_6:c820276-817760   GMORB2_6718
ORIG/output_file.sqn:Gene   GMORB2_2885 lcl|scaffold_15:58872-60434 GMORB2_2885
ORIG/output_file.sqn:Gene   GMORB2_2887 lcl|scaffold_15:c61373-60757    GMORB2_2887

These four genes were removed manually from the gff

/home/robert/Tools/table2asn_GFF.Linux-3.10.0-1062.9.1.el7.x86_64-x86_64 \
-M n -J -c efswx -euk -a r10k -X E -augustus-fix \
-t template.sbt \
-gaps-min 10 -l paired-ends -j "[organism=Geosmithia morbida] [isolate=1262]" -i \
refgenome/G.morbida.refgenome.curated.fa -f gmorb2_curate.gff -o output_file.sqn -Z

Submitted this again on 2/19/2020, treating as a resubmission of previous.

9.4 1st submission results

Following submission (~20 minutes later) The Discrepancy.txt report is generated and must be manually interpreted, line by line!

• Many inappropriate names, have to be dealt with sequentially

• Parent for CDS uses the wrong separator! _ not ;. I fixed this in the pipeline.

• SUSPECT_PRODUCT_NAMES: 124 features contain underscore… these I can fix in the pipeline also during description cleanup.

• many other name errors like above; extensively revised the cleanup script (line 269) to deal with these.

Example messages I cannot interpret:

I cannot see anything wrong here:
CDS_HAS_NEW_EXCEPTION: 39 coding regions have new exceptions
ORIG/output_file.sqn:CDS    trans-aconitate 3-methyltransferase [lcl|scaffold_0:c801355-801338, c801333-800371] GMORB2_0260

I cannot see anything wrong here:
These inconsistencies do not hold up upon inspection.  Perhaps the actual stop is before the end of the gene? I could manually repair; this is gene-calling inconsistency
FEATURE_LOCATION_CONFLICT: 35 features have inconsistent gene locations.
FEATURE_LOCATION_CONFLICT: RNA feature location does not match gene location
ORIG/output_file.sqn:mRNA   mRNA    lcl|scaffold_0:c944231-943218   GMORB2_0310
ORIG/output_file.sqn:Gene   GMORB2_0310 lcl|scaffold_0:c>944231-943215  GMORB2_0310

these make no sense at all:
GENE_PARTIAL_CONFLICT: mRNA feature partialness conflicts with gene on both ends
ORIG/output_file.sqn:mRNA   mRNA    lcl|scaffold_0:c250922-249185   GMORB2_0074
ORIG/output_file.sqn:Gene   GMORB2_0074 lcl|scaffold_0:c>250922-<249185 GMORB2_0074

GENE_PARTIAL_CONFLICT: mRNA feature partialness conflicts with gene on 5' end
ORIG/output_file.sqn:mRNA   mRNA    lcl|scaffold_0:c52647-51893 GMORB2_0018
ORIG/output_file.sqn:Gene   GMORB2_0018 lcl|scaffold_0:c>52647-51893    GMORB2_0018

GENE_PARTIAL_CONFLICT: mRNA feature partialness conflicts with gene on 3' end
ORIG/output_file.sqn:mRNA   mRNA    lcl|scaffold_0:c30780-28384 GMORB2_0010
ORIG/output_file.sqn:Gene   GMORB2_0010 lcl|scaffold_0:c30780-<28384    GMORB2_0010