Seqanswers Leaderboard Ad



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

  • DESeq2 Anomaly/Issue?


    I am trying to identify patterns of gene expression changes (to group genes by their regulating TFs) in different genetic backgrounds. However I am getting one (at least) anomaly that I can't explain.

    This gene has one transcript. here are the cufflinks generated (from mapped BAM files of biological replicates) mean FPKM values:

    Wt: 139.2 (+/- 16.0) n = 4
    GOF: 318.9 (+/- 60.9) n = 4
    LOF: 6.8 (+/- 1.6) n = 4
    2x LOF: 11.7 (+/-16.3) n = 4

    It is upregulated when it should be (GOF) and off when it should be (LOF, 2x LOF), a pretty good example of an activated target gene for this group. (yes I know there seems to be an outlier replicate in 2x LOF)

    Yet every DESeq2 comparison gives NA across the board, because it reports the base mean expression is 0????

    Is there some problem with annotation of this gene going from the mapped read BAM files to HTSeqCount output file to DESeq2?

    I am working with a very well annotated genome - C. elegans ce10 version.

    Also side question: what exactly is the base mean expression/how is it calculated? is it the mean expression between the groups you are comparing?

    Thanks for any help!
    Last edited by The_Corsair; 02-27-2015, 08:43 PM.

  • #2
    Ok so I have been doing more digging.

    The problem seems to be when going from the mapped read BAM files to the HTSeqCount output file.

    I checked the HTSeqCount output files for the 4 Wt replicates and 4 GOF replicates and every single one has 0 counts for this gene.

    Can anyone tell me what is going on here? Why is HTSeqCount (or cuff links) miscounting the reads for this gene??


    • #3
      The only thing I can think of is it has too much sequence similarity to other genes but of the 5 related genes in the genome, all have counts.


      • #4
        It sounds like an HTSeqCount issue, my first guess would be that the GTF/gff file contains overlapping gene regions. Quick thing to check, view the gene models file directly or look at it in a genome browser in that region of interest. Similarly, look at the Tophat2 coverage in this region, you should be able to see whether there is read coverage of zero or not.

        Apologies if you checked this line of thinking already...

        Perhaps subject to some debate, HTSeqCount chooses not to report counts in regions where multiple genes overlap, by design. In theory, after removing a portion of a gene, you'd still be able to detect differential changes using the remaining non-overlapping region. But sometimes gene models are completely overlapping (intended with other assumptions in mind) and so some genes potentially are completely lost during the count step. It isn't really for HTSeqCount to correct issues in gene annotations, but something to keep in mind. Easiest fix is to change the gtf file manually to remove overlapping regions (or merge them into one larger gene locus) and regenerate counts.


        • #5
          DEXseq has a script to help with this ( See Section 2.4:
          The Bioconductor project aims to develop and share open source software for precise and repeatable analysis of biological data. We foster an inclusive and collaborative community of developers and data scientists.

          I've also modified htseq-count to optionally include ambiguous reads (ie those that overlap multiple entries in your GTF file) - happy to provide the source if interested.


          • #6
            fanli, many thanks! I would be interested in the modified htseq-count just for the sake of adding another tool to the bag of tricks.

            I respect the intent of not counting those regions, and add caution that including counts could cause other cascade issues. So be aware.

            I have used that DEXseq script, I find it very useful. In fact, a good workflow may be to run that DEXseq script for the purpose of identifying any regions with overlapping genes. Then update the GTF file to resolve those regions in some way. Finally, run htseq-count using the updated GTF/GFF3 file which contains no overlapping regions.

            The GenomicFeatures R package has a function disjointExons() which accomplishes a very similar task, for what it's worth.

            All this, and we don't really know if this addresses the OP. Let us know!


            • #7
              Attached. See the --allow_ambiguous option. And second jmw's point about being cautious when including these counts
              Attached Files


              • #8
                Thanks for the suggestions!

                There are reads for the gene but not many and they must all be mapping to multiple genes, since there are a total of six related genes that should all be expressed in wild-type. The sequencing depth was lower because we had so many replicates and different samples - cost affecting parameters. I think this gene was just expressed at a low level and I didn't pick up any unique reads for it.

                I think I won't bother trying anything else with this gene. It's not needed for my analysis, was just curious as to why of the 6 it was the only one not expressed. Thanks for the helpful tips!


                Latest Articles


                • seqadmin
                  The Impact of AI in Genomic Medicine
                  by seqadmin

                  Artificial intelligence (AI) has evolved from a futuristic vision to a mainstream technology, highlighted by the introduction of tools like OpenAI's ChatGPT and Google's Gemini. In recent years, AI has become increasingly integrated into the field of genomics. This integration has enabled new scientific discoveries while simultaneously raising important ethical questions1. Interviews with two researchers at the center of this intersection provide insightful perspectives into...
                  02-26-2024, 02:07 PM
                • seqadmin
                  Multiomics Techniques Advancing Disease Research
                  by seqadmin

                  New and advanced multiomics tools and technologies have opened new avenues of research and markedly enhanced various disciplines such as disease research and precision medicine1. The practice of merging diverse data from various ‘omes increasingly provides a more holistic understanding of biological systems. As Maddison Masaeli, Co-Founder and CEO at Deepcell, aptly noted, “You can't explain biology in its complex form with one modality.”

                  A major leap in the field has
                  02-08-2024, 06:33 AM





                Topics Statistics Last Post
                Started by seqadmin, 02-28-2024, 06:12 AM
                0 responses
                Last Post seqadmin  
                Started by seqadmin, 02-23-2024, 04:11 PM
                0 responses
                Last Post seqadmin  
                Started by seqadmin, 02-21-2024, 08:52 AM
                0 responses
                Last Post seqadmin  
                Started by seqadmin, 02-20-2024, 08:57 AM
                0 responses
                Last Post seqadmin