[blancha]

featureCounts' lack of strict intersect mode for paired end reads

Why doesn't featureCounts have a strict intersect mode for paired end reads, like htseq-count?
It's a great software in many other aspects, but this is a mistake

I've put a very simple example in attachment where htseq-count will count one fragment, and featureCounts 2 fragments.

To me, the count of one for HSPB6 is correct. There is one fragment that aligns entirely in the gene, while only one read of the other fragment aligns to the gene, therefore it should not be counted.

I did set the minimum overlap to 50, even though it is inconvenient if I trim my reads. I tried with the minimum overlap to 100, but then I get a count of 0 for all the genes since it only applies to reads and not fragments.

[sdriscoll]

ah, with cufflinks did you use the -G or -g option? also how did you determine that cufflinks assigned more reads to features?

you could also give htseq-count a try just for a second opinion. that does basically what featureCounts does. another thing you can try is to align directly to transcripts with bowtie and see how many reads successfully align. that may be the best ground truth to how much of your data actually intersects the gene annotation.

[travelk]

Sorry, I worded it poorly. Indeed, I provided STAR with a fasta file to generate the genome files it aligns my reads to which would suggest it is genome based. To clarify, I provide the exact same gtf file for both FeatureCounts and Cufflinks. (I use the fasta and gtf files provided together by Ensembl.) When you recommend to run FeatureCounts with the gtf made by cufflinks, do you mean the the transcripts.gtf file it generates along with the genes.fpkm_tracking? I'm a bit confused about which file you are talking about...

[sdriscoll]

Hi Everyone,

I have a quick fundamental question about featureCounts.

I ran my reads through STAR and used featureCounts to get the counts of reads to that were successfully aligned.

Despite mapping my reads to the transcriptome in both STAR and featureCounts, I often get a lot of Unassigned_NoFeatures.

featureCounts -p -a $dirgtf -t exon -f -g gene_id -O -o ${sample}featurecounts.txt ${sample}star_sort.bam

Status	{sample}star_sort.bam
Assigned	3478951
Unassigned_Ambiguity	0
Unassigned_MultiMapping	992009
Unassigned_NoFeatures	3061428
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0
Unassigned_Nonjunction	0
Unassigned_Duplicate	0

In contrast, when I run the STAR bam files through cufflinks, it does pick up all the reads and so I am getting huge discrepancies between my featureCounts readout and my cufflinks readout.

I've also run the same fastq files through tophat and got a similar result.

Status	{sample}tophat_sort.bam
Assigned	3306323
Unassigned_Ambiguity	0
Unassigned_MultiMapping	1625816
Unassigned_NoFeatures	2938541
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0
Unassigned_Nonjunction	0
Unassigned_Duplicate	0

Could someone explain why despite being successfully aligned by either tophat or STAR, featureCounts returns almost 50% Unassigned_NoFeatures? What would be the reason a read is mapped by the aligner but not by featureCounts? I tried to look in the featureCounts user guide and original article but couldn't figure it out.

Thanks in advance.

Just to be sure, you say you aligned to the transcriptome with both STAR and featurecounts. Do you mean you aligned to a genome with STAR and then used featureCounts to count hits to transcripts with a GTF? That's how it should work and if that is what you did then the discrepancy is between then GTF you have and the potential extent of unannotated features in your actual samples. Cufflinks, remember, actually assembles a transcriptome and quantifies it in the same run so it is likely to use more of your aligned data in counts to genes. FeatureCounts is bound by your GTF. It seems likely that 50% of your data is mapping to in annotated features. That's all. Consider running featureCounts with the gtf made by cufflinks. You can annotate that gtf by using cuffcompare with the -r <gtf annotation> to pull in the information from the Gtf you are already using.

[Michael Love]

hi Shi and sindrle,

mcols(object)$basepairs when defined by the user should be the union of basepairs from all exons (so no double counting of basepairs which are in more than one exon), so the same as fc$annotation$Length.

[travelk]

Hi Everyone,

I have a quick fundamental question about featureCounts.

I ran my reads through STAR and used featureCounts to get the counts of reads to that were successfully aligned.

Despite mapping my reads to the transcriptome in both STAR and featureCounts, I often get a lot of Unassigned_NoFeatures.

featureCounts -p -a $dirgtf -t exon -f -g gene_id -O -o ${sample}featurecounts.txt ${sample}star_sort.bam

Status	{sample}star_sort.bam
Assigned	3478951
Unassigned_Ambiguity	0
Unassigned_MultiMapping	992009
Unassigned_NoFeatures	3061428
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0
Unassigned_Nonjunction	0
Unassigned_Duplicate	0

In contrast, when I run the STAR bam files through cufflinks, it does pick up all the reads and so I am getting huge discrepancies between my featureCounts readout and my cufflinks readout.

I've also run the same fastq files through tophat and got a similar result.

Status	{sample}tophat_sort.bam
Assigned	3306323
Unassigned_Ambiguity	0
Unassigned_MultiMapping	1625816
Unassigned_NoFeatures	2938541
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0
Unassigned_Nonjunction	0
Unassigned_Duplicate	0

Could someone explain why despite being successfully aligned by either tophat or STAR, featureCounts returns almost 50% Unassigned_NoFeatures? What would be the reason a read is mapped by the aligner but not by featureCounts? I tried to look in the featureCounts user guide and original article but couldn't figure it out.

Thanks in advance.

[shi]

I'm not sure what mcols(object)$basepairs represents, but fc$annotation$Length gives you the total number of exonic bases each gene has (each exonic base is counted exactly once). You probably will get more help if you post to bioconductor mailing list.

You may also have a look at a case study of using featureCounts in RNA-seq data analysis, which includes the step of generating RPKM values for genes if that's what you want to have:

http://bioinf.wehi.edu.au/RNAseqCaseStudy/

[sindrle]

Im actually struggling with the following, I want to easily create a DESeqDataSet from FeatureCounts. My goal is to use the
fpkm() function within DESeq2, not the one in FC.

In the DESeq2 help it says Ill get a matrix which is normalized per kilobase of the union of basepairs in the GRangesList or GRanges
of the mcols(object), and per million of mapped fragments, either using the robust median ratio
method (robust=TRUE, default) or using raw counts (robust=FALSE). Defining a column mcols(object)$basepairs takes precedence over internal calculation of the kilobases for each row.

My question is how I get the mcols(object)$basepairs from FC? Does it correspond to the fc$annotation$Length column?

So this is what I did, but is it correct?

mcols(ddsFullCountTable) <- fc$annotation$Length
mcols(ddsFullCountTable)$value
mcols(ddsFullCountTable) <- rename(mcols(ddsFullCountTable), c("value"="basepairs"))
mcols(ddsFullCountTable)$basepairs
fpkm <- fpkm(ddsFullCountTable)

[shi]

Since you specified '-O' option, which instructs featureCounts to assign the read to all its overlapping features (they are exons in this case since you specified '-f'), reads will be assigned to all their overlapping exons. That says a read will be counted multiple times if it overlaps multiple exons. You can see this in the assignment results for two of the exons:

ENSG00000223972 1 13221 14409 + 1189 4
ENSG00000223972 1 13225 14412 + 1188 4

If you do not specify the '-O' option, reads overlapping with multiple exons will not be counted. Thus, you possibly will get zero count for both the two exons above.

[ndaniel]

how exons/features are treated?

In case that one counts the reads mapping on exons, how does featureCounts treats the exons/feature which overlap (and they belong to the same gene)?

More exactly, if there are two exons which belong to the same gene such that:
(i) each exon belongs to different isoforms of the same gene, and
(ii) both exons have the same start position on the genome, and
(iii) second exon is longer than the first one
how the reads count for these exons is reported? Is reported the count of reads which (i) map on the intersection of these two exons, or (ii) map on the longest exon in this case (because the longest exons contains also the shortest one)? Are there reported two counts of reads for each exon?

Here is how the featureCounts is run for these questions (i.e. counting the reads mapping on exons)
# Program:featureCounts v1.4.4; Command:"featureCounts" "-T" "4" "-a" "/ensembl.gtf" "-t" "exon" "-f" "-g" "gene_id" "-O" "-o" "/output$
Geneid Chr Start End Strand Length /accepted_hits.sorted.bam
ENSG00000223972 1 11869 12227 + 359 0
ENSG00000223972 1 12613 12721 + 109 0
ENSG00000223972 1 13221 14409 + 1189 4
ENSG00000223972 1 11872 12227 + 356 0
ENSG00000223972 1 12613 12721 + 109 0
ENSG00000223972 1 13225 14412 + 1188 4
ENSG00000223972 1 11874 12227 + 354 0
ENSG00000223972 1 12595 12721 + 127 0
ENSG00000223972 1 13403 13655 + 253 3
ENSG00000223972 1 13661 14409 + 749 1
ENSG00000223972 1 12010 12057 + 48 0
ENSG00000223972 1 12179 12227 + 49 0
...

[shi]

Thanks!

Is this:


featureCounts \
-a organism.gtf \
-t exon \
-f \
-g exon_id \
-O \
-o counts_per_exons.txt \
accepted_hits.sorted.bam

or this (the difference is "-g gene_id" and "-g exon_id" ?


featureCounts \
-a organism.gtf \
-t exon \
-f \
-g gene_id \
-O \
-o counts_per_exons.txt \
accepted_hits.sorted.bam

Typically this should be "gene_id", but you may check your GTF annotation file to make sure it is correct (the 9th column).

Also, what are the command line parameters for running featureCounts for counting reads mapping on exons from command line in order to produce exactly the same output as this R version of featureCounts as shown below (taken from limma user guide which was updated last time on April 12, 2014, page 127):


fc_SE <- featureCounts(SE_bam_files, annot.ext="unique.gff",
+ isGTFAnnotationFile=TRUE,
+ GTF.featureType="exon", GTF.attrType="ID",
+ useMetaFeatures=FALSE, allowMultiOverlap=TRUE)

This should give you the same results as you got from sourceforge featureCounts. You need to adapt this R command to make it work for your data though. Also, please make sure you are using Rsubread 1.14.0 and sourceforget Subread 1.4.4 (these two are matched versions), otherwise you may get slightly different results.

[ndaniel]

featureCounts for counting reads mapping on exons

Replacing "-f exon" with just "-f" will instruct featureCounts to perform feature (eg exon) level summarization. The output will include number of reads assigned to each exon.

You may also want to turn on "-O" option, which will instruct featureCounts to assign a read to all its overlapping exons. This is often useful for exon-level read counting.

Wei

Thanks!

Is this:


featureCounts \
-a organism.gtf \
-t exon \
-f \
-g exon_id \
-O \
-o counts_per_exons.txt \
accepted_hits.sorted.bam

or this (the difference is "-g gene_id" and "-g exon_id" ?


featureCounts \
-a organism.gtf \
-t exon \
-f \
-g gene_id \
-O \
-o counts_per_exons.txt \
accepted_hits.sorted.bam

Also, what are the command line parameters for running featureCounts for counting reads mapping on exons from command line in order to produce exactly the same output as this R version of featureCounts as shown below (taken from limma user guide which was updated last time on April 12, 2014, page 127):


fc_SE <- featureCounts(SE_bam_files, annot.ext="unique.gff",
+ isGTFAnnotationFile=TRUE,
+ GTF.featureType="exon", GTF.attrType="ID",
+ useMetaFeatures=FALSE, allowMultiOverlap=TRUE)

[shi]

Replacing "-f exon" with just "-f" will instruct featureCounts to perform feature (eg exon) level summarization. The output will include number of reads assigned to each exon.

You may also want to turn on "-O" option, which will instruct featureCounts to assign a read to all its overlapping exons. This is often useful for exon-level read counting.

Wei

[ndaniel]

If this is used for counting reads mapping on genes using Ensembl GTF file:

featureCounts \
-a organism.gtf \
-t exon \
-g gene_id \
-o counts_per_gene.txt \
accepted_hits.sorted.bam

then how should be counted the reads mapping on exons using the same Ensembl GTF file?

Is it something like this?
featureCounts \
-a organism.gtf \
-t exon \
-f exon \
-g exon_id \
-o counts_per_exons.txt \
accepted_hits.sorted.bam

???

[roryk]

You could just count the third column in your SAM file. Aligning to the whole chromosomes is what you really want to do though.