[iris_aurelia]

Hi Wei,

I have tested your test-dataset and it indeed gave the same results for both HTseq-count as FeatureCounts.
I think the difference is made by the NH tag in my BAM/SAM file. When I tried a single read that mapped 5 times it gives a different result.
According the documentation of HTseq-count it excludes all reads with more than one reported alignment (whenever the NH tag is available), which I believe FeatureCounts does not.
Using the mapping quality filter in FeatureCounts, it will produce the same result, since reads with more than one reported alignment have a lower mapping quality.

I have used Tophat for the alignment of my reads, followed by HTseq-count as well as FeatureCounts. I've attached a test set with 5 multimapped hits from the data I got from TopHat.
The commands I used were:

HTseq-count:

htseq-count MMread.sam annotation.gtf -s no -q > genecounts_MMread.txt

FeatureCounts:

featureCounts -a annotation.gtf -t exon -g 'gene_id' -i MMread.sam -T 5 -o feature_counts_MMread.txt

Thanks for the quick responses!

Iris

[Nicolas Nalpas]

Dear Wei,
This is great, thanks a lot for being so helpful and available.
I will give this a try on our dataset.
Thanks a lot.
Regards,
Nicolas

[shi]

Dear Nicolas,

We have just modified featureCounts to make it be able to perform reverse-stranded read counting. The latest version can be found in Subread v1.3.3-p5 (http://subread.sourceforge.net).

Providing a parameter '-s 2' to featureCounts will instruct it to perform reverse-stranded read counting.

Thank you for your suggestion and I hope it works for your dataset.

Best wishes,

Wei

[Nicolas Nalpas]

Dear Wei,

I just had another question about featureCounts and the strand specificity parameter.
We have a strand-specific RNA-seq dataset, on which we initially used HTseq-count, so as you know HTseq-count has three options for the strand-specific parameter: --stranded no or --stranded yes or --stranded reverse and in our case the appropriate options is actually the latter: --stranded reverse. So now we want to use featureCounts (as we found it very quick and accurate) on this dataset, however the current standed options are: unstranded by default or -s for stranded. So after testing we understood that featureCounts unstranded correspond to HTseq-count --stranded no and that featureCounts -s correspond to HTseq-count --stranded yes. So as I said above the option that we want for this dataset is really --stranded reverse, which we were not able to find in featureCounts.
So I was wondering if I missed something in the featureCounts users guide and this option is already available? Or if it is not possible to do that, I was wondering if you have any plan to incorporate such option?

Thanks a lot for your help.
Regards,
Nicolas

[shi]

Dear Iris,

Thanks for your clarification and providing the comparison results.

However, we could not reproduce what you have observed. Attached is a toy example. In this example, a read was mapped to five different locations and all locations were reported in the SAM file. We found that featureCounts and HTseq-count gave the identical results when summarizing it.

Could you please provide the commands you used when running featureCounts and HTseq-count? It will also be very helpful if you can provide a runnable example dataset so that we can try to figure out what's going on.

Thanks,

Wei

[iris_aurelia]

Hi Wei,

I do mean the reads that are mapped to different chromosomal locations. So for example if a read maps to gene A on chr8 and maps to gene F on chr 9 the read is counted for both genes.
I've indeed read about the open and closed position what could result in a difference between the two tools, but that's not the case here.

I downloaded the newest version, subread-1.3.3-p4 and used the quality filter option to see if this will make a difference. When I only used the reads with the highest mapping qualities featureCounts reported exactly the same counts as HTseq-count does. (see attached figure). HTseq-count hower never reported more counts than featureCounts in my case.

Iris

[shi]

Hi Iris,

I'm a little bit confused here. Are you talking about the reads mapped to different chromosomal locations, or you are talking about the reads that have only one reported mapping location but were found to overlap with more than one gene?

If it is the latter, featureCounts does not count those reads overlapping with more than one gene in its default setting and that shouldn't be the reason why featureCounts counted more reads than HTseq-count. A possible reason for this is that HTseq-count takes the rightmost chromosomal location of the feature as a open position (ie. not included in the summarization) whereas featureCounts takes it as a closed position (ie. included in the summarization).

We have seen what you reported in our evaluation as well (http://arxiv.org/abs/1305.3347). In our evaluation, ~500 genes got more counts from featureCounts and ~150 genes got more counts from HTseq-count when summarizing single-end reads. Did you see HTseq-count give you more counts for some genes in your data?

Thanks for suggesting us to add a parameter to filter for low mapping quality reads. We are now adding this parameter and will let you know once this is done.

Best wishes,
Wei

[iris_aurelia]

Hi Wei,

I compared HTseq-count with featureCounts, using a 50 bp single-end dataset, and noticed that for many genes featureCounts reports a higher count than HTseq-count does. It looks like that featureCounts does count reads mapped to multiple locations twice (or more), while HTseq-count excludes these reads. In principle we cannot say to which gene it actually belongs, so I would rahter discard those reads instead of counting them for each gene.
Wouldn't it be an option to include a parameter that can for example filter for low mapping quality reads or only use reads that are mapped uniquely?

Iris

[shi]

Dear Iris and Nicolas,

Thanks for letting me know it works for you now.

It is a non-trivial task to add an option to Subread to make it support the processing of gzipped FASTQ file directly. We are currently working on some further developments to Subread which takes a higher priority, but you may subscribe to the Subread updates (http://subread.sourceforge.net) so that you can get notified when we add this option.

Best wishes,

Wei

[Nicolas Nalpas]

Hi Wei,

Thanks a lot for all your help, indeed it works perfectly now and so fast.
This is great.
Thanks again,
Nicolas

[iris_aurelia]

Hi Wei,

Thanks, this solution works! However, I was using the Human reference (version 66), which had this problem as well. Maybe newer versions of the GTF file do not have this problem..?

Another thing I bumped into was the lack of reading gzipped fastQ files. Is it possible to either read from standard in or to create an option to use gzipped files?

Iris

[shi]

Dear Nicolas and iris_aurelia,

Please replace -g gene-id with -g ' gene_id' in your command (note that there is a space right before gene_id in the new argument) and try it again.

I suspect the problem you encountered was because there is a space right before the gene_id attribute in the 9th column of your annotation file. I have seen this for the yeast GTF annotation file generated by Ensembl. This is not compliant with the GTF-format specification which requires all the columns to be tab-delimited.

However, the human GTF annotation generated by Ensembl is fine. So it seems this problem only occurred for some species. I'm going to contact Ensembl to let them know this issue.

Please let me know if the new argument does not work.

Best wishes,

Wei

[iris_aurelia]

I tried subread-align and featureCount as well and encountered the same 'problem'. I was hoping that it would give me a gene identifier in the output file instead of line numbers. I've also used the -t and the -g option.

[Nicolas Nalpas]

Dear Wei,
I was just wondering if you could give advise on a featureCounts job that I tried to perform, I have RNA-seq data aligned to bovine genome (UMD3.1.71 from ENSEMBL) and I used subread-1.3.3-p2 (which is the latest version I think).
The error I get is:

No meta-feature id is found on the 1-th line. If it is a GTF file, you may need to check the name of the gene_id field and specify a correct field name using a '-g' opt
ion.

My code is:

featureCounts -a /path/annotation.gtf -t exon -g gene_id -b -i /path/file.bam -o file.out -s -R -p -B

And if I understand correctly the outputs I obtained, the reads assignment was still performed against my feature (exon in this case, represented as LINE_*), however I thought the output should have been my counts per gene_id:

geneid  length  nreads
LINE_0000001    88      0
LINE_0000002    167     0
LINE_0000003    51      0
LINE_0000004    337     0
LINE_0000005    1399    0
LINE_0000006    107     0
LINE_0000007    110     0
LINE_0000008    110     0
LINE_0000009    133     0
LINE_0000010    169     0
LINE_0000011    204     0
LINE_0000012    211     0
LINE_0000013    229     0
LINE_0000014    45      0
LINE_0000015    98      0
LINE_0000016    183     0
LINE_0000017    139     0
LINE_0000018    1429    0
LINE_0000019    218     0
LINE_0000020    91      0
LINE_0000021    202     0
LINE_0000022    142     2
LINE_0000023    157     6

My gtf file looks like this:

4       protein_coding  exon    432160  432247  .       -       .        gene_id "ENSBTAG00000003625"; transcript_id "ENSBTAT00000004727"; exon_number "1"; gene_name "BT.66099"; gene_biotype "protein_coding"; transcript_name "BT.66099-201"; exon_id "ENSBTAE00000393952";
4       protein_coding  CDS     432160  432235  .       -       0        gene_id "ENSBTAG00000003625"; transcript_id "ENSBTAT00000004727"; exon_number "1"; gene_name "BT.66099"; gene_biotype "protein_coding"; transcript_name "BT.66099-201"; protein_id "ENSBTAP00000004727";
4       protein_coding  start_codon     432233  432235  .       -       0        gene_id "ENSBTAG00000003625"; transcript_id "ENSBTAT00000004727"; exon_number "1"; gene_name "BT.66099"; gene_biotype "protein_coding"; transcript_name "BT.66099-201";
4       protein_coding  exon    430223  430389  .       -       .        gene_id "ENSBTAG00000003625"; transcript_id "ENSBTAT00000004727"; exon_number "2"; gene_name "BT.66099"; gene_biotype "protein_coding"; transcript_name "BT.66099-201"; exon_id "ENSBTAE00000037981";
4       protein_coding  CDS     430223  430389  .       -       2        gene_id "ENSBTAG00000003625"; transcript_id "ENSBTAT00000004727"; exon_number "2"; gene_name "BT.66099"; gene_biotype "protein_coding"; transcript_name "BT.66099-201"; protein_id "ENSBTAP00000004727";
4       protein_coding  exon    429103  429153  .       -       .        gene_id "ENSBTAG00000003625"; transcript_id "ENSBTAT00000004727"; exon_number "3"; gene_name "BT.66099"; gene_biotype "protein_coding"; transcript_name "BT.66099-201"; exon_id "ENSBTAE00000369947";
4       protein_coding  CDS     429103  429153  .       -       0        gene_id "ENSBTAG00000003625"; transcript_id "ENSBTAT00000004727"; exon_number "3"; gene_name "BT.66099"; gene_biotype "protein_coding"; transcript_name "BT.66099-201"; protein_id "ENSBTAP00000004727";
4       protein_coding  exon    426030  426366  .       -       .        gene_id "ENSBTAG00000003625"; transcript_id "ENSBTAT00000004727"; exon_number "4"; gene_name "BT.66099"; gene_biotype "protein_coding"; transcript_name "BT.66099-201"; exon_id "ENSBTAE00000037983";

Any idea on how I can sort this out?
Thanks a lot for your help.
Regards,
Nicolas

[shi]

Thanks for letting me know, NateP.

Cheers,
Wei