Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • You are supposed to pass first the SAM file, then the GFF file.

    Comment


    • End location of read in HTSeq-count

      Hi Simon,

      I am a user of the HTSeq-count.
      It's a wonderful tool, and I really love it.

      Nevertheless, I have a small question about it.
      I am wondering how dose HTSeq-count define the end location (right site end) of a read from RNA-seq sam file?
      I am not very familiar with Pathon, so I cannot well understand the source code.

      For single-end strand-specific RNA-seq sam file, is the end location defined by the CIGAR?
      That is, End= Start + ( No. of Ms) + (No. of Ns) + (No. of Ds) ?

      For paired-end strand-specific RNA-seq sam file, how does HTSeq-count handle this?
      Is the start location (left site) defined as the smaller one of the two start loci within a pair?
      And the end location (right site end) defined as the larger one of the two end sites within a pair?


      Many thanks and best regards,
      Jerry

      Comment


      • Originally posted by sumanex
        However, the following command is not showing any coverage only the coordinates and a blue horizotal line at 0.
        >>> import matplotlib
        >>> from matplotlib import pyplot
        >>> pyplot.plot( list( cvg[ HTSeq.GenomicInterval( "chr2", 2000000, 3000000, "+" ) ] ) )
        [<matplotlib.lines.Line2D object at 0x104751ac>]
        >>> pyplot.show()
        If this is all you have, it's not surprising that there's no coverage information, because you haven't put any in. The GenomicInterval object defines a region of a chromosome.

        Do have to provide a gtf file also for showing the coverage.
        GTF files are also region definition objects, rather than coverage objects. You need a GenomicArray in there, together with something that populates the array. The HTSeq documentation has good code examples. You need to get your head around what's going on before diving into code:

        http://www-huber.embl.de/users/ander...-full-coverage

        Comment


        • matplot coverage

          Hi,
          I had to change the coordinates to chr1, 133363801,134605170

          of cvg vector to get a decent plot.

          thanks
          suman

          Comment


          • HtSeq TSS plots

            Hi,
            I am following the TSS plots code in HTSeq in ubunut 12.04 but getting error:

            1) in step p in tsspos answer is false. but getting a wincvg plot in the matplot window. the gtf file is for mouse release64 ensembl.

            Python 2.7.3 (default, Aug 1 2012, 05:16:07)
            [GCC 4.6.3] on linux2
            Type "help", "copyright", "credits" or "license" for more information.
            >>> import HTSeq
            >>> bamfile = HTSeq.BAM_Reader( "galaxy58.bam" )
            >>> bamfile
            <HTSeq.BAM_Reader object at 0xb721afac>
            >>> gtffile = HTSeq.GFF_Reader( "Mus_musculus.NCBIM37.64.gtf" )
            >>> coverage = HTSeq.GenomicArray( "auto", stranded=False, typecode="i" )
            >>> coverage
            <HTSeq._HTSeq.GenomicArray object at 0x9065b0c>
            >>> for almnt in bamfile:
            ... if almnt.aligned:
            ... coverage[ almnt.iv ] += 1
            ...
            >>> import itertools
            >>> for feature in itertools.islice( gtffile, 100):
            ... if feature.type == "exon" and feature.attr["exon_number"] == "1":
            ... print feature.attr["gene_id"], feature.attr["transcript_id"], feature.iv.start_d_as_pos
            ...
            ENSMUSG00000000702 ENSMUST00000105216 NT_166433:11954/+
            ENSMUSG00000078423 ENSMUST00000105217 NT_166433:28586/+
            ENSMUSG00000078424 ENSMUST00000105218 NT_166433:47744/+
            ENSMUSG00000078424 ENSMUST00000105219 NT_166433:47927/+
            ENSMUSG00000071964 ENSMUST00000096694 NT_166402:12542/+
            ENSMUSG00000091539 ENSMUST00000165255 18:3123464/-
            ENSMUSG00000063889 ENSMUST00000151311 18:3327588/-
            >>> tsspos = set()
            >>> for feature in gtffile:
            ... if feature.type == "exon" and feature.attr["exon_number"] == "1":
            ... tsspos.add( feature.iv.start_d_as_pos )
            ...

            >>>
            >>> p = HTSeq.GenomicPosition( "chr1", 133363801, "+" )
            >>> p
            <GenomicPosition object 'chr1':133363801, strand '+'>
            >>> p in tsspos
            False


            >>> halfwinwidth = 3000
            >>> window = HTSeq.GenomicInterval( p.chrom, p.pos - halfwinwidth, p.pos + halfwinwidth, "." )
            >>> window
            <GenomicInterval object 'chr1', [133360801,133366801), strand '.'>
            >>> list( coverage[window] )
            >>> import numpy
            >>> wincvg = numpy.fromiter( coverage[window], dtype='i', count=2*halfwinwidth )
            >>> wincvg
            array([0, 0, 0, ..., 0, 0, 0])
            >>> from matplotlib import pyplot
            >>> pyplot.plot( wincvg )
            [<matplotlib.lines.Line2D object at 0x10c8122c>]
            >>> pyplot.show()


            and 2) Index error start too small


            >>> profile = numpy.zeros( 2*halfwinwidth, dtype='i' )
            >>> profile
            array([0, 0, 0, ..., 0, 0, 0])
            ^

            >>> for p in tsspos:
            ... window = HTSeq.GenomicInterval( p.chrom, p.pos - halfwinwidth, p.pos + halfwinwidth, "." )
            ... wincvg = numpy.fromiter( coverage[window], dtype='i', count=2*halfwinwidth )
            ... if p.strand == "+":
            ... profile += wincvg
            ... else:
            ... profile += wincvg[::-1]
            ...
            Traceback (most recent call last):
            File "<stdin>", line 3, in <module>
            File "_HTSeq.pyx", line 520, in HTSeq._HTSeq.GenomicArray.__getitem__ (src/_HTSeq.c:10593)
            File "_HTSeq.pyx", line 370, in HTSeq._HTSeq.ChromVector.__getitem__ (src/_HTSeq.c:7538)
            IndexError: start too small

            any help is appreciated
            suman

            Comment


            • Originally posted by sumanex View Post
              >>> tsspos = set()
              >>> for feature in gtffile:
              ... if feature.type == "exon" and feature.attr["exon_number"] == "1":
              ... tsspos.add( feature.iv.start_d_as_pos )
              ...
              >>> p = HTSeq.GenomicPosition( "chr1", 133363801, "+" )
              >>> p
              <GenomicPosition object 'chr1':133363801, strand '+'>
              >>> p in tsspos
              False
              You've printed the type of p, but not the type of feature.iv.start_d_as_pos -- I find this omission a bit strange. The 'in' operator checks for membership (I think with a deep equality check), and it may be that the precise GenomicPosition object that you've declared is not in that set. Perhaps there's a position HTSeq.GenomicPosition( "chr1", 133363801, "."), or a position HTSeq.GenomicPosition( "chr1", 133363800, "+"). It's hard to know for sure without this detail.

              >>> list( coverage[window] )
              >>> import numpy
              ...
              >>> for p in tsspos:
              ... window = HTSeq.GenomicInterval( p.chrom, p.pos - halfwinwidth, p.pos + halfwinwidth, "." )
              ... wincvg = numpy.fromiter( coverage[window], dtype='i', count=2*halfwinwidth )
              ... if p.strand == "+":
              ... profile += wincvg
              ... else:
              ... profile += wincvg[::-1]
              ...
              Traceback (most recent call last):
              ...
              IndexError: start too small
              I guess it's possibly that wincvg and/or coverage[window] is empty, or the numpy.fromiter count value is too high for the size of coverage[window]. Could you print the first few entries of tsspos?

              Comment


              • tsspos

                the lines of tsspos are as follows:

                <GenomicPosition object '2':49762399, strand '-'>, <GenomicPosition object '2':26280379, strand '-'>, <GenomicPosition object '6':42236683, strand '+'>, <GenomicPosition object '6':47528008, strand '-'>, <GenomicPosition object '4':154836045, strand '+'>, <GenomicPosition object '2':136717343, strand '+'>, <GenomicPosition object '15':73649839, strand '-'>, <GenomicPosition object '7':148651406, strand '-'>, <GenomicPosition object '3':29912426, strand '-'>, <GenomicPosition object '15':91602775, strand '+'>, <GenomicPosition object '7':133709879, strand '-'>, <GenomicPosition object '4':101516627, strand '-'>, <GenomicPosition object '7':16456750, strand '-'>, <GenomicPosition object '3':90262860, strand '-'>,


                however, i am getting plot of wincvg showing all values are not zero.
                Attached Files
                Last edited by sumanex; 02-02-2013, 11:24 PM.

                Comment


                • Hi Simon,

                  It looks like the HTSeq_example_data.tgz link is broken...can you please fix it?

                  Thanks!
                  Tommy

                  Comment


                  • Originally posted by crazyhottommy View Post
                    It looks like the HTSeq_example_data.tgz link is broken...can you please fix it?
                    Done.

                    Thanks for reporting it.

                    Comment


                    • Hi Simon,

                      I have a question regarding my data and numbers of counts resulting from htseq-count.

                      I have ~45m SE 80bp Illumina reads aligning to reference. I then remove duplicates (using Picard MarkDuplicates) ending up with ~25m reads. When I run htseq-count on the original aligned SAM and the duplicate-removed SAM I get ~23m, 8.8m counts.

                      The proportion of 'retained' reads (following duplicate removal) is ~0.53, whereas the proportion of 'retained' counts is ~0.39. I see this as quite a large discrepancy. Another way of looking at this is sum of counts over total reads in original and duplicate-removed. This gives the same proportions, ~0.5 in original SAM, ~0.37 in duplicate-removed.

                      I cannot see why the duplicate-removed SAM should give less counts than the original: if anything I would think it would give more because the duplicates would only give a single count (for highest quality if I understand htseq-count correctly). Is there an obvious answer or have you seen this before? I am using the 'union' mode with same reference as I aligned to, non-stranded data, all other parameters are default. I looked at the SAM flag and counted 0, 16 which give the same ratios for aligned vs. duplicate-removed reads. I would be interested to know your thoughts.

                      Thanks,

                      Bruce.

                      Comment


                      • new cigar operator

                        --Hi,

                        i think htseq-count doesn't understand the new CIGAR operators X and =
                        I have this error when reading a sam file:

                        Error occured when reading first line of sam file.
                        Error: ("Unknown CIGAR code 'X' encountered.", 'line 29 of file AB4_AR013.sam'
                        [Exception type: ValueError, raised in _HTSeq.pyx:1163]

                        Laurent --

                        Comment


                        • Hi,
                          I'm having the same issue with one of my SAM files when using htseq-count.
                          For some reason, one of my SAM files raises the error:

                          Error occured in line 38586430 of file RNA7_sorted.sam.
                          Error: ("invalid literal for int() with base 10: ''", 'line 38586430 of file RNA7_sorted.sam')
                          [Exception type: ValueError, raised in _HTSeq.pyx:1166]

                          I am using HTSeq-0.5.3p9

                          Can you help please, thanks, alig

                          Comment


                          • library type and stranded parameter

                            Hello,

                            I'm trying to figure out the right "stranded" parameter to use for my RNA-seq data which was aligned using TopHat with the "--library-type fr-firststrand" parameter. I'm using paired-end reads.

                            From what I can see, the results of running "stranded=no" is similar to "stranded=reverse" which gives me about ~50% of the total fragments, the majority have no feature. But if I ran using "stranded=yes", I only get ~2% of total fragments as having a feature.

                            I'm thinking that the "stranded=reverse" is the way to go if I want to measure sense expression, since for the fr-firststrand protocol, the right most strand is sequenced first which is opposite to the coding strand. Is this correct?

                            Thanks,
                            Patrick

                            Comment


                            • Hello,

                              I've used Tophat 2.0.9 & then HTseq version 0.5.4p3 & just with 3 of my 28 SAM files I get this error.

                              Error occured in line 63841485 of file RNA8_sorted.sam.
                              Error: ("'seq' and 'qualstr' do not have the same length.", 'line 63841485 of file RNA8_sorted.sam')
                              [Exception type: ValueError, raised in _HTSeq.pyx:772]

                              Can anyone please help as it's holding up my analysis.

                              Thank you
                              alig

                              Comment


                              • I'm thinking that the "stranded=reverse" is the way to go if I want to measure sense expression, since for the fr-firststrand protocol, the right most strand is sequenced first which is opposite to the coding strand. Is this correct?
                                Yes. I've posted on this here:http://seqanswers.com/forums/showpos...8&postcount=50

                                Comment

                                Latest Articles

                                Collapse

                                • seqadmin
                                  Latest Developments in Precision Medicine
                                  by seqadmin



                                  Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.

                                  Somatic Genomics
                                  “We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...
                                  05-24-2024, 01:16 PM
                                • seqadmin
                                  Recent Advances in Sequencing Analysis Tools
                                  by seqadmin


                                  The sequencing world is rapidly changing due to declining costs, enhanced accuracies, and the advent of newer, cutting-edge instruments. Equally important to these developments are improvements in sequencing analysis, a process that converts vast amounts of raw data into a comprehensible and meaningful form. This complex task requires expertise and the right analysis tools. In this article, we highlight the progress and innovation in sequencing analysis by reviewing several of the...
                                  05-06-2024, 07:48 AM

                                ad_right_rmr

                                Collapse

                                News

                                Collapse

                                Topics Statistics Last Post
                                Started by seqadmin, 05-24-2024, 07:15 AM
                                0 responses
                                16 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 05-23-2024, 10:28 AM
                                0 responses
                                18 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 05-23-2024, 07:35 AM
                                0 responses
                                22 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 05-22-2024, 02:06 PM
                                0 responses
                                11 views
                                0 likes
                                Last Post seqadmin  
                                Working...
                                X