Announcement

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

  • DEXseq

    I have a 50X PE RNA-seq data and want to use DExseq. I understand that I need to get read counts from HTseq and normalize and count the expression. I have two questions:
    1. in case of no replicate how differential expression is estimated?
    2. For transcript expression do I have to go back and run cufflink first and then run DEseq, that will not make sense as cufflink will give me normalized count RPKM which cannot be used in DESEq. In summary how DExSeq estimate transcript expression. Any pointers will be great.


    Thanks

    Mathew

  • #2
    Dear Matthew,

    1. in case of no replicate how differential expression is estimated?

    The short answer is: no differential expression can be estimated without biological replicates. You can have a look around SEQAnswers about the necessity of estimating biological variability to determine differential expression.

    2. For transcript expression do I have to go back and run cufflink first and then run DEseq, that will not make sense as cufflink will give me normalized count RPKM which cannot be used in DESEq. In summary how DExSeq estimate transcript expression. Any pointers will be great.

    DEXSeq does not estimate transcript expression. It tests for differential usage of exons between conditions.

    For more details of discussion for both of your questions you can have a look into the DEXSeq paper:
    An international, peer-reviewed genome sciences journal featuring outstanding original research that offers novel insights into the biology of all organisms


    Hope it is useful!
    Alejandro
    Last edited by areyes; 10-22-2012, 01:31 AM. Reason: wording

    Comment


    • #3
      Hi!
      I'm trying to use the pipeline tophat2 -> HTSEq -> DEXSeq, but I'm stuck on the DEXSEq workflow!

      Here is what I did:

      1) I ran TopHat2 using a junc file:
      tophat2 -o tophat_out -i 10 -I 30000 -p 8 --library-type fr-unstranded -j combined.juncs -G my.gff --transcriptomeindex=transcriptome_data chr.fa myreads_1.fastq myreads_2.fastq
      2) I converted bam to the sam file:
      samtools sort -n accepted_hits.bam accepted_hits.nsorted
      samtools view accepted_hits.nsorted.bam > ./accepted_hits.nsorted.sam
      3) Then I produced the count table using htseq-count
      htseq-count -i ID accepted_hits.nsorted.sam my.gff -q --stranded=yes > count_HTSeq_exons.txt
      4)Now I'm trying to use the DEXSeq package to test differential use of exons

      #My code in R
      library("DEXSeq")
      samples = data.frame(
      condition = factor(c("sl","all","all","all","all")),
      replicate = factor(c(1:1,1:4)),
      type = c("single-read", "paired-end", "paired-end", "paired-end", "paired-end"),
      row.names = factor(c("cercaria_SL_Trapping_count_HTSeq_exons","Todos_cercaria_ERR022872_count_HTSeq_exons", "Todos_cercaria_ERR022875_count_HTSeq_exons", "Todos_cercaria_ERR022877_count_HTSeq_exons", "Todos_cercaria_ERR022878_count_HTSeq_exons")),
      stringsAsFactors = TRUE,
      check.names = FALSE
      )
      annotationfile = file.path("./GFF/v5.07.08.12.chado.HTSeq.gff")
      ecs = read.HTSeqCounts(countfiles = paste(rownames(samples), "txt", sep=".") , design = samples, flattenedfile = annotationfile)
      sampleNames(ecs) = rownames(samples)
      But then I get an error: Erro em x[[i]] : índice fora de limites

      If I try the last command without the GFF file, the object ecs is created without problems, what makes me think that the problem is with my GFF file, the same I used with the HTSeq-count.

      Here is a piece of my gff file:
      ##gff-version 3
      ##sequence-region Schisto_mansoni.Chr_1 1 65476681
      Schisto_mansoni.Chr_1 chado contig 1 19825 . + . ID=contig_14209;Name=null;
      Schisto_mansoni.Chr_1 chado gene 11159 12750 . + . ID=Smp_186980;
      Schisto_mansoni.Chr_1 chado exon 11159 11220 . + 0 ID=Smp_186980.1:exon:1;Parent=Smp_186980.1;
      Schisto_mansoni.Chr_1 chado exon 12411 12750 . + 0 ID=Smp_186980.1:exon:2;Parent=Smp_186980.1;
      Schisto_mansoni.Chr_1 chado mRNA 11159 12750 . + . ID=Smp_186980.1;Parent=Smp_186980;
      Schisto_mansoni.Chr_1 chado polypeptide 11159 12750 . + . ID=Smp_186980.1ep;Derives_from=Smp_186980.1;timelastmodified=14.10.2011+12:49:04+BST;feature_id=23195294;
      Schisto_mansoni.Chr_1 chado gene 16927 17315 . + . ID=Smp_197050;
      Schisto_mansoni.Chr_1 chado exon 16927 17082 . + 0 ID=Smp_197050.1:exon:1;Parent=Smp_197050.1;
      Schisto_mansoni.Chr_1 chado exon 17286 17315 . + 0 ID=Smp_197050.1:exon:2;Parent=Smp_197050.1;
      Schisto_mansoni.Chr_1 chado polypeptide 16927 17315 . + . ID=Smp_197050.1ep;Derives_from=Smp_197050.1;colour=2;timelastmodified=14.10.2011+12:50:09+BST;feature_id=23195325;
      Schisto_mansoni.Chr_1 chado mRNA 16927 17315 . + . ID=Smp_197050.1;Parent=Smp_197050;
      Schisto_mansoni.Chr_1 chado gap 19826 20025 . + . ID=Schisto_mansoni.Chr_1.embl:gap:19825-20025;
      Schisto_mansoni.Chr_1 chado contig 20026 60973 . + . ID=contig_14210;Name=null;
      Schisto_mansoni.Chr_1 chado gene 51849 114890 . + . ID=Smp_160500;
      Schisto_mansoni.Chr_1 chado exon 51849 51861 . + 0 ID=Smp_160500.1:exon:1;Parent=Smp_160500.1;
      Schisto_mansoni.Chr_1 chado exon 51901 51987 . + 0 ID=Smp_160500.1:exon:2;Parent=Smp_160500.1;
      Schisto_mansoni.Chr_1 chado exon 52024 52416 . + 0 ID=Smp_160500.1:exon:3;Parent=Smp_160500.1;
      Schisto_mansoni.Chr_1 chado exon 52448 52679 . + 0 ID=Smp_160500.1:exon:4;Parent=Smp_160500.1;
      Schisto_mansoni.Chr_1 chado exon 52715 52802 . + 0 ID=Smp_160500.1:exon:5;Parent=Smp_160500.1;

      Any sugestion??

      Thanks!

      Comment


      • #4
        Admittedly, the error message is not very helpful, but your problem is that you used htseq-count to get your count table. However, you are supposed to use the two Python scripts that come with DEXSeq instead. Please see the vignette (last section) for details.

        Comment


        • #5
          Hi
          I have used DEXSeq R package successfully to do alternative splicing detection (at least it gives output without error)
          My question is about the novel junction/splicing patterns. Since I notice DEXSeq uses annotation file, does DEXSeq detects novel splicing pattern, for example a transcript contains some of the intron, or any scenario it can be arise in the cancer abnormal genome?

          I did use TopHat with annotation, maybe in order to detect novel splicing junction I ought to use the de novo aligner to start with, instead of normal TopHat. Lastly now one can run TopHat 2 without annotation.

          Maybe my question is really an aligner question rather a DEXSeq question, in order to detect novel splicing one ought to use a de novo aligner. Any suggestion for a de novo aligner?
          Thanks
          Wenhong

          Comment


          • #6
            The counts script that comes with DEXseq that Dr. Anders referenced will produce counts per exon, which are then linked to annotated genes. No de novo detection is accomplished unless you provide a GTF file with de novo exons. Further, the junctions are not used, only the relative exon abundances are used.

            However, some people here have described positive results using the junctions output from TopHat2, as if they were exons. They get fewer overall reads (since reads wholly contained within an exon, and not spanning the splice boundary, are not used), but described getting consistently reliable, and verifiable results. And with 100-base PE reads, it seems more likely to get reads spanning junctions. I'm intending to try it soon, I'm interested in the outcome. It perhaps more directly addresses the question about differential isoforms, than using relative exon abundances which I feel only indirectly infers the alternative splicing. Still, for exons that go from zero to a zillion, it works great! :-)

            I have not created a GTF file from the junctions.bed output of TopHat2 -- I presume it's reasonable once you figure out the GTF details, linking junctions to parent genes. If you want to enable truly de novo work, you'd run CuffMerge (which runs CuffCompare) then use the merged.gtf file to define genes -- then link the junctions.bed entries together using those genes. That way if there is a completely novel gene locus, you'd be able to stitch together its junctions.bed entries appropriately. DEXseq will only test a gene if it has more than one child GTF feature, so you can't just feed it junctions.bed without grouping them by gene.

            Comment


            • #7
              Originally posted by jmw86069 View Post
              The counts script that comes with DEXseq that Dr. Anders referenced will produce counts per exon, which are then linked to annotated genes. No de novo detection is accomplished unless you provide a GTF file with de novo exons. Further, the junctions are not used, only the relative exon abundances are used.

              However, some people here have described positive results using the junctions output from TopHat2, as if they were exons. They get fewer overall reads (since reads wholly contained within an exon, and not spanning the splice boundary, are not used), but described getting consistently reliable, and verifiable results. And with 100-base PE reads, it seems more likely to get reads spanning junctions. I'm intending to try it soon, I'm interested in the outcome. It perhaps more directly addresses the question about differential isoforms, than using relative exon abundances which I feel only indirectly infers the alternative splicing. Still, for exons that go from zero to a zillion, it works great! :-)

              I have not created a GTF file from the junctions.bed output of TopHat2 -- I presume it's reasonable once you figure out the GTF details, linking junctions to parent genes. If you want to enable truly de novo work, you'd run CuffMerge (which runs CuffCompare) then use the merged.gtf file to define genes -- then link the junctions.bed entries together using those genes. That way if there is a completely novel gene locus, you'd be able to stitch together its junctions.bed entries appropriately. DEXseq will only test a gene if it has more than one child GTF feature, so you can't just feed it junctions.bed without grouping them by gene.

              Hi jmw86069. did you eventually try to run DEXseq with tophat junction.bed counts? I would be interested in knowing if it worked. I tried it with DESeq2, using each junction as separate identity and it sort of worked.

              Comment


              • #8
                Hi all, I just saw a poster in a conference in which they used the exon-exon junction counts from the STAR aligner as input to DEXSeq. It seemed to work nicely for human cancer genomes...

                Comment


                • #9
                  de novo?

                  Originally posted by areyes View Post
                  Hi all, I just saw a poster in a conference in which they used the exon-exon junction counts from the STAR aligner as input to DEXSeq. It seemed to work nicely for human cancer genomes...

                  Great! I never used STAR, but do you remember if they used de novo junctions detected from the data?

                  Maybe is just my data that is too funky. The dispersion estimates are completely out of whack:

                  Comment


                  • #10
                    I think STAR does looks for de novo junctions!

                    Yes, the plot look funky! There seems to be three clouds of points, the lowest one is what is found in the "normal" cases. The upper clouds are the funny ones... they even seem to have the opposite mean-variance relation to the lower cloud...

                    I am a bit out of ideas, maybe Mike or Simon could give more feedback... does the matrix count seem a bit sparsed? What kind of data is this?

                    Alejandro

                    Comment


                    • #11
                      Originally posted by areyes View Post
                      What kind of data is this?
                      The data comes from the junctions.bed files from tophat. I gave each junction a unique name (a combination of locus and junctions coordinates), created a matrix of counts by merging these junctions.bed for all my conditions and ran DESeq with it.

                      does the matrix count seem a bit sparsed?
                      A lot yes. I see that many junctions have counts for only 2 or 3 of the conditions. I do get ~400 significantly expressed junctions (FDR < 0.1), but these are *all* like:

                      c1 c1 c2 c2
                      0 0 100 80
                      900 850 0 0

                      Mock numbers but it gives you an idea. Other junctions with low FDR are similar.

                      So, yes, funky indeed. I could use only counts for annotated junctions, but since I am interested in global mis-regulation of splicing it would sort of defy the point.

                      Comment


                      • #12
                        It seems that your junctions fall into two classes, those which behave soemwhat normal and those which make up this strange upper cloud.

                        I guess the upper cloud is made by junctions with zero counts in all but one sample, and in this one sample, you get a lot of counts. The higher the counts in the one non-zero sample, the larger are both dispersion and mean, and this is why the cloud curves upwards rather than downwards.

                        Now, if you see counts only in a single sample, and not in all replicate samples of the same condition, you should not infer much from it. And if this happens a lot, this means that the differences between your samples are driven by something which does not go along sample groups (conditions), i.e., you have a huge non-reproducibilty in your experiment.

                        Maybe tell us more about your experimental design and biology.

                        Comment


                        • #13
                          If I was working with count data that produced the above plot, I would split out groups for the curve fitting and shrinkage of dispersion estimates, because the underlying assumption that the dispersions belong to a single cloud is violated.

                          Code:
                          dds= estimateSizeFactors(dds)
                          dds = estimateDispersionsGeneEst(dds)
                          # put here an empty vector for now
                          dispersions(dds) = numeric(nrow(dds))
                          
                          # define a cutoff by eye
                          cutoff = 1
                          upper = which(mcols(dds)$dispGeneEst > cutoff)
                          lower = which(mcols(dds)$dispGeneEst <= cutoff)
                          
                          ddsLower = dds[lower,]
                          ddsLower = estimateDispersionsFit(ddsLower)
                          ddsLower = estimateDispersionsMAP(ddsLower)
                          
                          # just use the MLE for the upper ones
                          dispersions(dds)[upper] = mcols(dds[upper,])$dispGeneEst
                          
                          # use the final estimates for the lower ones
                          dispersions(dds)[lower] = dispersions(ddsLower)
                          
                          # this is just necessary for the correct labeling in the plot
                          mcols(dds)$dispOutlier = FALSE
                          mcols(dds)$dispOutlier[mcols(dds)$dispersion == mcols(dds)$dispGeneEst] = TRUE
                          
                          plotDispEsts(dds)

                          Comment


                          • #14
                            Thanks for the reply Simon.

                            Originally posted by Simon Anders View Post
                            Now, if you see counts only in a single sample, and not in all replicate samples of the same condition, you should not infer much from it. And if this happens a lot, this means that the differences between your samples are driven by something which does not go along sample groups (conditions), i.e., you have a huge non-reproducibilty in your experiment.
                            Yep, this seems to be true. The samples are not mine, I am just the guy trying to make sense of the data.

                            Originally posted by Simon Anders View Post
                            Maybe tell us more about your experimental design and biology.
                            In brief, experimentally, my collaborators are trying to get a very small number of cells from a very defined population. This involves a very long protocol in the bench, and my opinion the longer it is the more prone to errors and sample degradation, and in the end, they need to pool several samples together. So the replicates are more technical then biological. To solve this issue, they are now changing the experimental design. This will involve fewer steps, but a more heterogeneous tissue composition.

                            Their hypothesis does make sense though. They are looking at a protein expressed only a in cell niche, that interacts with splicing factors. Together with other evidence, they hypothesised that it might be regulating AS (or mis-regulating it in the kd mice).

                            Comment


                            • #15
                              Well, if they see such strong differences between samples that are essentially technical replicates, then they are probably indeed well advised to try to improve their bench protocol.

                              Comment

                              Latest Articles

                              Collapse

                              • seqadmin
                                Advanced Tools Transforming the Field of Cytogenomics
                                by seqadmin


                                At the intersection of cytogenetics and genomics lies the exciting field of cytogenomics. It focuses on studying chromosomes at a molecular scale, involving techniques that analyze either the whole genome or particular DNA sequences to examine variations in structure and behavior at the chromosomal or subchromosomal level. By integrating cytogenetic techniques with genomic analysis, researchers can effectively investigate chromosomal abnormalities related to diseases, particularly...
                                09-26-2023, 06:26 AM
                              • seqadmin
                                How RNA-Seq is Transforming Cancer Studies
                                by seqadmin



                                Cancer research has been transformed through numerous molecular techniques, with RNA sequencing (RNA-seq) playing a crucial role in understanding the complexity of the disease. Maša Ivin, Ph.D., Scientific Writer at Lexogen, and Yvonne Goepel Ph.D., Product Manager at Lexogen, remarked that “The high-throughput nature of RNA-seq allows for rapid profiling and deep exploration of the transcriptome.” They emphasized its indispensable role in cancer research, aiding in biomarker...
                                09-07-2023, 11:15 PM
                              • seqadmin
                                Methods for Investigating the Transcriptome
                                by seqadmin




                                Ribonucleic acid (RNA) represents a range of diverse molecules that play a crucial role in many cellular processes. From serving as a protein template to regulating genes, the complex processes involving RNA make it a focal point of study for many scientists. This article will spotlight various methods scientists have developed to investigate different RNA subtypes and the broader transcriptome.

                                Whole Transcriptome RNA-seq
                                Whole transcriptome sequencing...
                                08-31-2023, 11:07 AM

                              ad_right_rmr

                              Collapse

                              News

                              Collapse

                              Topics Statistics Last Post
                              Started by seqadmin, 09-27-2023, 06:57 AM
                              0 responses
                              10 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 09-26-2023, 07:53 AM
                              0 responses
                              10 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 09-25-2023, 07:42 AM
                              0 responses
                              15 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 09-22-2023, 09:05 AM
                              0 responses
                              45 views
                              0 likes
                              Last Post seqadmin  
                              Working...
                              X