[bruce01]

Hi Wei,

that makes sense for s=0, of course. I suppose my real question is that if we know what strand the RNA originally comes from do we not always need to use s=1+2, seeing as some will be from '+' and some '-'? Is this something other people do?

Bruce.

[shi]

Hi @bruce01,

The '-s 0' option is likely to give you less counts compared to the sum of counts from '-s 1' and '-s 2'. This is because reads are more likely to be found to overlap with more than one gene when '-s 0' is used (genes with overlapping coordinates but on different strands will be taken as overlapping genes at unstranded counting mode, but not taken as overlapping genes at stranded counting mode).

Wei

[bruce01]

Strand issues

Hi Wei,

I have been using single end 50bp stranded library, and was getting low counts and losing quite a lot of data. I tried '-s =0, 1, 2' and found discrepancies with counts. I expected 1+2 to sum to output of 0, is this correct?

ENS_ID          0    1    2
ENSG00000000003 48   15   33
ENSG00000000005 44   17   27
ENSG00000000419 77   50   91
ENSG00000000457 235  72   193
ENSG00000000460 677  581  360

I have 23,038 matching, and 40,642 which do not sum to output from '-s 0'.

I have tried resorting (aligned with STAR giving 'sort-by-coord' option) which made no difference to featureCounts output, and my command is:

featureCounts -a ${HOME}/Homo_sapiens.GRCh37.75.saf -F SAF -t gene -o $out -s "$x" -T 8 -d 30 $1

Any ideas would be very welcome and I appreciate your continued support of you tools.

[Gordon Smyth]

At the moment I do take the R route to generate RPM or RPKM, yet this is labour intensive for non-model organisms as I need to create multiple annotation objects. Seeing as all the information is available to featureCounts, I think there would be a demand for this.

Hi Colin,

We seem to be in complete agreement. Our group probably uses RPKM much the same as you do, for descriptive plots and expression summaries.

We don't output RPKM from featureCounts just on the principle that it seems redundant information. featureCounts already outputs all the information one needs to compute RPKM. You can use the rpkm() function in the edgeR package to extract the RPKM (as in the example I linked to) or else just do it yourself in a couple of lines:

fc <- featureCounts(stuff)
RPKM <- fc$counts / fc$annotation$Length * 1e3
RPKM <- t(t(RPKM) / colSums(fc$counts) * 1e6)

In our own analyses, we generally like to filter unwanted rows (mitochondrial genes, Ig DNA segments, ribosomal genes, non-expressed genes and what not) and estimate normalization factors before computing RPKM. So there are few steps to do before we extract RPKM.

BTW, did you know that voom and limma can work directly with RPKM? Just do this:

v <- voom(counts, design, plot=TRUE)
v$E <- v$E - log2(GeneLength/1e6)

Then v$E contains log2(RPKM), and all the limma analyses will work correctly in terms of these expression values.

[colindaven]

Counts are needed for a proper statistical analyses, as is explained in the publications behind edgeR, voom and DESeq.

However, generating RPKM from featureCounts is easy, see:

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

for an example.

Dear Gordon,

thanks for your response and the link, I agree entirely. Incidentally I have been using edgeR and more recently voom/limma for statistical analysis for many years. However, my users are still very keen to see RPKM values for many, many samples.

RPKM, despite it's statistical weaknesses, is used widely. At the moment I do take the R route to generate RPM or RPKM, yet this is labour intensive for non-model organisms as I need to create multiple annotation objects. Seeing as all the information is available to featureCounts, I think there would be a demand for this.

Like you say, the downside is that people then might want to put RPKM into statistical packages like edgeR, but this topic has been dealt with repeatedly and might be avoided by an appropriate warning.

Thanks again,
Colin

[tamfernandez]

In feautureCounts built-in annotation, multicopy genes (with the same Entrez gene id) were combined to form one gene that contains exons from all its copies. Assigning 1 count for a read overlapping with this meta-gene should be equivalent to assigning a 1/n count to each of its n copies. You may try the inbuilt annotation to see it works for you.

Wei

Thank you Wei for your ansewer. This will solve part of the problem. Maybe for genes like tRNA-Pro isodecoders, which have almost the same sequence but still different ids, I can try to modify the id's so I can use what you suggested.

Thank you,
Tamara

[shi]

In feautureCounts built-in annotation, multicopy genes (with the same Entrez gene id) were combined to form one gene that contains exons from all its copies. Assigning 1 count for a read overlapping with this meta-gene should be equivalent to assigning a 1/n count to each of its n copies. You may try the inbuilt annotation to see it works for you.

Wei

[tamfernandez]

Hi Tamara,

Although this is something we are considering to add to featureCounts, I do not really think you will get much gain from doing this since this will add noise to your read count data, which may then compromise your downstream differential expression analysis.

If you worry about those reads overlapping with multiple genes, you should consider using a more conservative annotation such as RefSeq gene annotation which contains much less overlapping genes compared to other annotations such as Ensembl genes.

Wei

Hi Wei,

Thank you for your reply! I absolutely agree with you but I'm trying to count short reads that map to some multicopy genes and some conserved genes like tRNAs among others. I thought that assigning a fraction to each gene copy could be a good solution to count reads that map specifically to these genes. Do you have any suggestions for dealing with this using featureCounts?

Again, thank you for your attention.
Tamara

[shi]

Dear Wei,

I would like to use featureCounts for counting multiple aligned reads, but instead of adding 1 to the count for each location y would like to add fractional counts (1/number of locations where the read mapped). Is this posible?

Thank you in advance for your reply.
Tamara

Hi Tamara,

Although this is something we are considering to add to featureCounts, I do not really think you will get much gain from doing this since this will add noise to your read count data, which may then compromise your downstream differential expression analysis.

If you worry about those reads overlapping with multiple genes, you should consider using a more conservative annotation such as RefSeq gene annotation which contains much less overlapping genes compared to other annotations such as Ensembl genes.

Wei

[Gordon Smyth]

what are your opinions on using RPKM or another alternative (unique counts, expected counts).

Can featureCounts generate RPKM, at least for comparison? Are there good reasons why it should not ?

Counts are needed for a proper statistical analyses, as is explained in the publications behind edgeR, voom and DESeq.

However, generating RPKM from featureCounts is easy, see:

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

for an example.

[tamfernandez]

Dear Wei,

I would like to use featureCounts for counting multiple aligned reads, but instead of adding 1 to the count for each location y would like to add fractional counts (1/number of locations where the read mapped). Is this posible?

Thank you in advance for your reply.
Tamara

[colindaven]

Hi everybody,

what are your opinions on using RPKM or another alternative (unique counts, expected counts).

Can featureCounts generate RPKM, at least for comparison? Are there good reasons why it should not ?

Thanks,
Colin

[emmanouela]

Hi Shi,

Its because the reads for the rest of the genes are reversely stranded (!!!). Actually there has been an issue with possible contamination of that sample, because this gene is not supposed to be expressed at all in this sample and yet there are quite a few reads there. Could those reads be coming from another sample and how would that be??? (possibly the wrong place to be asking that but anyway, if anyone has an answer it would be nice to know)

Thanks for your answer! It helped a lot!

[shi]

Why did you use the option '-s 2' (reversely stranded), given that your reads are located on the same strand as that of the target gene?

[emmanouela]

unassigned reads

Hi all,

I have a question regarding some reads that have been excluded as "Unassigned_NoFeatures" but I have no idea why, so any clues would be really helpful!

Here is my example. I have a gtf file and looking at this gene:

grep Hba-x genes_ucsc_mm9.gtf | grep exon
chr11	unknown	exon	32176600 32176758	.	+	.	gene_id "Hba-x"; transcript_id "NM_010405"; gene_name "Hba-x"; p_id "P12021"; tss_id "TSS5558";
chr11	unknown	exon	32177569 32177773	.	+	.	gene_id "Hba-x"; transcript_id "NM_010405"; gene_name "Hba-x"; p_id "P12021"; tss_id "TSS5558";
chr11	unknown	exon	32177880 32178115	.	+	.	gene_id "Hba-x"; transcript_id "NM_010405"; gene_name "Hba-x"; p_id "P12021"; tss_id "TSS5558";

I run featureCounts with this command:

featureCounts -T 5 -a genes_ucsc_mm9.gtf -t exon -g gene_id -R -o ../features_counted_rev_try -p  -s 2 my_bam.bam

Extracting some random reads that I know that are in the gene (from igv):

$samtools view my_bam.bam | grep 113:8147:70638
HISEQ2000-06:419:C4HB0ACXX:6:1113:8147:70638	73	chr11	32177640	50	100M	*	0	0	TGCGGGCCCACGGCTTCAAGATCATGACCGCCGTAGGGGATGCGGTTAAGAGCATCGACAACCTCTCTAGTGCTTTGACTAAGCTGAGCGAGCTGCATGC@@@AD@DDDGHHHIIGHGFHIIHE@FGIIHIIGGEAEEBBDDECB;?8:ACACCCABABDBBDACDDDDA4>C>@::@>>CDCDDCDDDDBBDBDDCDCD	AS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:100	YT:Z:UU	XS:A:+	NH:i:1
$samtools view my_bam.bam |  grep 1307:1932:2934
HISEQ2000-06:419:C4HB0ACXX:6:1307:1932:2934	99	chr11	32176639	50	100M	=	32176684	955	TCCAACCCTCACCACCACCACCACCATGTCTCTGATGAAGAATGAGAGAGCTATCATCATGTCCATGTGGGAGAAGATGGCTGCTCAGGCCGAGCCCATT	@<@FFDDFHFHFHIGFFEHIIGHGIIEGD>DDEHFG@?GD@G>DHIGDFHGGHHHIEFIGG=FEIICHG@=AAEFF@6@CCEBD3@@>??<=@BBD?B9(	AS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:100	YT:Z:UU	XS:A:+	NH:i:1
HISEQ2000-06:419:C4HB0ACXX:6:1307:1932:2934	147	chr11	32176684	50	75M810N25M	=	32176639	-955	GAGAGCTATCATCATGTCCATGTGGGAGAAGATGGCTGCTCAGGCCGAGCCCATTGGCACTGAGACTCTAGAGAGGCTCTTCTGCAGCTACCCCCAGACG	:>>>@::>CCAAAACCCCDCB?CAC@:CCDCCABB@>@A;;6=AAB=.9AIFBFHCB9*?89*?9?9??9JIGEIIHGGIIGIIJIIHFCHHDFFFFBCC	AS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:100	YT:Z:UU	XS:A:+	NH:i:1
$samtools view my_bam.bam |  grep 2305:5148:22998
HISEQ2000-06:419:C4HB0ACXX:6:2305:5148:22998	73	chr11	32176646	50	100M	*	0	0	CTCACCACCACCACCACCATGTCTCTGATGAAGAATGAGAGAGCTATCATCATGTCCATGTGGGAGAAGATGGCTGCTCAGGCCGAGCCCATTGGCACTGCCCFFDFFFHHHGJJJIGGIBHHIGJJJJJJJJIJIJJJJJIJJIIJJJJIIIIHIIJIGHIJJGIJJIGHHFHFBDFFEECEDDDDDDDDDCACDCDD9	AS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:100	YT:Z:UU	XS:A:+	NH:i:1

I see that according to their coordinates and the gtf they should be counted.

But in the report (-R) I get them as unassigned no features:

$ grep 113:8147:70638 my_bam.bam.featureCounts 
HISEQ2000-06:419:C4HB0ACXX:6:1113:8147:70638	Unassigned_NoFeatures	*	*
$ grep 1307:1932:2934 my_bam.bam.featureCounts 
HISEQ2000-06:419:C4HB0ACXX:6:1307:1932:2934	Unassigned_NoFeatures	*	*
$  grep 2305:5148:22998 my_bam.bam.featureCounts 
HISEQ2000-06:419:C4HB0ACXX:6:2305:5148:22998	Unassigned_NoFeatures	*	*

What am I missing???
Thank you for your help!