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
                      Exploring the Dynamics of the Tumor Microenvironment
                      by seqadmin




                      The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
                      07-08-2024, 03:19 PM
                    • seqadmin
                      Exploring Human Diversity Through Large-Scale Omics
                      by seqadmin


                      In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
                      06-25-2024, 06:43 AM

                    ad_right_rmr

                    Collapse

                    News

                    Collapse

                    Topics Statistics Last Post
                    Started by seqadmin, Today, 11:09 AM
                    0 responses
                    16 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 07-19-2024, 07:20 AM
                    0 responses
                    148 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 07-16-2024, 05:49 AM
                    0 responses
                    124 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 07-15-2024, 06:53 AM
                    0 responses
                    111 views
                    0 likes
                    Last Post seqadmin  
                    Working...
                    X