[also.casey]

Can I modify the linked gtf in any way to allow for summarization of the mapped reads as aligned?

[roryk]

The ensemble file you posted are chromosome coordinates but you aligned to a FASTA file of transcript sequences, not the whole chromosomes.

[also.casey]

Thanks for your reply roryk. I checked that reads were aligned using:

samtools view foo.bam | awk '$3 != "*"' | head -1

which returns:

"HWI-ST913:209:C2V7JACXX:6:1101:1587:2173 83 ENST00000394020 4460 0 58M = 4417 -101 CCTTATAGTTCTCTGCACACAATAGGTGCTTCTTGGATGTTCTTGGGTTTGGAAATAA IGHGGCCFHHFCHIIHIGJIJIJIHHAEB3EFD<CGGIJIIIIIHGIGGIHHIHIGJH NM:i:0 MD:Z:58 AS:i:58 XS:i:58"

Counting reads mapped to reference sequence names such as the one above in the example bam files using:

samtools view foo.bam | awk '$3 != "*"' | wc -l

indicates ~58% of reads were mapped in total, with ~3.5% mapping to ERCC reference sequences. Reads were aligned, but aren't being counted by featureCounts it seems.

[roryk]

I'm attempting to evaluate featureCounts to summarize a collection of bam files from an RNA-Seq experiment aligned to ensembl sequences. FASTA and GTF files for ensemble transcripts were obtained from the ensembl FTP site, modified by appending ERCC sequences and descriptors, respectively, and aligned using bwa-mem.

featureCounts happily processes the resulting bam files when using the ERCC GTF info only (and count totals seem to match that obtained using summarizeOverlaps in R), however using the ensembl GTF alone results in 0 "successfully assigned reads", while the combined ERCC and ensembl transcript GTF file results in only ERCC mapping reads being quantified (i.e. all ensembl transcripts appear in the output file, but all 0 counts). I think the issue may be because of some sort of formatting issue in the ensembl GTF? Can anyone shed any light on this issue?

The head of ensembl gtf I'm using is attached.

Thanks!

Is your combined FASTA file or the index you made for bwa-mem messed up? You would get this behavior if the reads were not aligned to the genome at all, only the ERCC sequences.

[also.casey]

Issues using ensembl (version 75) GTF

I'm attempting to evaluate featureCounts to summarize a collection of bam files from an RNA-Seq experiment aligned to ensembl sequences. FASTA and GTF files for ensemble transcripts were obtained from the ensembl FTP site, modified by appending ERCC sequences and descriptors, respectively, and aligned using bwa-mem.

featureCounts happily processes the resulting bam files when using the ERCC GTF info only (and count totals seem to match that obtained using summarizeOverlaps in R), however using the ensembl GTF alone results in 0 "successfully assigned reads", while the combined ERCC and ensembl transcript GTF file results in only ERCC mapping reads being quantified (i.e. all ensembl transcripts appear in the output file, but all 0 counts). I think the issue may be because of some sort of formatting issue in the ensembl GTF? Can anyone shed any light on this issue?

The head of ensembl gtf I'm using is attached.

Thanks!

[shi]

We have just released a new version (1.4.4). Please try that one to see if the problem persists.

Cheers,
Wei

[roryk]

Hi Wei,

featureCounts is awesome-- thanks for writing such nice software. We use it in the bcbio-nextgen pipeline for generating counts and it is insanely fast.

We have encountered this error on some alignment files when running featureCounts, version 1.4.3-p:

*** glibc detected *** /n/hsphS10/hsphfs1/chb/local/bin/featureCounts: double free or corruption (!prev): 0x0d7312f0 ***
======= Backtrace: =========
[0x80cfb0c]
[0x80d2e3b]
[0x808e405]
======= Memory map: ========
08048000-08133000 r-xp 00000000 00:20 366245842982 /net/hsphfs1/srv/export/hsphfs1/share_root/chb/local/bin/featureCounts
08133000-08134000 rw-p 000eb000 00:20 366245842982 /net/hsphfs1/srv/export/hsphfs1/share_root/chb/local/bin/featureCounts
08134000-08138000 rw-p 00000000 00:00 0
084e6000-0d744000 rw-p 00000000 00:00 0 [heap]
55555000-55556000 r-xp 00000000 00:00 0 [vdso]
55556000-55855000 rw-p 00000000 00:00 0
55857000-5644c000 rw-p 00000000 00:00 0
5644c000-57333000 rw-p 00000000 00:00 0
57333000-58391000 rw-p 00000000 00:00 0
58400000-58421000 rw-p 00000000 00:00 0
58421000-58500000 ---p 00000000 00:00 0
589f7000-589f8000 ---p 00000000 00:00 0
589f8000-59f83000 rw-p 00000000 00:00 0
fff96000-fffac000 rw-p 00000000 00:00 0 [stack]

The command is:

/n/hsphS10/hsphfs1/chb/local/bin/featureCounts -a /net/hsphfs1/srv/export/hsphfs1/share_root/chb/biodata/genomes/Mmusculus/mm10/rnaseq/ref-transcripts.gtf -o /net/hsphfs1/srv/export/hsphfs1/share_root/chb/projects/ljg_cardiac_rnaseq/project/work/htseq-count/tx/tmpJu9kkp/CRUSH2673-active.counts -s 0 -p -B -C /net/hsphfs1/srv/export/hsphfs1/share_root/chb/projects/ljg_cardiac_rnaseq/project/work/align/CRUSH2673-active/5_2014-03-12_project_star/CRUSH2673-active.sorted.bam

The header looks fine and the BAM and GTF files are also all good. Any thoughts?

[sindrle]

Its only because me, and maybe more, do not have a background in programming at all.

Its very simple to follow an explained guide, but with FeatureCounts I have to read and search for every step I take. Taking time which may be saved using HTSeq.

[shi]

Thank you, you have been very helpful.

Im currently struggling with how to implement featureCounts to DEXseq, its much easier to use HTSeq, since its described in the manual.

I really want to be able to use featureCounts, since with HTSeq I had to "hack" the python script to work with my annotation file, whereas I don't need that with featureCounts.

Also, I want to be able, in the end, to present expression values of transcripts as counts/library size/gene length, which also is easier with featureCounts, since it outputs gene length!

Could you elaborate a bit the difficulty you encountered with using featureCounts output with dexseq? Basically, featureCounts just outputs a read count table and it can be readily used by downstream DE analysis programs such as limma/voom and edgeR. The analysis can be performed at both gene level and exon level with these programs.

[shi]

hi yangliao,
Thanks for the reply. Ok, so I understand that in non-PE mode, sorting won't be performed. But in that case neither -B and -C flags would have any meaning.
Both the flags are really good for stringency/ QC. Moreover in non-PE mode, reads where R1 is in one meta-feature and R2 in another, would both get counted. That is to say the multi-overlap filter would become less effective.

I guess I would bear with the extra time. Thank you for a really good tool!

I agree it is better to count read pairs (PE mode) rather than individual reads for paired end read data.

[shi]

hello featureCounts users.
I have a bunch of "name-sorted" SAM files which I used for HTSeq. I'm trying featureCounts now but I cant find an option to turn off sorting.

Though much faster than HTSeq, the program spends major portion of its run time in sorting. If I could turn off sorting I could get results like superfast..

Is there a way?

thanks

featureCounts only performs read reordering when it is necessary, ie. when reads are not properly paired and/or there are missing mates. Note that sorting a SAM file by name using a tool like samtools does not guarantee that reads are properly paired after sorting. For example, reads that were reported multiple times in the SAM file may still not be properly paired after being sorted by read names.

If all the reads in a SAM/BAM file are properly paired, featureCounts will not perform read reordering at all although it will spend a very small amount of time to check if reads are properly paired or not. It is very important to check this, otherwise the read summarization results might be incorrect and the users may not be aware of it.

[yangliao]

Both reads in a pair are counted on the Single-end mode even if they are assigned to different meta-features.

Thank you for using our tools

hi yangliao,
Thanks for the reply. Ok, so I understand that in non-PE mode, sorting won't be performed. But in that case neither -B and -C flags would have any meaning.
Both the flags are really good for stringency/ QC. Moreover in non-PE mode, reads where R1 is in one meta-feature and R2 in another, would both get counted. That is to say the multi-overlap filter would become less effective.

I guess I would bear with the extra time. Thank you for a really good tool!

[amitm]

hi yangliao,
Thanks for the reply. Ok, so I understand that in non-PE mode, sorting won't be performed. But in that case neither -B and -C flags would have any meaning.
Both the flags are really good for stringency/ QC. Moreover in non-PE mode, reads where R1 is in one meta-feature and R2 in another, would both get counted. That is to say the multi-overlap filter would become less effective.

I guess I would bear with the extra time. Thank you for a really good tool!

[yangliao]

Dear Amitm,

Thank you for using featureCounts. If you wish to count the reads on the paired-end mode, and if the the input SAM or BAM files are not sorted by the read names, featureCounts has to sort the reads according to their names, so that it can count both reads in a pair at the same time. If, however, the reads are already sorted by names (like SAM or BAM files from subread, bowtie, bwa and etc), subread can count the reads without re-sorting them.

You can bypass the sorting step by turning off paired-end counting, or using subtools (provided in the /bin/utilities directory in subread-1.4.3 and latter versions) to sort the SAM file by read names before running featureCounts.

Cheers,

Yang

hello featureCounts users.
I have a bunch of "name-sorted" SAM files which I used for HTSeq. I'm trying featureCounts now but I cant find an option to turn off sorting.

Though much faster than HTSeq, the program spends major portion of its run time in sorting. If I could turn off sorting I could get results like superfast..

Is there a way?

thanks

[amitm]

Using pre-sorted SAM in featureCounts?

hello featureCounts users.
I have a bunch of "name-sorted" SAM files which I used for HTSeq. I'm trying featureCounts now but I cant find an option to turn off sorting.

Though much faster than HTSeq, the program spends major portion of its run time in sorting. If I could turn off sorting I could get results like superfast..

Is there a way?

thanks