Hello, I am running the pipeline with a reference genome that was automatically annotated using RAST annotation server (fna and gff files). The pipeline runs well until FeatureCount where I got this error. Any idea how to go through it?

The pipeline expects gff(3) format. Actually, it's the featureCount step that is sensitive and it seems to have noticed that you passed a gtf file but failed anyway. We do not test the pipeline with other formats so my guess this is why it failed. Maybe you can try to convert your file. There's a discussion here: https://www.biostars.org/p/9545464/

Sorry, will do next time!

It looks quite good in principle. The key things are there: the coordinates (from, to, strand) and the ID. Does the gff file contain only CDS features? By default, the pipeline counts CDS, rRNA, tRNA and tmRNA features. If any of those types are missing from the gff file, you might get an error. Try specifying

--features CDS

as a parameter.

In the file there is also tRNA and RNA features, and the RNA features are rRNA. Maybe I have to rename the RNA to rRNA?

It should work without renaming if you specify:

--features CDS,RNA,tRNA

.

Ok, and I should put it in the command line

nextflow run ... --feaures CDS,RNA,tRNA -params-file ...

or in a nextflow config file?

Either on the command line, or in a params yml or json file, not in a config file.

ok, I will try and let you know. Thanks!

We also get the

Unknown annotation format

message in our tests. Try dropping

RNA

from your features, keeping only

CDS

, to see if this works.

I also tried to convert my .gff to a prokka style gff with

agat_convert_sp_gxf2gxf.pl

as suggested earlier but I got a bunch of exon and it seems weird since it's a bacteria

That was likely not the problem. I got mislead by the warning message, until I realized we also get that.

Having

--features CDS

only seems to work:

Right, then we at least know where to look, i.e. the format of your file for

RNA

features.

As

features

we have CDS, RNA and tRNA. The RNA rows look like:

Scaffold_3	FIG	RNA	566940	568416	.	-	.	ID=fig|6666666.1056986.rna.31;Name=SSU rRNA
Scaffold_3	FIG	RNA	566940	568416	.	-	.	ID=fig|6666666.1056986.rna.31;Name=SSU rRNA
Scaffold_3	FIG	RNA	566940	568416	.	-	.	ID=fig|6666666.1056986.rna.31;Name=SSU rRNA

Looks like duplicates. I have no idea if that's a problem.

I just noticed that all

tRNA

and

RNA

rows looks like this, always 3 identical lines

Possibly. I'd also try removing the triplicate lines.

So triplicate lines doesn't seem to be an issue as with

--features CDS,tRNA

the pipeline was successful. I renamed the

RNA

into

rRNA

, do you think I should remove the triplicate or should I keep the file like the original RAST output?

I'm not sure. I don't know what happens when you have triplicates and I have no idea what they're for. Since you're anyway updating the file, I'd make sure entries are unique.

Even by renaming the

RNA

into

rRNA

I end up with the same error

And the triplicates? Strange. To investigate further, you can

cd

into the work directory (typically

work/nn/lotsofletters/

, visible in the log). There's a file called

.command.sh

that contains the actual command run. If you have Subread installed, you should be able to run the command, fiddle with featureCount params and see if you can get it to work.

this is what's there is in the

.command.sh

:

#!/usr/bin/env bash -C -e -u -o pipefail
featureCounts \
    -g ID -t rRNA -F gtf -M   \
    -p \
    -T 6 \
    -a genomes.gz \
    -s 0 \
    -o p96_4_cu.rRNA.featureCounts.tsv \
    p96_4_cu.sorted.bam

# capture process environment
set +u
set +e
cd "$NXF_TASK_WORKDIR"

nxf_eval_cmd() {
    {
        IFS=$'\n' read -r -d '' "${1}";
        IFS=$'\n' read -r -d '' "${2}";
        (IFS=$'\n' read -r -d '' _ERRNO_; return ${_ERRNO_});
    } < <((printf '\0%s\0%d\0' "$(((({ shift 2; "${@}"; echo "${?}" 1>&3-; } | tr -d '\0' 1>&4-) 4>&2- 2>&1- | tr -d '\0' 1>&4-) 3>&1- | exit "$(cat)") 4>&1-)" "${?}" 1>&2) 2>&1)
}

echo '' > .command.env
#
nxf_eval_cmd STDOUT STDERR bash -c "featureCounts -v 2>&1 | sed 's/featureCounts v//'"
status=$?
if [ $status -eq 0 ]; then
  echo nxf_out_eval_23="$STDOUT" >> .command.env
  echo /nxf_out_eval_23/=exit:0 >> .command.env
else
  echo nxf_out_eval_23="$STDERR" >> .command.env
  echo /nxf_out_eval_23/=exit:$status >> .command.env
fi

I'm installing subread now. How could I run this command using it?

You just need to run the top part:

featureCounts \
    -g ID -t rRNA -F gtf -M   \
    -p \
    -T 6 \
    -a genomes.gz \
    -s 0 \
    -o p96_4_cu.rRNA.featureCounts.tsv \
    p96_4_cu.sorted.bam

(That could be specified on one line; if you do, remove the

\

.)

and I should run this in the work folder where the

.command.sh

is right?

Correct.

So I get the same error :

ERROR: no features were loaded in format GTF. The annotation format can be specified by the '-F' option, and the required feature type can be specified by the '-t' option..

But when I run the same command with

-t tRNA

it is running, so it's really the

rRNA

part that causes the issue

This is really strange. I can only see one difference in the format for

RNA

entries compared to `CDS`: The 7th column are numbers for the latter, but

.

for the former. Prokka (the tool the pipeline uses) has the same format so it shouldn't be that. This is how my Prokka files look:

(prokka) dl@sequoia:/data/gtdb/genomes$ zgrep rRNA GCF_029167045.1_ASM2916704v1_genomic.prokka.d/GCF_029167045.1_ASM2916704v1_genomic.gff.gz |head
gnl|Prokka|CJHEBLHE_22  Prodigal:002006 CDS     9052    9837    .       +       0       ID=CJHEBLHE_00472;Parent=CJHEBLHE_00472_gene;eC_number=2.1.1.-;db_xref=COG:COG0565;inference=ab initio prediction:Prodigal:002006,similar to AA sequence:UniProtKB:P44676;locus_tag=CJHEBLHE_00472;product=putative tRNA/rRNA methyltransferase;protein_id=gnl|Prokka|CJHEBLHE_00472
gnl|Prokka|CJHEBLHE_66  barrnap:0.9     rRNA    88753   88869   1.9e-14 -       .       ID=CJHEBLHE_01095;Parent=CJHEBLHE_01095_gene;locus_tag=CJHEBLHE_01095;product=5S ribosomal RNA
gnl|Prokka|CJHEBLHE_66  barrnap:0.9     rRNA    89034   91925   0       -       .       ID=CJHEBLHE_01096;Parent=CJHEBLHE_01096_gene;locus_tag=CJHEBLHE_01096;product=23S ribosomal RNA
gnl|Prokka|CJHEBLHE_66  barrnap:0.9     rRNA    92274   93742   3.9e-226        -       .       ID=CJHEBLHE_01098;Parent=CJHEBLHE_01098_gene;locus_tag=CJHEBLHE_01098;product=16S ribosomal RNA
(prokka) dl@sequoia:/data/gtdb/genomes$ zgrep CDS GCF_029167045.1_ASM2916704v1_genomic.prokka.d/GCF_029167045.1_ASM2916704v1_genomic.gff.gz |head -n 4
gnl|Prokka|CJHEBLHE_1   Prodigal:002006 CDS     47      988     .       +       0       ID=CJHEBLHE_00001;Parent=CJHEBLHE_00001_gene;Name=scpB;gene=scpB;inference=ab initio prediction:Prodigal:002006,protein motif:HAMAP:MF_01804;locus_tag=CJHEBLHE_00001;product=Segregation and condensation protein B;protein_id=gnl|Prokka|CJHEBLHE_00001
gnl|Prokka|CJHEBLHE_1   Prodigal:002006 CDS     1022    1417    .       +       0       ID=CJHEBLHE_00002;Parent=CJHEBLHE_00002_gene;inference=ab initio prediction:Prodigal:002006;locus_tag=CJHEBLHE_00002;product=hypothetical protein;protein_id=gnl|Prokka|CJHEBLHE_00002
gnl|Prokka|CJHEBLHE_1   Prodigal:002006 CDS     1454    2296    .       +       0       ID=CJHEBLHE_00003;Parent=CJHEBLHE_00003_gene;inference=ab initio prediction:Prodigal:002006;locus_tag=CJHEBLHE_00003;product=hypothetical protein;protein_id=gnl|Prokka|CJHEBLHE_00003
gnl|Prokka|CJHEBLHE_1   Prodigal:002006 CDS     2443    3459    .       +       0       ID=CJHEBLHE_00004;Parent=CJHEBLHE_00004_gene;inference=ab initio prediction:Prodigal:002006;locus_tag=CJHEBLHE_00004;product=hypothetical protein;protein_id=gnl|Prokka|CJHEBLHE_00004

The RAST gff is very similar but it's not in the 7th column but in the 8th. Here from the original file before renaming into

rRNA

Scaffold_8	FIG	CDS	50383	52389	.	-	1	ID=fig|6666666.1056986.peg.4348;Name=Protein translocase subunit SecA
Scaffold_3	FIG	tRNA	562290	562366	.	-	.	ID=fig|6666666.1056986.rna.26;Name=tRNA-Met-CAT
Scaffold_3	FIG	RNA	562420	562541	.	-	.	ID=fig|6666666.1056986.rna.27;Name=5S rRNA

My counting was off.

Yes I will use mostly the CDSs as it is part of a project where they also did proteomics. rRNA depletion has been done on the samples so it should not be a lot attributed to them, but I guess it would be nice to have all the info

If depletion was done, the rRNA is useless I'm afraid. So, rely on the pipeline run using only CDS as feature.

Ok, I wasn't sure if not having the rRNA would affect the rest of the data

No. What you typically do is calculate TPMs over the sample sums for the feature of interest, usually CDS.

That's what I did with my other project. I just looked at the

magmap.overall.stats

and it the number of mapped index is only half of the

n_post_trimming

, but mostly

CDS

. Do you think it's mostly due because the genome annotated by RAST might not be a reference genome and not so complete?

sample  n_post_trimming idxs_n_mapped   idxs_n_unmapped CDS     tRNA
p96_1_ctr       19147680        9181044 916988  8908731 261
p96_1_cu        14965590        9470574 875426  9209124 292
p96_2_ctr       31405522        9963401 769139  9770194 260
p96_2_cu        19838230        10837155        985873  10550935        374
p96_3_ctr       18940164        9334830 1010148 9067306 240
p96_3_cu        16695126        9540558 882514  9289162 329
p96_4_ctr       17218266        9475844 979006  9198694 237
p96_4_cu        23705058        13001861        1198417 12674375        476

You should compare the

idxs_n_mapped

with

CDS

. It might be that

n_post_trimming

counts all reads, but that the other two counts pairs, i.e. half as many. The other alternative is that less than half actually map to the contigs. I need to investigate this... If I were you, I wouldn't worry about this.

Ok. Thank you very much for all the help :)

Glad you're using our pipeline! 🙂