[chrisbala]

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

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

Feature all.3.1 does not contain a 'gene_id' attribute

Is the strand the problem?

[chrisbala]

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?

[Cole Trapnell]

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.

[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

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]

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

[Siva]

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

[Siva]

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

[steven]

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..

[lpachter]

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.

[mgogol]

one way to do it

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.

[Simon Anders]

Hi Siva

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

[Siva]

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

[Simon Anders]

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).

[lpachter]

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.

[Simon Anders]

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.