RNA-seq and normalization numbers

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

[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

[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]

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]

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

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

[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?

[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]

yep, seems that strand needs to be +/- for HT-Seq, but Cufflinks produces some transcripts without strand info (which seems reasonable?)

[Siva]

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

[cek]

I order to use DESeq to test for differential expression between my RNAseq conditions, I calculate raw read counts by transcript using cufflinks output with the following formulae (as proposed by RockChalkJayhawk) :
raw1 = FPKM * length (kb) * number of mapped reads (million).

However, on another seqanswers post (http://seqanswers.com/forums/showthr...links+coverage), Cole Trapnell suggests to calculate raw read counts like this :
raw2 = coverage * length (from transcripts.expr file)

However, these two calculations do not lead to the same result. Have someone notice the same difference in their data ?

[answersseq]

I do learn a lot from discussion here.
Any opinions on microRNA sequencing data? Their length is similar but many reads can be mapped to multiple locations (or mature miRs).
To compare differential expression between cell lines, tissue, I guess we would expect big difference, as well as no house keeping miRs..

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

[hypatia]

normalization with all or uniquely map reads

Hi Zee
I was wondering If you got the answer to this question. Is it 3??
should I elminate the unaligned or ambiguous maps out of the normalization?


"I've read about people doing counts as reads per million and log transforming these values to fit Poisson distribution, but it's sprung multiple ideas in my mind. Would this be as simple as dividing my counts for each experiment by
1) 1 Million
2) the total number of reads sequenced
3) the total number of uniquely mapped reads

I'm inclined to option (3) because that represents the amount of usable sequence data."

[zee]

I would go with uniquely mapped reads because it's a more accurate representation of how much sequence data you obtained from your runs.
You could get a bit more stringent by using Picard to filter out possible PCR duplicates from the alignments in BAM format.

[syambmed]

Hi guys,

I have trancriptome data from Illumina and am using CLC Genomic workbench for data analysis. I dont know or not familiar with other programs for transcriptome analysis. the data are from 1 sample of control cells and 1 sample of treated cells (no replicate for each sample) and I am looking for differently express genes.

The problem is normalization step. There are 3 types of normalization method offered by the software 1) scaling [option for normalization value= mean or median, baseline = median mean or median median] 2) quantile 3) total reads per 1million.

I dont know which one to choose..T_T Help me..

Then there are statistical tests on Gaussion data or on proportions. How to know that my data is suitable for which test..? I read that mostly people use Baggerley's.

A thing with Baggerley test is that the test outcome have p-value and false discovery rate (FDR) p-value correction. which one is used for determining differentially express genes..?

Thank you.