Calculating small RNA expressions from Solexa Sequencing with Cufflinks

[DrD2009]

I am in need of your collective bioinformatic brains.

Here is my situation:
I have data sets of small RNAs from Solexa sequencing. I'm wanting to calculate the expressions of these small RNAs, but having some issues. To calculate expressions I used SAM output from Bowtie and used Cufflinks to calculate the expressions. The expressions calculated people in my lab are having concerns about. The small RNAs should be <35bps, but the calculations from Cufflinks involve regions that are much larger than 35bps. People in my lab wish to calculate expression for only the regions <35bps and I can only think of doing this by doing counts. Is Cufflinks an inappropriate expression calculation tool for small RNAs (miRNAs, siRNAs, etc.)? Or is there another program I could used to achieve this? I think Cufflinks is calculating the expressions correctly, but I thought I would reach out to all of you and ask.

Thanks in advance.


-Brandon

[minghui]

first,i do not understand what your mean to the ground. maybe you mean that you only want to calculate the read whose length > 35 bp ? if so ,you can write a script to filter your dataset before calculation.

[DrD2009]

I'll try to clear that up. I guess I said it in a confusing way.

We calculated the expressions in Cufflinks for the small RNAs that were mapped, but the region size for the expression calculations range from ~24 to ~7,000 bps. I think the following is occuring when Cufflinks calculates the expressions:

(small RNAs aligned to genome)
_______ _______
___ _______ ________
________ _______

[------------------------] <--- Size of region grouped by Cufflinks and expression value calculated for.

What the members of my lab want is the following:

(small RNAs aligned to genome)
_______ _______
___ _______ ________
________ _______

[-------] [-------]
[--] [------] [--------]
[-------] [-------] <----Size of regions calculated for expression.


Basically, they want a way to measure the expression of each read so they remain the size of small RNAs and are not lengthened, which seems to happen in Cufflinks, but Cufflinks is built for genes as opposed to small RNA analysis. The only way I can currently think of how to do expression calculations of each read would be to simply perform unique read counts.

[minghui]

ok,i understand what you mean exactly . Cufflinks measures transcript abundances in Fragments Per Kilobase of exon per Million fragments mapped (FPKM). it is fit for mRNA-seq .et , but, not smallRNA. Because of the length of smallRNA, one read one molecule. when we use FPKM , the same expression may have different FRKM, because the different length of different smallRNA. So, the exact regions you want may be useless. In my opinion:

two solutions :
1.use rpm (reads per million reads)
(OR)2.use software align the reads to genome , define each region yourself, count the reads in each region.

""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
(small RNAs aligned to genome)
_______ _______
___ _______ ________
________ _______

[-------] [-------]
[--] [------] [--------]
[-------] [-------] <----Size of regions calculated for expression.
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

[DrD2009]

Minghui,

Thanks for getting back with me on this.

I think we are going to do counts. Do you know if there is a preferred way to calculate expressions for small RNAs? Counts versus RPM? And do you know a software that calculates RPM?

I'm using Bowtie as the aligner.


Thanks again,
Brandon

[minghui]

hi,Brandon!

I am sorry ,i do not know any software that calculates RPM. Because I calculated RPM myself. There is no difference between counts and RPM,when you calculate in a single library. counts/million=RPM. RPM, is a normalization between two or many libraries. I think that you could try to align all reads to genome ,and calculate the numbers of reads in each region.

(small RNAs aligned to genome)
_______ _______
___ _______ ________
________ _______

[-------] [-------]
[--] [------] [--------]
[-------] [-------] <----Size of regions calculated for expression.

For example: ATGCATGCATGC ATGGATGCATGC TGCACGATCGAT (3 reads)

alignment :----------------------------------1----------2--3----------4
-----------(genome sequence) GGGGGGTAGCGATGCATGCATGCACGATCGAT
-----------(read)--------------------------- ATGCATGCATGC
-----------(read)--------------------------- ATGGATGCATGC
-----------(read)---------------------------------------TGCACGATCGAT

Calculation:"1-3" region : expression level :2
"2-4" region : expression level :1

Advantage: when couts unique sequences ,these thress have the same expression level 1; but that may be wrong,because RNA edit or sequencing errors.

[DrD2009]

Minghui,

Thanks for explaining that to me it makes sense. From the papers I've been reading lately they have all been using counts as well so I suppose that is how we will go about it too.

Thanks again.