[Simon Anders]

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

[Jerry_Zhao]

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

[gringer]

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/anders/HTSeq/doc/tss.html#using-the-full-coverage

[sumanex]

matplot coverage

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

of cvg vector to get a decent plot.

thanks
suman

[sumanex]

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

[gringer]

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

[sumanex]

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.

[crazyhottommy]

Hi Simon,

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

Thanks!
Tommy

[Simon Anders]

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

Done.

Thanks for reporting it.

[bruce01]

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.

[mslider]

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 --

[alig]

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

[ppatrickt]

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

[alig]

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

[dvanic]

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