[shi]

Dear bruce01,

The featureCounts results just tell you that after trimming you got a lot more fragments mapped to regions where there is more than one gene residing. In fact, for the untrimmed reads you got about 23 million fragments successfully assigned (the groups 'ACCEPTED_GENE' and 'ACCEPTED_2VOTE_GENE'), whereas for the trimmed reads you got only about 15 million fragments assigned. So you lost 8 million assigned fragments. I bet these 8 million reads were almost all affected by your trimming and their mapping locations were changed. But the strange thing here is why they were all remapped to the regions where there are overlapping genes? You will have to ask the developers of the aligner you used about this.

Using Q20 to filter out read bases seems to result in the removal of a lot of good quality bases. This will make the mapping of RNA-seq data harder because reads become shorter, especially for those exon-spanning reads, from my point of view. We do not perform read trimming when mapping RNA-seq reads. Part of the reason for this is that we use the Subread aligner which can tolerate quite a few sequencing errors.


Hope this helps.

Best wishes,

Wei

[bruce01]

Sorry that I misunderstood the issue.

Still, I will try the Q filtering of featureCounts as it will remove a preprocessing step of the data. Thanks for your input.

[yangliao]

Sorry that I misunderstood the issue.

[bruce01]

Hi Wei,

I am wondering about some results I have had with tests I am running. I have been looking at appropriate quality cutoff values using Trimmomatic to trim reads based on quality score. I have tested multiple Q scores. Run parameters in featureCounts:

featureCounts -a /home/bin/my.gtf \
-i /home/data/quality/Q20/Aligned.out.sam \ 
-o /home/data/quality/Q20/Q20.featco.trimmo.counts \
-T 12 -p -P -C -R -D 500000

I run:

cut -f 2 *.counts.reads | sort | uniq -c > *.diags

giving me:

::::::::::::::
Q20.diags
::::::::::::::
  87805 ACCEPTED_2VOTE_GENE
14507634 ACCEPTED_GENE
4345122 MULTI_MAPPING
2315513 NOTFOUND_GENE
8709356 OVERLAPPED_GENES
  20001 PAIR_DISTANCE
::::::::::::::
Q40.diags
::::::::::::::
 149517 ACCEPTED_2VOTE_GENE
22640462 ACCEPTED_GENE
3425942 MULTI_MAPPING
4112768 NOTFOUND_GENE
  96425 OVERLAPPED_GENES
   2464 PAIR_DISTANCE

Can you see any reason that the untrimmed (ie Q40) alignment would have such a vastly reduced number of OVERLAPPED_GENES? Both input SAMs made in the same way using STAR with same parameters.

Anyone else with any ideas also very welcome to contribute!

[shi]

Dear Nicolas,

Thanks for letting me know it worked for you. That's great. We really appreciate you putting up with our bugs and helping us to improve the program. I hope you find it useful for your analyses in the future.

Best wishes,

Wei

[Nicolas Nalpas]

Dear Wei,
I have just tried the new version of featureCounts, it works perfectly (see below).

## UNSTRANDED ##	accepted_gene	multi_mapping	Notfound_gene	Overlapped_genes
featureCounts with SAM containing only unique	 3,181,215 	 -   	 794,232 	 11,469 
Htseq-count with SAM containing only unique	 3,174,243 	 -   	 810,359 	 2,314
				
## STRANDED ##	accepted_gene	multi_mapping	Notfound_gene	Overlapped_genes
featureCounts with SAM containing only unique	 2,784,181 	 -   	 1,198,726 	 4,009 
Htseq-count with SAM containing only unique	 2,783,943 	 -   	 1,202,275 	 698
				
## STRANDED REVERSE ##	accepted_gene	multi_mapping	Notfound_gene	Overlapped_genes
featureCounts with SAM containing only unique	 421,127 	 -   	 3,564,919 	 870 
Htseq-count with SAM containing only unique	 420,840 	 -   	 3,566,004 	 72

I am definitely switching to featureCounts for my new analyses, since it really is an amazing software.
Thanks again for all your help and for always being available.
Best wishes,
Nicolas

[shi]

Dear Nicolas,

We have just fixed this issue. Please try on the latest version 1.3.5-p3 (http://subread.sourceforge.net).

We have tested it and it worked fine in our end. Let me know if you ran into other problems.


Best wishes,

Wei

[sdriscoll]

I've written several utilities for post-processing on sam/bam files and I've learned that aligners definitely don't all report things in the same way. Another one is how they report singleton mates. Seems like there are a couple conventions for this.

[shi]

Dear Nicolas,

Many thanks for sending us your read data which helped us to identify the problem. The problem was that in your SAM file about half the fragments had their forward read appearing before their reverse read and the other half had their forward read appearing after their reverse read. This is not what featureCounts expects to see since most aligners always put the forward read before the reverse read in their mapping output.

The change of read order does not affect unstranded read counting with featureCounts, but it resulted in around half of the fragments not being counted when stranded read counting was performed, which is what you observed.

Nonetheless, we are now modifying featureCounts to deal with this situation and a patch will be released soon.

Best wishes,

Wei

[Nicolas Nalpas]

Dear Wei,

Thanks for your reply.
Enclosed is a file containing my first 5,000 mapped reads (forgot to say in my earlier post that I used STAR to align my reads against Bos taurus 3.1.71 genome).
I checked again featureCounts on this file and I was able to reproduce my results of yesterday (using Bos taurus 3.1.71 annotation).
Also here is my commands for unstranded, stranded and reverse stranded:

featureCounts -a /workspace/storage/genomes/bostaurus/UMD3.1.71/annotation_file/Bos_taurus.UMD3.1.71.gtf -t exon -g gene_id -i 5000_reads.sam -o test.unstranded -s 0 -R -p -B
featureCounts -a /workspace/storage/genomes/bostaurus/UMD3.1.71/annotation_file/Bos_taurus.UMD3.1.71.gtf -t exon -g gene_id -i 5000_reads.sam -o test.stranded -s 1 -R -p -B
featureCounts -a /workspace/storage/genomes/bostaurus/UMD3.1.71/annotation_file/Bos_taurus.UMD3.1.71.gtf -t exon -g gene_id -i 5000_reads.sam -o test.reverse -s 2 -R -p -B

I hope this help, thanks a lot for all your help.
Regards,
Nicolas

[shi]

Dear Nicolas,

Thanks for reporting the problem in much details. I certainly agree with your explanation of the stranded and unstranded read summarization and this is actually the counting scheme we implemented in featureCounts program. So I'm a bit confused with what you have observed.

Is it possible that you can provide us the first 10 reads you used for summarization and also your command for running featureCounts? This will help us to figure out what went wrong. Thanks.


Best wishes,

Wei

[Nicolas Nalpas]

Dear Wei,


We have now tried featureCounts v1.3.5-p1 to perform unstranded, stranded and reverse-stranded read counting on our new strand-specific dataset. While the unstranded parameter for featureCounts gives results very similar to HTseq-count, it seems that when it comes to stranded or reverse-stranded parameters the results between featureCounts and HTseq-count are quite divergent (see pdf file enclosed). I am actually expecting results similar to the one provided by HTseq-count.

I may be wrong but I believe that at the moment, featureCounts is assigning reads only to genes present on the + strand when using the stranded (s 1) parameter. And when using the reverse stranded (s 2) parameter, reads will be assigned only to genes present on the – strand.

However, it does not correspond to what I would call strand-specific read counts (I cannot talk for other RNA-seq users, but other RNA-seq users in our lab agree with me). Indeed, s 1 and s 2 parameters should in fact be one single parameter since users want read counts summarisation for genes present on + and – strands. Strand-specific is not the difference between genes present on the + or – strand.

What strand-specific actually means is that:
- (1) A read mapping to the + strand (for example read with flag 0 or 99 [its mate having flag of 147]) can be assigned to a gene present on the + strand (in this case it is a sense read)
- (2) A read mapping to the – strand (for example read with flag 16 or 83 [its mate having flag of 163]) can be assigned to a gene present on the - strand (in this case it is also a sense read)
- (3) A read mapping to the + strand (for example read with flag 0 or 99 [its mate having flag of 147]) can be assigned to a gene present on the - strand (in this case it is an antisense read)
- (4) A read mapping to the – strand (for example read with flag 16 or 83 [its mate having flag of 163]) can be assigned to a gene present on the + strand (in this case it is also an antisense read)

So I believe that, firstly when using the unstranded parameter, gene counts summarisation should contain summarisation for all reads corresponding to these 4 possibilities listed above; secondly when using the stranded (s 1) parameter, gene counts summarisation should contain summarisation for all reads corresponding to possibility (1) and (2); and thirdly when using the reverse-stranded (s 2) parameter, gene counts summarisation should contain summarisation for all reads corresponding to possibility (3) and (4).

So it also means that if a read is overlapping a gene present on the + strand and a gene present on the – strand, this read should be considered as overlapped_genes (and not counted unless specified by –O parameter on command line). In other terms, when you add number of reads with accepted-genes obtained from s 1 and s 2, you should obtain very close to the number obtained from s 0, indeed a read counted in s 1 should not also be counted in s 2.

Thanks a lot for all your help with this, I hope it was understandable (not the easiest to explain via email, otherwise let me know). And if I may suggest, it will be great if you can check these unstranded, stranded and reverse-stranded parameters for read summarisation.
Thanks a lot,
Nicolas

[iris_aurelia]

Hi Wei,

Thanks, I'm definitely gonna try this out, since featureCounts is much faster than HTseq-count.
I'll let you know if I run into any other problems!

Iris

[shi]

Dear Iris,

Thanks for looking into this deeply and providing a reproducible example.

You were correct that featureCounts did not exclude multi-mapping reads from the summarization. We have just added a '-M' option to featureCounts to enable it to disregard multi-mapping reads during read summarization. When '-M' is provided, featureCounts will use the 'NH' tag to identify multi-mapping reads and then exclude them.

These changes have been incorporated into the latest version (1.3.4) of the Subread package. You can download it from http://subread.sourceforge.net . We have tested it and it worked fine.

Please let me know if you run into any problems.

Best wishes,

Wei

Dear Iris,

Sorry, we have just flipped the meaning of the '-M' option. So now when running featureCounts in its default setting, multi-mapping reads will be excluded. The '-M' option is only useful if you want to count those multi-mapping reads (this is certainly not your case here).

The latest version is still 1.3.4 (the usage info for featureCounts and users guide has been updated).

Best wishes,

Wei

[shi]

Dear Iris,

Thanks for looking into this deeply and providing a reproducible example.

You were correct that featureCounts did not exclude multi-mapping reads from the summarization. We have just added a '-M' option to featureCounts to enable it to disregard multi-mapping reads during read summarization. When '-M' is provided, featureCounts will use the 'NH' tag to identify multi-mapping reads and then exclude them.

These changes have been incorporated into the latest version (1.3.4) of the Subread package. You can download it from http://subread.sourceforge.net . We have tested it and it worked fine.

Please let me know if you run into any problems.

Best wishes,

Wei