How to deal with NA data for non model species?

[cagenet]

Hi Federico,
Thank you for this really nice and helpful package. It is really useful.
I used it for agronomical species (pig and bovine).
I encountered a problem when trying to use enhance_table function. I have the following error :

p <- enhance_table(res_enrich_pig,

•                res,
  

•                n_gs = 30,
  

•                gtl<-gtl_pig,
  

•                annotation_obj = resannot_Pig,
  

•                chars_limit = 60)
  

Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, :
les noms de lignes contiennent des valeurs manquantes
I guess the trouble came from the fact that I use ENSEMBLID, convert it in HGNC symbol but as not all genes are identified , I have some missing data (valeurs manquantes in french). Could it be possible to allow missing values in this function ?
Thank you in advance
Regards
Carine

[federicomarini]

Hi Carine!
Quick note on the way you call the function - it should be enough to use the gtl object, as the res/res_enrich/annotation are all picked from that.
Re: the missing values problem - does the error persist if you

• keep the ids as ensembl in the result
• have the conversion handled internally with the annotation object?

[cagenet]

Hi Federico,
Thank you so much for your quick reply.
First, I didn't catch what you mean by "the way I call the function ". I follow your user guide (step 5.1).
For the missing values problem :
I think that there is two issues : the first one is that for the pig specie there is no mapping between Ensembl and symbol in org.Ss.eg.db (see https://support.bioconductor.org/p/9137574/)...
So I used bovine species (I have 2 datasets I applied a treatment in bovine and in pig and want to compare the results obtained in the two species).
Now with the bovine species, the enhance_table work with symbol but not with ENSEMBL ids.
I do have a conversion handled internally in the annotation object (a two table file with gene_id and gene_symbol (call HGNC in my code) and I have in the RowData in dds object an annotation with Ensemblid and gene symbol (call HGNC in my code see below).
What I don't understand is that I have an conversion table in the annotation object and also in the RowData (dds).

Thank you for your help.
Cheers

Here is my code :

library(GeneTonic)
load(file="DESEQ2_dds1_Bov.rda")#dds2

dds2
class: DESeqDataSet
dim: 15779 20
metadata(1): version
assays(4): counts mu H cooks
rownames(15779): ENSBTAG00000000005 ENSBTAG00000000009 ... ENSBTAG00000055312 ENSBTAG00000055315
rowData names(60): ENSEMBLID HGNC ... deviance maxCooks
colnames(20): BE1_CTRL BE2_CTRL ... BE14_BMP15 BE16_BMP15
colData names(3): Cow Treatment sizeFactor

load(file="resBta_04juin.rda")

summary(resBta)
out of 15779 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up) : 2323, 15%
LFC < 0 (down) : 2460, 16%
outliers [1] : 0, 0%
low counts [2] : 306, 1.9%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

resannot_Bov<-read.table(file="annot_bov.txt", header = TRUE, sep="\t")

head(resannot_Bov)
gene_id gene_name
1 ENSBTAG00000000005 GRK3
2 ENSBTAG00000000008 KCNJ1
3 ENSBTAG00000000009 FOXF1
4 ENSBTAG00000000010 UBL7
5 ENSBTAG00000000011 TDH
6 ENSBTAG00000000012 TTC33

de_symbol_bov <- deseqresult2df(resBta, FDR = 0.05)$SYMBOL
bg_bov_sym <- rowData(dds2)$HGNC[rowSums(counts(dds2)) > 0]

resBta$ensemblid<-(row.names(resBta))
bg_ensembl_bov <- rowData(dds2)$ENSEMBLID[rowSums(counts(dds2)) > 0]
de_ensembl_bov <- deseqresult2df(resBta, FDR = 0.05)$ensemblid
library("topGO")
topgoBP_bov <- pcaExplorer::topGOtable(de_symbol_bov, bg_bov_sym,
ontology = "BP",
mapping = "org.Bt.eg.db",
geneID = "symbol",
topTablerows = 500)

topgoBPensembl_bov <- pcaExplorer::topGOtable(de_ensembl_bov, bg_ensembl_bov,
ontology = "BP",
mapping = "org.Bt.eg.db",
geneID = "ENSEMBL",
topTablerows = 500)
res_enrich_bov<-shake_topGOtableResult(topgoBP_bov)
res_enrich_bov_ensembl<-shake_topGOtableResult(topgoBPensembl_bov)
res_enrich_bov <- get_aggrscores(res_enrich = res_enrich_bov, res_de = resBta,
annotation_obj = resannot_Bov,
aggrfun = mean)
res_enrich_bov_ensembl<-get_aggrscores(res_enrich = res_enrich_bov_ensembl, res_de = resBta,annotation_obj = resannot_Bov,aggrfun = mean)
gtl_bov <- GeneTonic_list(
dds = dds2,
res_de = resBta,
res_enrich = res_enrich_bov,
annotation_obj = resannot_Bov
)

gtl_bov_ensembl<-GeneTonic_list(dds=dds2, res_de = resBta, res_enrich = res_enrich_bov_ensembl, annotation_obj = resannot_Bov)
p <- enhance_table(res_enrich_bov,
resBta,
n_gs = 30,
gtl<-gtl_bov,
annotation_obj = resannot_Bov,
chars_limit = 60)
p1<-enhance_table(res_enrich_bov_ensembl, resBta, n_gs=30, gtl = gtl_bov_ensembl, annotation_obj = resannot_Bov, chars_limit = 60)
Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, :
les noms de lignes contiennent des valeurs manquantes

[federicomarini]

Hi Federico,
Thank you so much for your quick reply.
First, I didn't catch what you mean by "the way I call the function ". I follow your user guide (step 5.1).

Hi Carine,
in theory you can call

enhance_table(gtl = gtl_bov_ensembl, n_gs=30, chars_limit = 60)

as gtl should have all it needs 😉

For the missing values problem :
I think that there is two issues : the first one is that for the pig specie there is no mapping between Ensembl and symbol in org.Ss.eg.db (see https://support.bioconductor.org/p/9137574/)...
So I used bovine species (I have 2 datasets I applied a treatment in bovine and in pig and want to compare the results obtained in the two species).
Now with the bovine species, the enhance_table work with symbol but not with ENSEMBL ids.
I do have a conversion handled internally in the annotation object (a two table file with gene_id and gene_symbol (call HGNC in my code) and I have in the RowData in dds object an annotation with Ensemblid and gene symbol (call HGNC in my code see below).
What I don't understand is that I have an conversion table in the annotation object and also in the RowData (dds).

Thank you for your help.
Cheers

Here is my code :

library(GeneTonic)
load(file="DESEQ2_dds1_Bov.rda")#dds2

dds2
class: DESeqDataSet
dim: 15779 20
metadata(1): version
assays(4): counts mu H cooks
rownames(15779): ENSBTAG00000000005 ENSBTAG00000000009 ... ENSBTAG00000055312 ENSBTAG00000055315
rowData names(60): ENSEMBLID HGNC ... deviance maxCooks
colnames(20): BE1_CTRL BE2_CTRL ... BE14_BMP15 BE16_BMP15
colData names(3): Cow Treatment sizeFactor

load(file="resBta_04juin.rda")

summary(resBta)
out of 15779 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up) : 2323, 15%
LFC < 0 (down) : 2460, 16%
outliers [1] : 0, 0%
low counts [2] : 306, 1.9%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

resannot_Bov<-read.table(file="annot_bov.txt", header = TRUE, sep="\t")

head(resannot_Bov)
gene_id gene_name
1 ENSBTAG00000000005 GRK3
2 ENSBTAG00000000008 KCNJ1
3 ENSBTAG00000000009 FOXF1
4 ENSBTAG00000000010 UBL7
5 ENSBTAG00000000011 TDH
6 ENSBTAG00000000012 TTC33

de_symbol_bov <- deseqresult2df(resBta, FDR = 0.05)$SYMBOL
bg_bov_sym <- rowData(dds2)$HGNC[rowSums(counts(dds2)) > 0]

resBta$ensemblid<-(row.names(resBta))
bg_ensembl_bov <- rowData(dds2)$ENSEMBLID[rowSums(counts(dds2)) > 0]
de_ensembl_bov <- deseqresult2df(resBta, FDR = 0.05)$ensemblid
library("topGO")
topgoBP_bov <- pcaExplorer::topGOtable(de_symbol_bov, bg_bov_sym,
ontology = "BP",
mapping = "org.Bt.eg.db",
geneID = "symbol",
topTablerows = 500)

Ok, now that you elaborate a little more on this (and since there's no ensembl->symbol conversion) it can be that the gs_genes column in the res_enrich...

topgoBPensembl_bov <- pcaExplorer::topGOtable(de_ensembl_bov, bg_ensembl_bov,
ontology = "BP",
mapping = "org.Bt.eg.db",
geneID = "ENSEMBL",
topTablerows = 500)
res_enrich_bov<-shake_topGOtableResult(topgoBP_bov)
res_enrich_bov_ensembl<-shake_topGOtableResult(topgoBPensembl_bov)

... here ☝️
does not have the gene symbols "as expected".

I am not so familiar with species other than human, mouse & fly, as these are the almost entirety of the catalogue our cooperation partners work with.

If you want to share with me the gtl object, I can have a look (will of course be handled confidentially).

Federico

res_enrich_bov <- get_aggrscores(res_enrich = res_enrich_bov, res_de = resBta,
annotation_obj = resannot_Bov,
aggrfun = mean)
res_enrich_bov_ensembl<-get_aggrscores(res_enrich = res_enrich_bov_ensembl, res_de = resBta,annotation_obj = resannot_Bov,aggrfun = mean)
gtl_bov <- GeneTonic_list(
dds = dds2,
res_de = resBta,
res_enrich = res_enrich_bov,
annotation_obj = resannot_Bov
)

gtl_bov_ensembl<-GeneTonic_list(dds=dds2, res_de = resBta, res_enrich = res_enrich_bov_ensembl, annotation_obj = resannot_Bov)
p <- enhance_table(res_enrich_bov,
resBta,
n_gs = 30,
gtl<-gtl_bov,
annotation_obj = resannot_Bov,
chars_limit = 60)
p1<-enhance_table(res_enrich_bov_ensembl, resBta, n_gs=30, gtl = gtl_bov_ensembl, annotation_obj = resannot_Bov, chars_limit = 60)
Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, :
les noms de lignes contiennent des valeurs manquantes

[cagenet]

Hi ;-)
First issue: I guess that it was something like that BUT I still obtained the same error ...
p<-enhance_table(gtl=gtl_bov_ensembl, n_gs = 30, chars_limit = 60)
Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, :
les noms de lignes contiennent des valeurs manquantes

Second issue : Thank you so much for your help ... Of course I can send you my gtl object ...
Will try to import and give you access to gtf object in my repository (GitHub newbie and never used it !)
Thank you

[federicomarini]

Sure, no prob.
You can reach me at the email address in the maintainer field of the package (trying to make bots' work to scrape my email a little harder, sorry for the extra step 😛 )