SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
RNA seq data normalization question slny Bioinformatics 35 10-19-2016 05:32 AM
RNA Seq normalization harshinamdar Bioinformatics 39 03-16-2013 01:12 AM
RNA-Seq: GC-Content Normalization for RNA-Seq Data. Newsbot! Literature Watch 0 12-20-2011 02:00 AM
RNA-seq library DSN normalization megalodon RNA Sequencing 2 09-01-2011 01:44 AM
clarification of rna-seq normalization frymor Bioinformatics 3 08-06-2011 06:07 PM

Reply
 
Thread Tools
Old 04-06-2010, 09:28 AM   #21
Siva
Member
 
Location: Ames, IA

Join Date: Nov 2008
Posts: 57
Default

Hi RCJ
Thanks. Btw, I did single end sequencing and not paired ends. I don't think that that should be a problem in calculating integer counts per transcript. I have contacted our statistics collaborator regarding FPKM and will get back after I get a reply.

best
Siva
Siva is offline   Reply With Quote
Old 04-06-2010, 10:01 AM   #22
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

In case this got lost in my lengthy post #12:

The reason why raw counts are advantageous to FPKM values for statistical inference is discussed in this thread, from post #6 onwards: http://seqanswers.com/forums/showthread.php?t=4349
Simon Anders is offline   Reply With Quote
Old 04-07-2010, 07:56 AM   #23
lpachter
Member
 
Location: Berkeley, cA

Join Date: Feb 2010
Posts: 40
Default

I'd like to point out again that in fact raw counts are not advantageous to FPKM.

FPKM is a unit for reporting an estimate of expression. Reports of expression estimates (whether in FPKM or other units) are by definition based on statistical inference from raw read counts.
lpachter is offline   Reply With Quote
Old 04-07-2010, 08:54 AM   #24
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Of course, FPKM is a proper unit to report expression, more so than raw counts.

However, if you want to perform a statistical test to decide the question whether expression in two conditions is significantly different, then a test that works on the level of raw counts is more powerful than one that works on the level of normalized expression, for the reasons that we discuss in our paper.

Of course, if I know transcript lengths and library sizes, I can always convert between one type of information and the other. If I have two test procedures, one that only asks for raw read counts, or one that only asks for FPKM values, with neither asking for transcript lengths, this two procedures will necessarily do something different.

Actually, not only DESeq but also cuffdiff works on the level of raw counts, as Cole commented here: http://seqanswers.com/forums/showpos...59&postcount=9

Finally, on the risk of being nit-picking: If one divides counts by transcript length or total reads, one not only changes the unit but also the type of quantity reported. If I tell you either that I traveled 10 miles or 13 km, this is the same quantity in different units. If it took me 2 hours and I only tell you that made 5 miles per hour, this is a different quantity, namely velocity rather than distance. Similarly, FPKM is a unit of expression, i.e., it is deemed to be proportional to molar concentration of transcript molecules. Counts is not a unit of expression. Nevertheless, it is advantageous to compare counts rather then expression values to do inference on differential expression.
Simon Anders is offline   Reply With Quote
Old 04-07-2010, 10:00 AM   #25
lpachter
Member
 
Location: Berkeley, cA

Join Date: Feb 2010
Posts: 40
Default

Perhaps I'll risk being nitpicky as well regarding a sentence in your nitpick.

You write that "Nevertheless, it is advantageous to compare counts rather then expression values to do inference on differential expression."

Now before you wrote _raw_ counts and now just counts. I agree with the latter (counts) but not with your statement before of using raw counts.

Consider the following simple example:
A gene with two isoforms A and B (with equal lengths), with the property that in experiment 1 isoform A is expressed at double the amount of B, vs. experiment 2 in which B is double A. The number of read counts from this gene may be identical in both experiments. However it may be possible to infer the changes in the expression levels of the individual transcripts. That is because of where the reads occur (e.g. one transcript may contain a splice junction and reads there may help to infer its abundance).

Now given _expression values_ rather than the raw counts, it may be clear that there has been differential expression. But looking at counts alone there would be no way to tell.

Having said this, I'm aware of your paper and I think your approach is valuable (and I agree with your argument that its better than the Poisson assumption). But what needs to be used is the probabilistic assignment of read counts to transcripts, rather than just raw read counts. And these _can_ be obtained from estimates of expression together with transcript lengths (as you point out). But again, they cannot be obtained from raw read counts alone.
lpachter is offline   Reply With Quote
Old 04-07-2010, 12:51 PM   #26
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Maybe our discussion revolved a bit on a misunderstanding regarding what we are aiming for. I want to get expression of genes, summing over transcripts, while you want to distinguish isoforms. (I may have worked too much with yeast and so got used to not having much alternative splicing going on in my data. )

If I don't care about isoforms or think that my coverage is too low to distinguish isoforms anyway, I expect to get optimal power by simply summing everything up.

Our (and edgeR's) negative binomial approach depends on the counts being "raw" (not normalized by anything) to get the shot noise right. So, if we want to use DESeq to distinguish isoforms, we would have to count reads for each exon and work on this level. (Actually, this might not work too well, as we would run into trouble with exons having different lengths in different isoforms.)

Cuffdiff is, as I understand it, designed to deal with such issues, while our approach ignores them. I expect that DESeq, in compensation for being unsuitable to detect differences in isoform proportions as in your example, achieves much better detection power for differences in total expression (per gene, summing over isoforms), especially at very low counts.

As I am not clear on how biological noise is taken into account by cuffdiff I cannot be fully sure whether this expectation will hold (and I'm quite curious to learn more about cuffdiff once your paper is out).
Simon Anders is offline   Reply With Quote
Old 04-07-2010, 09:08 PM   #27
Siva
Member
 
Location: Ames, IA

Join Date: Nov 2008
Posts: 57
Default

We are all learning a great deal from your discussions. Thanks, I now clearly understand many aspects of RNA-seq data normalization. Unfortunately, I am unable to convince our statistician to use cufflinks FPKM values even with the confidence intervals.

So I have two questions:
1. I would like to know if we can get transcript length in Cufflinks outputs so I can convert FPKM values to absolute read counts? I would assume this would need a script to automate.

2. Has someone got a script that uses the Bowtie SAM output to directly report counts per transcript and transcript length? Is there a software that I can use?

Siva
Siva is offline   Reply With Quote
Old 04-08-2010, 01:39 AM   #28
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Hi Siva

Quote:
Originally Posted by Siva View Post
2. Has someone got a script that uses the Bowtie SAM output to directly report counts per transcript and transcript length?
I do, and I've just written up the documentation for it yesterday: http://www-huber.embl.de/users/ander...doc/count.html

It is part of HTSeq, a Python framework to process HTS data, which I am working on currently. I'm still rounding some corners and will announce its release soon, but you can already use it, of course. Just tell me if you find any bugs or problems. For full documentation, see here: http://www-huber.embl.de/users/anders/HTSeq

Simon
Simon Anders is offline   Reply With Quote
Old 04-08-2010, 05:26 AM   #29
mgogol
Senior Member
 
Location: Kansas City

Join Date: Mar 2008
Posts: 197
Default one way to do it

Code:
bamToBed -i accepted_hits.sorted.bam > accepted_hits.sorted.bed
coverageBed -a accepted_hits.sorted.bed -b mm9.bed | sort -r -n -k7 > counts.txt
This will get the per gene counts from a sam file, but doesn't consider introns/exons.

bamToBed and coverageBed are from bedtools and mm9.bed is a bed file describing gene positions. I'm using this for genes though, not transcripts.

On second thought, maybe this is too crude.

Last edited by mgogol; 04-08-2010 at 06:01 AM.
mgogol is offline   Reply With Quote
Old 04-08-2010, 05:54 AM   #30
lpachter
Member
 
Location: Berkeley, cA

Join Date: Feb 2010
Posts: 40
Default

My personal advice to biologists performing RNA-Seq experiments is to care about isoforms. But that issue aside, as the paper by Li et al. from the Dewey lab recently pointed out, probabilistic assignment of reads is necessary even for genes, because many reads map to multiple genes. I therefore contest the statement that DESeq outperforms other methods simply because it uses raw read counts, even in yeast.
lpachter is offline   Reply With Quote
Old 04-08-2010, 10:33 AM   #31
steven
Senior Member
 
Location: Southern France

Join Date: Aug 2009
Posts: 269
Default

Quote:
Originally Posted by lpachter View Post
My personal advice to biologists performing RNA-Seq experiments is to care about isoforms.
I agree. Unless you believe that in your organism of interest the good old assumption "1 gene => 1 protein" is still valid..
steven is offline   Reply With Quote
Old 04-09-2010, 10:12 AM   #32
Siva
Member
 
Location: Ames, IA

Join Date: Nov 2008
Posts: 57
Default

Hi Simon
Thanks, I will get back to you after I use htseqcount.

Hi L Pachter
As a biologist, I take your advice seriously and understand your point about isoforms. I will do both Cuffdiff and DEseq for DE on my data. I will keep you posted.

thanks everyone
Siva

Last edited by Siva; 04-09-2010 at 10:15 AM. Reason: typo
Siva is offline   Reply With Quote
Old 04-28-2010, 11:02 AM   #33
Siva
Member
 
Location: Ames, IA

Join Date: Nov 2008
Posts: 57
Default Cufflinks...GTF...DESeq....analysis

I have done unpaired RNAseq on maize samples. I used cufflinks for transcript assembly on a series of 48 SAM files generated by Tophat. My aim is to do differential expression analysis over all these samples. (fyi I had used a Bowtie generated index of the genome for multiple alignment). I did not use any reference annotation. Cufflinks generates gene_iD or transcript ID like CUFF1.0 or CUFF1.1 etc. However I have observed that for each sample the same CUFF ID corresponds to different coordinates on the chromosome. I think Cuffdiff measures transcript abundance by tracking FPKM changes in transcripts sharing a common transcription start site and also by tracking different transcripts of a gene.

Programs like DESeq, EdgeR that look for differential expression needs read counts per gene_ID. I used htseqcount (check: http://www-huber.embl.de/users/ander...unt.html#count) using a SAM file and its corresponding Cufflinks generated GTF file. This generated read counts per cuff ID. If each CUFF ID corresponds to different portions on the chromosome how would DEseq or any other DE script sort this out?

any thoughts/suggestions appreciated.

Thanks
Siva

Last edited by Siva; 04-28-2010 at 01:17 PM. Reason: typing errors
Siva is offline   Reply With Quote
Old 04-28-2010, 06:54 PM   #34
Siva
Member
 
Location: Ames, IA

Join Date: Nov 2008
Posts: 57
Lightbulb

Quote:
Originally Posted by Siva View Post
I have done unpaired RNAseq on maize samples. I used cufflinks for transcript assembly on a series of 48 SAM files generated by Tophat. My aim is to do differential expression analysis over all these samples. (fyi I had used a Bowtie generated index of the genome for multiple alignment). I did not use any reference annotation. Cufflinks generates gene_iD or transcript ID like CUFF1.0 or CUFF1.1 etc. However I have observed that for each sample the same CUFF ID corresponds to different coordinates on the chromosome. I think Cuffdiff measures transcript abundance by tracking FPKM changes in transcripts sharing a common transcription start site and also by tracking different transcripts of a gene.

Programs like DESeq, EdgeR that look for differential expression needs read counts per gene_ID. I used htseqcount (check: http://www-huber.embl.de/users/ander...unt.html#count) using a SAM file and its corresponding Cufflinks generated GTF file. This generated read counts per cuff ID. If each CUFF ID corresponds to different portions on the chromosome how would DEseq or any other DE script sort this out?

any thoughts/suggestions appreciated.

Thanks
Siva
I have thought about this and I feel I am fundamentally wrong in using Cufflinks generated GTF output of one sample and the corresponding SAM file in htseqcount. Quoting from Cufflinks website; "Cufflinks constructs a parsimonious set of transcripts that explain the reads observed in an RNAseq expt." We know that it uses a pseudo ID to refer to various transcript bundles. If I use htseqcount using a SAM file against a cufflinks GTF I will only get read counts per pseudo ID. If I need read counts per gene or transcript i should have used a std reference annotation (GFF or GTF) to output counts per gene or transcript. However for many organisms including maize (see http://ftp.maizesequence.org/current/working-set/) the GFF3 format does not have the id attributes "gene_id" or "transcript id".

Am I right in assuming that I should not have used cufflinks generated GTF to output absolute counts per gene/transcript?

Can anyone help me with converting the std GFF to a GTF containing the necessary id attributes (viz., "gene_id" or "transcript_id")?

thanks much
Siva

Last edited by Siva; 04-28-2010 at 07:29 PM. Reason: errors
Siva is offline   Reply With Quote
Old 04-29-2010, 08:07 AM   #35
Siva
Member
 
Location: Ames, IA

Join Date: Nov 2008
Posts: 57
Default

Quote:
Originally Posted by Siva View Post
I have thought about this and I feel I am fundamentally wrong in using Cufflinks generated GTF output of one sample and the corresponding SAM file in htseqcount. Quoting from Cufflinks website; "Cufflinks constructs a parsimonious set of transcripts that explain the reads observed in an RNAseq expt." We know that it uses a pseudo ID to refer to various transcript bundles. If I use htseqcount using a SAM file against a cufflinks GTF I will only get read counts per pseudo ID. If I need read counts per gene or transcript i should have used a std reference annotation (GFF or GTF) to output counts per gene or transcript. However for many organisms including maize (see http://ftp.maizesequence.org/current/working-set/) the GFF3 format does not have the id attributes "gene_id" or "transcript id".

Am I right in assuming that I should not have used cufflinks generated GTF to output absolute counts per gene/transcript?

Can anyone help me with converting the std GFF to a GTF containing the necessary id attributes (viz., "gene_id" or "transcript_id")?

thanks much
Siva
Thanks Simon....as you said all I had to do was set --idattr=Parent in htseqcount and change the first column of the GFF file from 1 to chr1 and 2 to chr 2 etc..... It now seems to work!!!
more later
Siva
Siva is offline   Reply With Quote
Old 04-29-2010, 08:39 AM   #36
Cole Trapnell
Senior Member
 
Location: Boston, MA

Join Date: Nov 2008
Posts: 212
Default

Quote:
Originally Posted by Siva View Post
I have thought about this and I feel I am fundamentally wrong in using Cufflinks generated GTF output of one sample and the corresponding SAM file in htseqcount. Quoting from Cufflinks website; "Cufflinks constructs a parsimonious set of transcripts that explain the reads observed in an RNAseq expt." We know that it uses a pseudo ID to refer to various transcript bundles. If I use htseqcount using a SAM file against a cufflinks GTF I will only get read counts per pseudo ID. If I need read counts per gene or transcript i should have used a std reference annotation (GFF or GTF) to output counts per gene or transcript. However for many organisms including maize (see http://ftp.maizesequence.org/current/working-set/) the GFF3 format does not have the id attributes "gene_id" or "transcript id".

Am I right in assuming that I should not have used cufflinks generated GTF to output absolute counts per gene/transcript?

Can anyone help me with converting the std GFF to a GTF containing the necessary id attributes (viz., "gene_id" or "transcript_id")?

thanks much
Siva
Using cuffcompare on your individual Cufflinks sample assembles will produce a file called stdout.consensus.gtf, which will unify all of your transcripts under a common naming scheme.
Cole Trapnell is offline   Reply With Quote
Old 04-29-2010, 09:18 AM   #37
chrisbala
Member
 
Location: North Carolina

Join Date: Jan 2010
Posts: 82
Default

Hi Siva,

IN addition to Cole's answer, I have one other suggestion/idea.

Rather than predicting transcripts for each condition independently, I have piled all of my reads (2 conditions in my case) into one big file and used that for transcript prediction with cufflinks. The idea being that by doubling my data and improving my coverage, the predictions would be better. I don't know if this has any negative consequences for transcripts that might be unique to each condition?

Then I used the HT-Seq counter, and the mappings for each condition, to derive read counts for each condition/transcript.

Does anyone see any problems with this?
chrisbala is offline   Reply With Quote
Old 04-29-2010, 09:28 AM   #38
chrisbala
Member
 
Location: North Carolina

Join Date: Jan 2010
Posts: 82
Default GTF-GFF - HT-Seq

ON a different note, and relating to Siva's other question, for HT-Seq, does the GFF have to have + or - for the strand?

I've been trying to figure out what is wrong with my gff

Code:
chr1	taeGut1_ensGene	exon	8373168	8373328	.	+	.	gene_id=ENSTGUT00000007148.3;Name=;Parent=ENSTGUT00000007148
chr1	taeGut1_ensGene	exon	16249	16331	.	+	.	gene_id=ENSTGUT00000007148.2;Name=;Parent=ENSTGUT00000007148
chr1	CufflinksGene	exon	17859	18021	.	.	.	gene_id=all.3.1;Name=;Parent=all.3
line 1 reads fine, and is from a gff that works fine, line 2 only works if i have + (and not .) in the strand column.

line 3 gets an error
Code:
Feature all.3.1 does not contain a 'gene_id' attribute
Is the strand the problem?
chrisbala is offline   Reply With Quote
Old 04-29-2010, 09:30 AM   #39
chrisbala
Member
 
Location: North Carolina

Join Date: Jan 2010
Posts: 82
Default

yep, seems that strand needs to be +/- for HT-Seq, but Cufflinks produces some transcripts without strand info (which seems reasonable?)
chrisbala is offline   Reply With Quote
Old 04-29-2010, 09:39 AM   #40
Siva
Member
 
Location: Ames, IA

Join Date: Nov 2008
Posts: 57
Default

Quote:
Originally Posted by chrisbala View Post
yep, seems that strand needs to be +/- for HT-Seq, but Cufflinks produces some transcripts without strand info (which seems reasonable?)
Hi Chris
Yes, if you use GTF file from Cufflinks, you should set --stranded=no in htseqcount. Cufflinks does not give you strand information in all cases.

thanks
Siva
Siva is offline   Reply With Quote
Reply

Tags
normalization, rna-seq

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:08 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO