Question

feature count command error

0

Entering edit mode

18 months ago

Adyasha • 0

Hi everyone ,

I have ran this feature count command :

featureCounts  -T 4 -p -a gencode.v43.basic.annotation.gtf -o featurecount.txt *.bam

this gave me this error :

  Process BAM file UI_E2_sorted.bam.bam...                                    ||
||    Paired-end reads are included.                                          ||
||    The reads are assigned on the single-end mode.                          ||
||    Total alignments : 49308686                                             ||
||    Successfully assigned alignments : 0 (0.0%)                             ||
||    Running time : 0.26 minutes

this error coming for all bam files.

I am using gencode.v43.basic.annotation.gtf as the annotation file but the feature count command is giving zero assigned alignments .Please help on this , I am not sure if it's the right gtf annotation file , I got it from NCBI.

I am using human ref genome hg38 for my analysis of RNA seq data .

This is how my bam file header looks by using this command :

samtools_0.1.18 view RA1_E2_sorted.bam.bam | head
NB551648:44:HK5KLBGXG:3:23507:16037:10429       419     NC_000001.11    13158   1       76M     =       183744  150203  GAAGGGGATGCACTGTTGGGGAGGCAGCTGTAACTCAAAGCCTTAGCCTCTGTTCCCACGAAGGCAGGGCCATCAG    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEE  AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:76 YS:i:0  YT:Z:CP NH:i:5
NB551648:44:HK5KLBGXG:3:23507:16037:10429       419     NC_000001.11    13158   1       76M     =       13225   143     GAAGGGGATGCACTGTTGGGGAGGCAGCTGTAACTCAAAGCCTTAGCCTCTGTTCCCACGAAGGCAGGGCCATCAG    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEE  AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:76 YS:i:0  YT:Z:CP NH:i:5
NB551648:44:HK5KLBGXG:3:23507:16037:10429       339     NC_000001.11    13225   1       76M     =       13158   -143    GGCCATCAGGCACCAAAGGGATTCTGCCAGCATAGTGCTCCTGGACCAGTGATACACCCGGCACCCTGTCCTGGAC    EEEEEEEEAEEEAEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA  AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:76 YS:i:0  YT:Z:CP NH:i:5
NB551648:44:HK5KLBGXG:3:12601:18920:17906       355     NC_000001.11    14360   1       2S74M   =       14414   132     GTCATCCTGCACAGCTAGAGATCCTTTATTAAAAGCACACTGTTGGTTTCTGCTCAGTTCTTTATTGATTGGTGTG    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEE  AS:i:-2 ZS:i:-2 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:74 YS:i:0  YT:Z:CP NH:i:2
NB551648:44:HK5KLBGXG:1:23104:13420:13760       355     NC_000001.11    14362   1       2S74M   =       14403   119     CGTCCTGCACAGCTAGAGATCCTTTATTAAAAGCACACTGTTGGTTTCTGCTCAGTTCTTTATTGATTGGTGTGCC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  AS:i:-2 ZS:i:-2 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:74 YS:i:0  YT:Z:CP NH:i:2
NB551648:44:HK5KLBGXG:3:22601:17672:18350       355     NC_000001.11    14366   1       1S75M   =       14517   228     AGCACAGCTAGAGATCCTTTATTAAAAGCACACTGTTGGTTTCTGCTCAGTTCTTTATTGATTGGTGTGCCGTTTT    AAAAAEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  AS:i:-1 ZS:i:-1 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:75 YS:i:0  YT:Z:CP NH:i:5
NB551648:44:HK5KLBGXG:3:22601:17672:18350       355     NC_000001.11    14366   1       1S75M   =       185038  150290  AGCACAGCTAGAGATCCTTTATTAAAAGCACACTGTTGGTTTCTGCTCAGTTCTTTATTGATTGGTGTGCCGTTTT    AAAAAEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  AS:i:-1 ZS:i:-1 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:75 YS:i:0  YT:Z:CP NH:i:5
NB551648:44:HK5KLBGXG:4:23409:8298:7644 355     NC_000001.11    14367   0       76M     =       184930  150180  CACAGCTAGAGATCCTTTATTAAAAGCACACTGTTGGTTTCTGCTCAGTTCTTTATTGATTGGTGTGCCATTTTCT    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  AS:i:-5 ZS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:69G6       YS:i:-10        YT:Z:CP NH:i:5
NB551648:44:HK5KLBGXG:2:13208:13922:5858        355     NC_000001.11    14399   1       76M     =       14437   114     GTTGGTTTCTGCTCAGTTCTTTATTGATTGGTGTGCCGTTTTCTCTGGAAGCCTCTTAAGAACACAGTGGCGCAGG    A//AAEEEEEEAA/A/AEE<EEEEE6AE6//E/AA//E/EEAEEEEE/EEE/6E6E<AA<AEAE6A////EE/AAE  AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:76 YS:i:0  YT:Z:CP NH:i:5
NB551648:44:HK5KLBGXG:3:22508:3274:10922        99      NC_000001.11    14402   1       76M     =       14492   166     GGTTTCTGCTCAGTTCTTTATTGATTGGTGTGCCGTTTTCTCTGGAAGCCTCTTAAGAACACAGTGGCGCAGGCTG    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:76 YS:i:0  YT:Z:CP NH:i:5

and this is how my gtf file looks :

##description: evidence-based annotation of the human genome (GRCh38), version 43 (Ensembl 109)
##provider: GENCODE
##contact: gencode-help@ebi.ac.uk
##format: gtf
##date: 2022-11-29
chr1    HAVANA  gene    11869   14409   .       +       .       gene_id "ENSG00000290825.1"; gene_type "lncRNA"; gene_name "DDX11L2"; level 2; tag "overlaps_pseudogene";
chr1    HAVANA  transcript      11869   14409   .       +       .       gene_id "ENSG00000290825.1"; transcript_id "ENST00000456328.2"; gene_type "lncRNA"; gene_name "DDX11L2"; transcript_type "lncRNA"; transcript_name "DDX11L2-202"; level 2; transcript_support_level "1"; tag "basic"; tag "Ensembl_canonical"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    11869   12227   .       +       .       gene_id "ENSG00000290825.1"; transcript_id "ENST00000456328.2"; gene_type "lncRNA"; gene_name "DDX11L2"; transcript_type "lncRNA"; transcript_name "DDX11L2-202"; exon_number 1; exon_id "ENSE00002234944.1"; level 2; transcript_support_level "1"; tag "basic"; tag "Ensembl_canonical"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    12613   12721   .       +       .       gene_id "ENSG00000290825.1"; transcript_id "ENST00000456328.2"; gene_type "lncRNA"; gene_name "DDX11L2"; transcript_type "lncRNA"; transcript_name "DDX11L2-202"; exon_number 2; exon_id "ENSE00003582793.1"; level 2; transcript_support_level "1"; tag "basic"; tag "Ensembl_canonical"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    13221   14409   .       +       .       gene_id "ENSG00000290825.1"; transcript_id "ENST00000456328.2"; gene_type "lncRNA"; gene_name "DDX11L2"; transcript_type "lncRNA"; transcript_name "DDX11L2-202"; exon_number 3; exon_id "ENSE00002312635.1"; level 2; transcript_support_level "1"; tag "basic"; tag "Ensembl_canonical"; havana_transcript "OTTHUMT00000362751.1";

NGS featurecount RNA sequencing linux • 1.4k views

ADD COMMENT • link updated 18 months ago by GenoMax 147k • written 18 months ago by Adyasha • 0

GenoMax · Answer 1 · 2023-05-04

1

Entering edit mode

18 months ago

Carlo Yague 8.9k

Your annotation and reference genome are not compatible ! Looks like you have been mapping on the transcriptome (not the genome) with contigs names like this: NC_000001.11 , whereas you are trying to count features on a completely different annotation with contigs names chr1.

You can fix this by using compatible reference and annotations. Download them from the same source : UCSC, ensembl, or illumina iGenome for instance.

ADD COMMENT • link 18 months ago by Carlo Yague 8.9k

0

Entering edit mode

Hi, thank you so much for your reply, I am really struggling with this issue since a quite long time now .

Do you have any idea which annotation file will work for my bam .

ADD REPLY • link 18 months ago by Adyasha • 0

1

Entering edit mode

Since you appear to be using NCBI's genome you should find the corresponding GTF and GFF files here: https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/

ADD REPLY • link 18 months ago by GenoMax 147k

0

Entering edit mode

but still its giving this error

ERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
The specified gene identifier attribute is 'gene_id'
An example of attributes included in your GTF annotation is 'ID=exon-NR_046018.2-1;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2'.

ADD REPLY • link updated 18 months ago by GenoMax 147k • written 18 months ago by Adyasha • 0

0

Entering edit mode

Add option -g gene to your featureCounts command.

ADD REPLY • link 18 months ago by GenoMax 147k

0

Entering edit mode

Hi,

now I am using the compatible gtf file

that looks like this

##gff-version 3

#!gff-spec-version 1.21

#!processor NCBI annotwriter

#!genome-build GRCh38.p14

#!genome-build-accession NCBI_Assembly:GCF_000001405.40

#!annotation-date 03/15/2023

#!annotation-source NCBI RefSeq GCF_000001405.40-RS_2023_03

##sequence-region NC_000001.11 1 248956422

##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=9606

NC_000001.11    RefSeq  region  1       248956422       .       +       .       ID=NC_000001.11:1..248956422;Dbxref=taxon:9606;Name=1;chromosome=1;gbkey=Src;genome=chromosome;mol_type=genomic DNA
NC_000001.11    BestRefSeq      pseudogene      11874   14409   .       +       .       ID=gene-DDX11L1;Dbxref=GeneID:100287102,HGNC:HGNC:37102;Name=DDX11L1;description=DEAD/H-box helicase 11 like 1 (pseudoge
NC_000001.11    BestRefSeq      transcript      11874   14409   .       +       .       ID=rna-NR_046018.2;Parent=gene-DDX11L1;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;Name=NR_046018.2;gbke
NC_000001.11    BestRefSeq      exon    11874   12227   .       +       .       ID=exon-NR_046018.2-1;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX
NC_000001.11    BestRefSeq      exon    12613   12721   .       +       .       ID=exon-NR_046018.2-2;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX
NC_000001.11    BestRefSeq      exon    13221   14409   .       +       .       ID=exon-NR_046018.2-3;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX
NC_000001.11    BestRefSeq      pseudogene      14362   29370   .       -       .       ID=gene-WASH7P;Dbxref=GeneID:653635,HGNC:HGNC:38034;Name=WASH7P;description=WASP family homolog 7%2C pseudogene;gbkey=Ge
NC_000001.11    BestRefSeq      transcript      14362   29370   .       -       .       ID=rna-NR_024540.1;Parent=gene-WASH7P;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;Name=NR_024540.1;gbkey=mi
NC_000001.11    BestRefSeq      exon    29321   29370   .       -       .       ID=exon-NR_024540.1-1;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    24738   24891   .       -       .       ID=exon-NR_024540.1-2;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    18268   18366   .       -       .       ID=exon-NR_024540.1-3;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    17915   18061   .       -       .       ID=exon-NR_024540.1-4;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    17606   17742   .       -       .       ID=exon-NR_024540.1-5;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    17233   17368   .       -       .       ID=exon-NR_024540.1-6;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    16858   17055   .       -       .       ID=exon-NR_024540.1-7;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    16607   16765   .       -       .       ID=exon-NR_024540.1-8;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    15796   15947   .       -       .       ID=exon-NR_024540.1-9;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7P
NC_000001.11    BestRefSeq      exon    14970   15038   .       -       .       ID=exon-NR_024540.1-10;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7
NC_000001.11    BestRefSeq      exon    14362   14829   .       -       .       ID=exon-NR_024540.1-11;Parent=rna-NR_024540.1;Dbxref=GeneID:653635,GenBank:NR_024540.1,HGNC:HGNC:38034;gbkey=misc_RNA;gene=WASH7
NC_000001.11    BestRefSeq      gene    17369   17436   .       -       .       ID=gene-MIR6859-1;Dbxref=GeneID:102466751,HGNC:HGNC:50039,miRBase:MI0022705;Name=MIR6859-1;description=microRNA 6859-1;gbkey=Gen

ADD REPLY • link updated 18 months ago by GenoMax 147k • written 18 months ago by Adyasha • 0