0
$\begingroup$

I am doing some RNA analysis and am having issues trying to generate count data. I mapped my reads to a reference genome fasta file (genbank fasta file from ncbi) using bbmap and .sam files as the output.

I am now trying to use featurecounts to generate count data but am running into issues with my gtf file. The gtf file is also from ncbi and corresponds to the alignment fasta file used for mapping. This is it's structure:

#gtf-version 2.2
#!genome-build ASM14559v1
#!genome-build-accession NCBI_Assembly:GCA_000145595.1
CP002120.1      Genbank gene    1       1362    .       +       .       gene_id "SAA6008_00001"; transcript_id ""; gbkey "Gene"; gene "dnaA"; gene_biotype "protein_coding"; locus_tag "SAA6008_00001";
CP002120.1      Genbank CDS     1       1359    .       +       0       gene_id "SAA6008_00001"; transcript_id "unassigned_transcript_1"; gbkey "CDS"; gene "dnaA"; locus_tag "SAA6008_00001"; product "chromosomal replication initiator protein DnaA"; protein_id "ADL64089.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank start_codon     1       3       .       +       0       gene_id "SAA6008_00001"; transcript_id "unassigned_transcript_1"; gbkey "CDS"; gene "dnaA"; locus_tag "SAA6008_00001"; product "chromosomal replication initiator protein DnaA"; protein_id "ADL64089.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank stop_codon      1360    1362    .       +       0       gene_id "SAA6008_00001"; transcript_id "unassigned_transcript_1"; gbkey "CDS"; gene "dnaA"; locus_tag "SAA6008_00001"; product "chromosomal replication initiator protein DnaA"; protein_id "ADL64089.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank gene    1640    2773    .       +       .       gene_id "SAA6008_00002"; transcript_id ""; gbkey "Gene"; gene "dnaN"; gene_biotype "protein_coding"; locus_tag "SAA6008_00002";
CP002120.1      Genbank CDS     1640    2770    .       +       0       gene_id "SAA6008_00002"; transcript_id "unassigned_transcript_2"; gbkey "CDS"; gene "dnaN"; locus_tag "SAA6008_00002"; product "DNA-directed DNA polymerase beta subunit"; protein_id "ADL64090.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank start_codon     1640    1642    .       +       0       gene_id "SAA6008_00002"; transcript_id "unassigned_transcript_2"; gbkey "CDS"; gene "dnaN"; locus_tag "SAA6008_00002"; product "DNA-directed DNA polymerase beta subunit"; protein_id "ADL64090.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank stop_codon      2771    2773    .       +       0       gene_id "SAA6008_00002"; transcript_id "unassigned_transcript_2"; gbkey "CDS"; gene "dnaN"; locus_tag "SAA6008_00002"; product "DNA-directed DNA polymerase beta subunit"; protein_id "ADL64090.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank gene    3163    3399    .       +       .       gene_id "SAA6008_00003"; transcript_id ""; gbkey "Gene"; gene_biotype "protein_coding"; locus_tag "SAA6008_00003";
CP002120.1      Genbank CDS     3163    3396    .       +       0       gene_id "SAA6008_00003"; transcript_id "unassigned_transcript_3"; gbkey "CDS"; locus_tag "SAA6008_00003"; product "conserved hypothetical protein"; protein_id "ADL64091.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank start_codon     3163    3165    .       +       0       gene_id "SAA6008_00003"; transcript_id "unassigned_transcript_3"; gbkey "CDS"; locus_tag "SAA6008_00003"; product "conserved hypothetical protein"; protein_id "ADL64091.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank stop_codon      3397    3399    .       +       0       gene_id "SAA6008_00003"; transcript_id "unassigned_transcript_3"; gbkey "CDS"; locus_tag "SAA6008_00003"; product "conserved hypothetical protein"; protein_id "ADL64091.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank gene    3396    4508    .       +       .       gene_id "SAA6008_00004"; transcript_id ""; gbkey "Gene"; gene "recF"; gene_biotype "protein_coding"; locus_tag "SAA6008_00004";
CP002120.1      Genbank CDS     3396    4505    .       +       0       gene_id "SAA6008_00004"; transcript_id "unassigned_transcript_4"; gbkey "CDS"; gene "recF"; locus_tag "SAA6008_00004"; note "Required for DNA replication; binds preferentially to single-stranded, linear DNA"; product "DNA repair and genetic recombination protein"; protein_id "ADL64092.1"; t$
CP002120.1      Genbank start_codon     3396    3398    .       +       0       gene_id "SAA6008_00004"; transcript_id "unassigned_transcript_4"; gbkey "CDS"; gene "recF"; locus_tag "SAA6008_00004"; note "Required for DNA replication; binds preferentially to single-stranded, linear DNA"; product "DNA repair and genetic recombination protein"; protein_id "ADL640$
CP002120.1      Genbank stop_codon      4506    4508    .       +       0       gene_id "SAA6008_00004"; transcript_id "unassigned_transcript_4"; gbkey "CDS"; gene "recF"; locus_tag "SAA6008_00004"; note "Required for DNA replication; binds preferentially to single-stranded, linear DNA"; product "DNA repair and genetic recombination protein"; protein_id "ADL640$
CP002120.1      Genbank gene    4518    6452    .       +       .       gene_id "SAA6008_00005"; transcript_id ""; gbkey "Gene"; gene "gyrB"; gene_biotype "protein_coding"; locus_tag "SAA6008_00005";
CP002120.1      Genbank CDS     4518    6449    .       +       0       gene_id "SAA6008_00005"; transcript_id "unassigned_transcript_5"; gbkey "CDS"; gene "gyrB"; locus_tag "SAA6008_00005"; note "TIGRFAM: DNA gyrase, B subunit PFAM: DNA gyrase, subunit B domain protein; ATP-binding region, ATPase domain protein domain protein; TOPRIM

When I try to run featurecounts in R using this code:

featureCounts("s1_mapped.sam", annot.ext = "GCA_000145595.1_ASM14559v1_genomic.gtf", isGTFAnnotationFile = T)

I get the following error:

The chromosome name contains unexpected characters: "CP002120.1 Staphylococcus aureus subsp. aureus str. JKD6008, complete genome" (77 chars)
featureCounts has to stop running

FATAL Error: The program has to terminate and no counting file is generated.

Error in featureCounts("s1_mapped.sam", annot.ext = "GCA_000145595.1_ASM14559v1_genomic.gtf",  : 
  No counts were generated.

Do I need to format my gtf file in some way? I have also tried the code with GTF.featureType = "gene" but I get the same error.

$\endgroup$
2
  • 1
    $\begingroup$ The error specifies that the Chromosome name is not correct (I think the "," or the spaces are causing the problems). According to the GTF the chomosome is just called "CP002120.1" without the additional text. Having all this text is really quite unusual. Can you get the header of the sam file? samtools -H your_file.sam should do this or look at the first lines of the sam file in a text editor. The chromosome names start with @SQ $\endgroup$
    – PPK
    Dec 1, 2021 at 12:18
  • $\begingroup$ Also I just realized that the GTF you posted seems truncated. The line with "Genbank CDS 3396" ends with "$" as do the next two lines. This does cause the whole third line from the bottom to be quoted. If this is the case in your file and not just some copy paste error your GTF is damaged $\endgroup$
    – PPK
    Dec 1, 2021 at 12:25

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.