Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • #16
    Originally posted by Simon Anders View Post
    Please take a few of these many reads in location 2, hover your mouse over them in IGV to get the read name, then grep for the read name to find the corresponding lines in the SAM file. Maybe they are all non-unique mapping reads.

    The fact that htseq-count does not count these reads is deliberate, of course. See earlier threads for an explanation.
    Many thanks Simon! I think I need to get my basic understanding of mapping solid.

    Forgive my naivety, If multi-mapping becomes an issue when we are translating mapping reads to count data (expression level)...what than would be an acceptable standard to extract count data?? (is it using htseq-count ,bedtools, cuffdiff)

    Also....during the mapping stage..what bias would be introduced, if we discard multimap reads?


    Edit - you are correct in saying that the reads at location 2 is multimapped, many thanks for the insight!

    Comment


    • #17
      See my post #4 in this thread: http://seqanswers.com/forums/showthread.php?t=9129

      Comment


      • #18
        Dear Zaki,

        You ask an excellent question. This paper answers it in detail:

        Lior

        Comment


        • #19
          Thanks Simon and Lior for the suggested readings!

          Taking a step back and reading the fundamentals behind each approach is really essential. I think its something I overlooked at first.

          Since the protocol and vignette of each approach is very well described and easy to follow, I guess I jumped the gun and started running commands without understanding much behind it (I think I still have a lot yet to understand... statistics...mapping..biology...etc...)

          Thanks again for your comments and direction! It helped a lot.



          Zaki
          Last edited by zaki; 04-01-2013, 12:47 AM.

          Comment


          • #20
            Originally posted by Simon Anders View Post
            Please take a few of these many reads in location 2, hover your mouse over them in IGV to get the read name, then grep for the read name to find the corresponding lines in the SAM file. Maybe they are all non-unique mapping reads.

            The fact that htseq-count does not count these reads is deliberate, of course. See earlier threads for an explanation.
            Indeed. I was having the same problem. I was getting crazy because I had already checked what Simon pointed about the possibility of having non-unique mapping reads. I think I have found out why.
            If you use htseq-count with -i gene_id, it will construct the set where to add up reads based on the ENSEMBL GENE ID, which will sum reads that could come from different ENSEMBL transcript ID. You may lose the information about transcripts at this level (maybe you can use next DEXSeq approach to look inside...) but at least, you do not lose reads mapping to different identifiers(transcript id), as was my case because I used the UCSC gtf file where gene_id=trascritpt_id.

            However, I still encountered IGV regions in my BAM files with many reads that htseq count anyhow flagged as ambigous, but were not multi aligned. I think I have found out why by looking at the locus with the ensembl genome browser and looking for the "gene_name", that is, the SYMBOL. I realized that the Ensembl annotation (present in the GTF) contained for the same locus, with the same gene_name (SYMBOL), two different ENSEMBLGENE Ids for what it called: "protein coding" and "non-sense-mediated decay", as gene biotype. Both lying in the same locus with high overlap. I guessed then that htseq-count had found reads aligning to both ENSEMBLGENE identifiers and discarded them as ambiguous.

            While discarding reads mapping to different locus or different genes which overlap (hence also different gene Symbols) may have sense, I do not see the sense in discarding reads merely because we have a maybe too exhaustive overlapping annotation. We might also have other features like this, antisense ncRNAs, pseudogenes, for instance, which may have a big overlapping with its parental, and hence, all counts in regions with the same sequence will be discarded.
            For the pseudogenes, unfortunately, if they have a huge overlap, all those reads will be discarded by htseq-count -i gene_name either, but if we are not interested in looking at pseudogenes we could choose to use a more limited annotation like RefSeq.

            So, what do you think of using as a global strategy with htseq-count when interested in having also all non coding, pseudogenes et. al a command line like this:

            htseq-count -m union -t exon -i gene_name - $GTF > $out

            It is always possible to try and deconvolve counts thereafter,
            What do you think about it?
            Thanks
            Best
            Jose
            Last edited by Jose Garcia; 07-17-2013, 07:28 AM.

            Comment


            • #21
              In order to get gene level expression which is a better tool to use? HTSeq or Bedtools. Or has anyone tried any other tool for this purpose?

              Comment


              • #22
                Don't use bedtools, the results will be wrong for what you want to do. Use htseq-count or featureCounts (from subread).

                Comment


                • #23
                  Thank you for your reply.
                  For gene expression measurement, htseq values I believe could be conservative as only unique mapping events are used along with non-overalapping reads. Feature counts is fast and through there is substantial concordance of count values with htseq count. About 20 features per sample have count differences of greater than 10000.

                  Which one does one believe?

                  Comment


                  • #24
                    As a general rule, htseq-count. You might use the debug options ("-o" for htseq-count, I don't recall what it is for featureCounts) to see why there's such a big difference between the two. Given the same thresholds (i.e., MAPQ filtering) and intersection model, the counts should be the same.

                    Comment


                    • #25
                      Originally posted by vyellapa View Post
                      Thank you for your reply.
                      For gene expression measurement, htseq values I believe could be conservative as only unique mapping events are used along with non-overalapping reads. Feature counts is fast and through there is substantial concordance of count values with htseq count. About 20 features per sample have count differences of greater than 10000.

                      Which one does one believe?
                      The following paper explains in details the difference between featureCounts and htseq-count:

                      http://www.ncbi.nlm.nih.gov/pubmed/24227677

                      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
                      10 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, Yesterday, 06:07 PM
                      0 responses
                      10 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 03-22-2024, 10:03 AM
                      0 responses
                      51 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