Hello all! I'm trying to utilize the prokaryotic profile. I got the .fna file for the bacterial strain of interest and created a gff file with Prokka (prokka --outdir prokka_RSM97 --prefix RSM97 --kingdom Bacteria --locustag RSM97 RSM97.fna). I ran the pipeline but wasn't able to get the biotype QC to work. I skipped that step and the pipeline completed nicely but the client is now asking for information about RNA types. Do you think the reason biotype QC had to be skipped is because of my gtf production method? Is there an alternate way to get that info for the strain? Also, many genes had identifiers given by the Prokka tool but didn't have a gene name / traditional ID. I'm also wondering if this is 'normal' for bacterial RNAseq or if it indicates an issue with my methodology. Some genes had names but the majority didn't. Any ideas? Thanks!

Hi! The prokaryotic profile is relatively new, and I'm not a prokaryotic expert myself, so I don't know how to answer your question. But the following is what Claude says having done a bit of research:

Biotype QC: That step (featureCounts) defaults to grouping by a gene_biotype attribute (controlled by --featurecounts_group_type, default gene_biotype). Prokka's GFF doesn't emit gene_biotype at all - its feature class lives in column 3 (CDS / tRNA / rRNA / tmRNA), and the attributes it writes
  are ID, locus_tag, product, gene (when matched), EC_number, inference, etc. The pipeline docs do mention this case briefly and recommend --skip_biotype_qc for non-Ensembl annotations, which is what you ended up doing.

  To get the RNA-type breakdown back, a couple of options:

  1. Patch the GFF. Add gene_biotype=... to each feature's attributes, derived from column 3 (CDS → protein_coding, tRNA → tRNA, rRNA → rRNA, tmRNA → tmRNA). Heads-up though: the prokaryotic profile also sets featurecounts_feature_type=CDS, so to actually see a breakdown you'd want to widen that
  (e.g. via --featurecounts_feature_type gene after also patching parent gene records into the GFF), otherwise featureCounts only looks at CDS rows.
  2. Use a RefSeq/NCBI annotation if your strain (or a very close one) is there - those tend to come with biotype-equivalent fields and integrate more cleanly.

  Gene names: Totally normal for bacteria. Prokka only assigns gene=... when there's a confident hit against its curated DBs (HAMAP and friends). For a non-model strain it's expected that the majority of CDSs come back with just a locus tag and no gene name. The product field is usually the more
  informative one. If you want richer functional annotation you can run eggNOG-mapper or InterProScan against the Prokka proteins and merge those results back in, but it won't change the biotype QC story.

I was not aware of Prokka (have not worked in microbiology for 15 years), but I looked it up, and apparently it is no longer maintained. They suggest to use its successor, Bacta, which may have more up-to-date databases that could help in your situation. The basic problem is that to use the

gene_biotype

statitistics in #CE8SSJV3N you need that field in the GTF, and it can only be there is the (basic) function of the gene is know. Since you are annotating from scratch (i.e. you not using a preannotated genome reference, presumably because there is none), it really depends on how well each gene matches any homologs in the databases from well-studied organisms, and how well the algorithm of your tool is able to detect that homology. More comprehensive databases should help here, but if your strain does not have well-characterized "relatives", the information you can attach to genes may be limited. (Additionally, the tool has to actually find the correct ORFs in the first place, before they can be compared to the database. Any mistakes there can further limit success. Might be worthwhile seeing if Bacta gives you different results.) The best you can do (beyond trying a different tool/different databases) is to BLAST the identified ORFs yourself and see if the output reveals anything beyond what the annotation tools give you (the tools will probably tend to be conservative to avoid being misleading when finding a tentative match with low confidence). That said, I am not very convinced that there is a lot of interesting information in the

gene_biotypes

. Bacteria have a lot fewer of those compared to eukaryotes, so besides the obvious ones, like rRNAs, tRNAs, tmRNAs, etc., I would expect most of them to just be "protein coding".

Thank you guys for your thorough responses! I'll try to spend some time working on adding annotations, perhaps using Bacta. I'll also spend some time thinking of the relevance of gene_biotypes analysis. A separate question is this: We found some differential genes that report a p-value of 0. I'll attach a snapshot of one example of this. What could cause a p-value to be reported as 0? Is it just a rounding thing? Is there a way to get a p-value for this gene so that we can report on it? It looks like a real difference but I doubt we could confidently report the change without a p-value.

What could cause a p-value to be reported as 0? Is it just a rounding thing?

I cannot tell from the picture if that number was was just lower than whatever precision they normally display (e.g. 3 or 4 digits after the decimal point) and that's why it looks like zero, if it is actually "zero" (if there is some text or Excel table with the full information, you can tell). I had the latter case in an experiment and PI had the same complaint as you. The reason for this is that the smallest number you can store in a 64-bit "double-precision" floating point number (which is commonly used for standard math in computer programs) is about 5*10^-324. If the p-value calculation (with rounding errors, etc.) yields a value lower than this, it will be rounded to and stored as zero.

Is there a way to get a p-value for this gene so that we can report on it?

If you don't want to show the "zero", you could say something like "<10^-300". Giving an upper bound is often done for p-values in scientific literature (most people don't care if a p-value is 10^-10, 10^-30, 10^-300, or less).

Perfect! I figured it was related to rounding and that the p-value is truly just that small. I'll relay that to the team. Thank you!

P-values are "the probability that rejecting the null hypothesis is wrong" given the test result for some given data. If you have a simple, mathematically tractable distribution, you can calculate that with some precision. But that is still based on the assumption that this distribution is actually true for the data (including the parameters that are estimated from the data), and often the math cannot be done precisely and one has to attempt a numeric estimate. (And that's before you even get into multiple-testing correction, etc.). It helps to think of p-value calculations as "voodoo" and not put too much stock into a precise value.