No announcement yet.
  • 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=estimateExonFoldChanges(dxd, fitExpToVar="condition")



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


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



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

        ## 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")


        • #5
          you could try dropbox or similar tools!


          • #6
            is it solved?

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


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


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



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


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


                    Latest Articles


                    • 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...
                      Yesterday, 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





                    Topics Statistics Last Post
                    Started by seqadmin, Today, 06:57 AM
                    0 responses
                    1 view
                    Last Post seqadmin  
                    Started by seqadmin, Yesterday, 07:53 AM
                    0 responses
                    Last Post seqadmin  
                    Started by seqadmin, 09-25-2023, 07:42 AM
                    0 responses
                    Last Post seqadmin  
                    Started by seqadmin, 09-22-2023, 09:05 AM
                    0 responses
                    Last Post seqadmin