Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • dexseq results understanding

    Hey, Now I am trying to find splicing differences between two tumor and two normal samples(so only one condition: N and T factors), and prepare data following your vignette, everything goes well, but after the last step(bellow), I got the dxr dataframe I find amost half lines with exon usage coefficient as NA, when I extract the corresponding gene I find the exons well covered (one example bellow), so, my question is why so many NAs? could anyone please give some possible explanations which I would check.
    Many Many Thanks.

    ######
    my processes:
    dxd = DEXSeqDataSetFromHTSeq(countFiles, sampleData=sampleTable,design=~sample + exon + condition:exon, flattenedfile=flattenedFile)
    dxd=estimateSizeFactors(dxd)
    dxd=estimateDispersions(dxd)
    dxd=testForDEU(dxd)
    dxd=estimateExonFoldChanges(dxd, fitExpToVar="condition")
    dxr=DEXSeqResults(dxd)

    example:

    head(dxr)

    LRT p-value: full vs reduced

    DataFrame with 6 rows and 13 columns
    groupID featureID exonBaseMean dispersion stat pvalue padj R S log2fold_R_S genomicData countData transcripts
    <character> <character> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <GRanges> <matrix> <list>
    ENSG00000000003:E001 ENSG00000000003 E001 259.50 0.005893988 0.7015350362 0.4022684 1 NA NA NA chrX:-:[99883667, 99884983] 199 201 471 ... ########
    ENSG00000000003:E002 ENSG00000000003 E002 112.00 0.001052854 1.7011996938 0.1921312 1 NA NA NA chrX:-:[99885756, 99885863] 98 96 183 ... ########
    ENSG00000000003:E003 ENSG00000000003 E003 92.75 0.005572201 0.0007320773 0.9784143 1 NA NA NA chrX:-:[99887482, 99887537] 97 76 145 ... ########
    ENSG00000000003:E004 ENSG00000000003 E004 71.50 0.010739013 0.0385784467 0.8442862 1 NA NA NA chrX:-:[99887538, 99887565] 76 64 112 ... ########
    ENSG00000000003:E005 ENSG00000000003 E005 81.00 0.007045578 0.7014820356 0.4022862 1 NA NA NA chrX:-:[99888402, 99888438] 90 61 131 ... ########
    ENSG00000000003:E006 ENSG00000000003 E006 111.75 0.002231303 1.9566454114 0.1618725 1 NA NA NA chrX:-:[99888439, 99888536] 117 80 192 ... ########

    exon counts:
    head(counts(dxd))
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
    ENSG00000000003:E001 199 201 471 167 855 737 1587 511
    ENSG00000000003:E002 98 96 183 71 956 842 1875 607
    ENSG00000000003:E003 97 76 145 53 957 862 1913 625
    ENSG00000000003:E004 76 64 112 34 978 874 1946 644
    ENSG00000000003:E005 90 61 131 42 964 877 1927 636
    ENSG00000000003:E006 117 80 192 58 937 858 1866 620

  • #2
    I'm having a similar issue where the log2fold value for many exons is also NA, while they show sufficient counts in the data. Did you find a solution/explanation for this behavior?

    Comment


    • #3
      Thanks for reporting this! The current implementation of DEXSeq estimates fold changes only for those exons for which a p-adjustes value was calculated (e.g. not being an outlier).

      How many of such cases of having a p-value but not a fold change do you see in your data? If they are many, could you maybe send me your object so I could have a closer look to what is happening?

      Alejandro

      Comment


      • #4
        Thanks for your response.

        Out of the total 716 genes that are significant (FDR < 0.05), 101 have NA for the log2fold column, but they all have an adjusted p-value.

        I am not able to attach the .Rdata object of the dxd variable created using the documentation as the file is too large. Can I send it to you in any other way or is there a more specific object that would be useful for you to debug what is going on? I generated the dxd object using the following code.

        Code:
        ## makeTranscriptDbFromGFF
        gffFile <- makeTranscriptDbFromGFF("/Genome_files/Mus_musculus/UCSC/mm10/Annotation/Genes/genes.gtf", format="gtf")
        
        ## preparing exonic parts
        exonicParts <- disjointExons(gffFile, by="exon", aggregateGenes=FALSE)
        
        
        align <- "/Samples"
        
        files <- list.files(path=align, pattern="*.bam", full.names=T, recursive=FALSE)
        
        bf1 <- BamFileList(c(files),index=character(), asMates=TRUE)
        
        genehits <- summarizeOverlaps(exonicParts, bf1, mode="IntersectionStrict", ignore.strand=FALSE, singleEnd=FALSE, inter.feature=TRUE, fragments=TRUE)
        
        colData(genehits)$condition <- c("WT", "MT", "WT", "MT", "WT", "MT")
        raw_dat <- assays(genehits)$counts
        
        conds <- c("WT1", "MT1", "WT2", "MT2", "WT3", "MT3")
        colnames(raw_dat) <- conds
        ## reorder columns
        raw_data <- cbind(raw_dat[,c(1,3,5)], raw_dat[,c(2,4,6)]) 
        
        
        geneID <- exonicParts$gene_id #-> contains gene id
        g <- unlist(geneID)
        exonID <- exonicParts$exonic_part # contains exon number
        
        nam <- paste(g,exonID, sep=":")
        
        rownames(raw_dat) <- nam
        
        
        dxd <- DEXSeqDataSet(raw_data,sampleTable, design, featureID= as.character(exonID), groupID= g)
        
        dxd <- estimateSizeFactors(dxd)
        dxd <- estimateDispersions(dxd)
        
        dxd <- testForDEU(dxd)
        dxd <- estimateExonFoldChanges(dxd, fitExpToVar="condition")
        dxr1 <- DEXSeqResults( dxd )
        
        save(dxd, file="DEXseq_v1.Rdata")

        Comment


        • #5
          you could try dropbox or similar tools!

          Comment


          • #6
            is it solved?

            I have the same problem. I wondering if the problem was solved, if yes, how? could anybody reply?
            Thanks!

            Comment


            • #7
              Just wondering if this has ever been solved -I'm having the same issue..

              Comment


              • #8
                Sorry for the very late reply. Could you confirm if this happens when you have a moderate number of samples (around 10-12) and for genes with lots of exonic regions?

                I think this has to do with the GLMs that are calculated for each gene to estimate the exon fold changes. The way this works is that a model frame is created where the rows is number of exonic regions times number of samples. I set up a threshold in 3000 to the number of rows of the model frame such the fit would not be done if the model frame for a gene passes this threshold. I added an option in the latest development version such that users can increase this number, but consider that the larger this value, the larger it will take to compute. The parameter is the maxRowsMF of the function estimateExonFoldChanges.

                Consider that this is a temporary solution, we need to think in a smarter way to deal with this!

                Alejandro

                Comment


                • #9
                  Thanks for the the update-ill give the development version a go, I have over 20 samples and the analysis is on human genes.

                  Comment


                  • #10
                    Can DEXSeq handle more than 70 samples to analyse? I have different human tissue to analyse which not same as normal vs disease as shown in vignette.

                    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
                    67 views
                    0 likes
                    Last Post seqadmin  
                    Working...
                    X