Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • #31
    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..

    Comment


    • #32
      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, 10:15 AM. Reason: typo

      Comment


      • #33
        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, 01:17 PM. Reason: typing errors

        Comment


        • #34
          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, 07:29 PM. Reason: errors

          Comment


          • #35
            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

            Comment


            • #36
              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.

              Comment


              • #37
                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?

                Comment


                • #38
                  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?

                  Comment


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

                    Comment


                    • #40
                      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

                      Comment


                      • #41
                        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 ?
                        Last edited by cek; 05-12-2010, 02:18 AM.

                        Comment


                        • #42
                          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..


                          Originally posted by Simon Anders View Post
                          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

                          Comment


                          • #43
                            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."

                            Comment


                            • #44
                              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.

                              Comment


                              • #45
                                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.

                                Comment

                                Latest Articles

                                Collapse

                                • seqadmin
                                  Strategies for Sequencing Challenging Samples
                                  by seqadmin


                                  Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                                  03-22-2024, 06:39 AM
                                • seqadmin
                                  Techniques and Challenges in Conservation Genomics
                                  by seqadmin



                                  The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                                  Avian Conservation
                                  Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                                  03-08-2024, 10:41 AM

                                ad_right_rmr

                                Collapse

                                News

                                Collapse

                                Topics Statistics Last Post
                                Started by seqadmin, Yesterday, 06:37 PM
                                0 responses
                                8 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, Yesterday, 06:07 PM
                                0 responses
                                8 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 03-22-2024, 10:03 AM
                                0 responses
                                49 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 03-21-2024, 07:32 AM
                                0 responses
                                66 views
                                0 likes
                                Last Post seqadmin  
                                Working...
                                X