Showing posts with label RPKM. Show all posts
Showing posts with label RPKM. Show all posts

Thursday, November 29, 2012

measuring expression of smallRNA genes

As the early note mentioned, it's not a good idea to use bias-correction based method like cufflinks to calculate the RPKM/FPKM values for small RNA (simply because it's not designed for that purpose).

Since the smallRNA is generally shorter, the usual RNA-seq protocol does not capture them properly (only the 'tail' of the distribution of fragment length). A good solution is to use smallRNA RNA-seq protocol. I meant using short reads length and without size selection.

Based on that, we may be able to use tools like cufflinks/DESeq to calculate the RPKM/FPKM. But again, I am not sure cufflinks will still apply the bias-correction to the small RNA genes when the reads length is short enough. As this post discussed, "The Best way to normalize miRNA HTS data", methods can be:

- RPM ( Read Per Million )
- RPKM (Reads per kilobase per million mapped) : Mortazi et al, Nat. Methods, 2008 (http://www.ncbi.nlm.nih.gov/pubmed/18516045)
- Trimmed mean of M-values : Robinson, Oshlack, Genome Biology 2010 (http://genomebiology.com/2010/11/3/R25)
- upper-quantile : http://www.ncbi.nlm.nih.gov/pubmed/20167110

Simon (one of the authors of DESeq) suggested to use TMM or DESeq (which are reads count-based methods) for differential expression analysis, at least for small RNA genes.

As the DEseq document described, to get reads count you need alignment (SAM/BAM) and the annotation file as input for some tools, like
  • the function summarizeOverlaps in the BioConductor GenomicRanges package,
  • the geneCounts function in the Bioconductor easyRNASeq package,
  • the htseq-count script distributed with the HTSeq Python package (Note:  htseq-count produces one count file for each sample. DESeq offers the function newCountDataSetFromHTSeqCount to get an analysis started from these files)  
  • or the bedtools multicov -count -bams sample1.bam sample2.bam sample3.bam -bed genes.bed
Note the naming issue in Refseq or UCSC that same locus has been assigned with multiple gene names (e.g. SMN1/2), and htseq-count requires unique gene_id in the annotation. You might want to make sure of this before using htseq-count. One solution to this is to use Ensembl GTF.

Also, be aware of the difference between bedtools multicov and htseq-count, esp. when you use multi-mappers. Read details here: http://seqanswers.com/forums/showthread.php?t=24903