htseq-general — Questions and discussions about HTSeq

[Htseq-general] Re ferred

[<da...@gm...>]

<style>
	.nisema {
   visibility: hidden;
}
</style>
	<div id="marofero" >
       <span class="nisema">be giant eagles with cruel steelcolored beaks and large brilliantly orange eyes The talons on 
</span>
</div>
<p></p><i>i've sent you a Offline message</i><p></p><a href="http://104.151.219.128/?dipu=razo&meyi=ZGF2bWxhd0BnbWFpbC5jb20=&id=aHRzZXEtZ2VuZXJhbEBsaXN0cy5zb3VyY2Vmb3JnZS5uZXQ=&lohizu=aHRzZXEtZ2VuZXJhbA==" >VIEW YOUR NOTIFICATION HERE</a>

[Htseq-general] Does --stranded=no in htseq-count affect resulting counts when using paired end reads?

[Komal R. <ko...@ma...>]

Hi everyone,

I have 64 paired-end RNAseq samples (human). I used cuffmerge to obtain a merged.gtf file but there are 3700 entries in my merged.gtf file that contain '.' in the strand field instead of '+' or '-'. I checked these entries and each of these have a class-code "u" which stands for "unknown". 

I used the merged.gtf file as input in my htseq-count command to get the counts (I eventually plan to do a differential gene expression analysis). However, I got the following error:

    Error occured when processing GFF file (line 360833 of file ./merged_asm/merged.gtf):
    Feature XLOC_003190 at chr1:[1285003,1285358)/. does not have strand information but you are running htseq-count in stranded mode. Use '--stranded=no'.
    [Exception type: SystemExit, raised in count.py:59]

I know the error is due to missing strand information. The only options that I can think of are:

    1) Use --stranded=no
    2) Remove all the entries in the merged.gtf corresponding to missing strand information

Does anybody have a suggestion? Do you think removing entries with missing information is better than using --stranded=no? I am afraid that using --stranded=no will affect the resulting counts (may increase the number of ambiguous counts when there is an overlap between exons/genes).

Thanks,

Komal S Rathi
Bioinformatics Application Developer
Perelman School of Medicine 
University of Pennsylvania 

[Htseq-general] HTseq count for miRNA

[Aggarwal, P. <pag...@mc...>]

Hi,

We did some miRNA sequencing (3 WT vs. 3 condition) and I would like to carry out differential expression analysis using DESeq (DESeq2). The alignment was done using the RNA-Star aligner (part of the ENCODE project). I know that HTSeq-count discards multi-mapped reads, however, for miRNA counts, discarding these reads would not capture the entire picture. A lot of the mature miRNA sequences are identical (especially when the miRNAs belong to the same family) and thus sequencing reads would be mapped to all the corresponding locations. Based on this, if we ignore these reads we are going to be left with zero to very few read counts for differential expression analysis.

HTseq count uses the NH column in the SAM output to look for multi mapped reads and RNA-STAR outputs this column (not all aligners do). So is there a way to turn the multi-mapped filtering OFF in HTSeq or should I use a custom counting script (instead of HTSeq) to generate my count tables. Also, do you think this would be a reasonable way to then use DESeq (we are currently focusing on validating an experimental idea/design and some of the interesting DE miRNAs would then be checked for using qPCR).

Any comments/suggestions on this issue will be greatly appreciated.

Thank you,
Praful

[Htseq-general] Bug fix fork of HTSeq

[David L. <da...@gm...>]

Hi guys, I really like HTSeq!

I found a few bugs, and sent them to the maintainer a few weeks ago
but have heard no response (holidays? Changed labs?)

I've put up a fork of the SVN code base here, until the patches are
put into the main branch:

https://bitbucket.org/sacgf/htseq


New Features:

htseq-count works with bam files
Wrote htseq-count-combine.py to merge multiple htseq-count outputs to
one big CSV file

Bug Fixes:

GFF Fix by David Eccles (gringer) on this mailing list
Fixed extend_to_include and GenomicInterval_from_directional attribute errors
Single ended bam files had paired_end incorrectly set
Sam/Bam mate start chrom incorrect default

[Htseq-general] GFF file read bug -- HTSeq crashes on GFF3 files with FASTA sections

[David E. (gringer) <bio...@gr...>]

[sending bug report here, because I couldn't find the usual reporting place on the sourceforge 
project page]

I'm using the HTSeq module for GFF file reading to read through the first stop codon (preliminary 
python script is attached), and came across an error right at the end of my file:

$ ./readthrough.py > ./protein_readthrough_S288C_R64-1-1_20110203.fa

Features read: 16405Traceback (most recent call last):
    File "./readthrough.py", line 41, in <module>
      for feature in gtfFile:
    File "/usr/local/lib/python2.7/dist-packages/HTSeq/__init__.py", line 214, in __iter__
      strand, frame, attributeStr ) = line.split( "\t", 8 )
ValueError: need more than 1 value to unpack


Looking at the GFF file, I realised/remembered that there was an extra section in the file:

$ tail -n +16423 
/home/bioinf/illuminadata/Genome/scer/S288C_reference_genome_R64-1-1_20110203/saccharomyces_cerevisiae_R64-1-1_20110208.gff 
| head
2-micron	SGD	CDS	3271	3816	.	+	0 
Parent=R0030W;Name=R0030W;gene=RAF1;Alias=RAF1;Ontology_term=GO:0003674,GO:0005634,GO:0006276;Note=Anti-repressor%20that%20increases%202%20micron%20plasmid%20copy%20number%20by%20relieving%20repression%20of%20the%20FLP1%20site-specific%20recombinase%20caused%20by%20the%20Rep1-Rep2p%20trascription%20regulator%3B%20also%20itself%20repressed%20by%20the%20Rep1p-Rep2p%20complex;dbxref=SGD:S000029674;orf_classification=Verified
2-micron	SGD	gene	5308	6198	.	-	. 
ID=R0040C;Name=R0040C;gene=REP2;Alias=REP2;Ontology_term=GO:0003674,GO:0005634,GO:0030543;Note=Master%20regulator%20that%20acts%20in%20concert%20with%20Rep1p%20to%20regulate%20transcript%20levels%20of%20the%20FLP1%20gene%20that%20promotes%20plasmid%20copy%20amplification%3B%20also%20autoregulates%20levels%20of%20its%20own%20transcript;dbxref=SGD:S000029676;orf_classification=Verified
2-micron	SGD	CDS	5308	6198	.	-	0 
Parent=R0040C;Name=R0040C;gene=REP2;Alias=REP2;Ontology_term=GO:0003674,GO:0005634,GO:0030543;Note=Master%20regulator%20that%20acts%20in%20concert%20with%20Rep1p%20to%20regulate%20transcript%20levels%20of%20the%20FLP1%20gene%20that%20promotes%20plasmid%20copy%20amplification%3B%20also%20autoregulates%20levels%20of%20its%20own%20transcript;dbxref=SGD:S000029676;orf_classification=Verified
###
##FASTA
  >chrI
CCACACCACACCCACACACCCACACACCACACCACACACCACACCACACCCACACACACACATCCTAACACTACCCTAAC
ACAGCCCTAATCTAACCCTGGCCAACCTGTCTCTCAACTTACCCTCCATTACCCTGCCTCCACTCGTTACCCTGTCCCAT
TCAACCATACCACTCCGAACCACCATCCATCCCTCTACTTACTACCACTCACCCACCGTTACCCTCCAATTACCCATATC
CAACCCACTGCCACTTACCCTACCATTACCCTACCATCCACCATGACCTACTCACCATACTGTTCTTCTACCCACCATAT

HTSeq seemed to be getting caught up by the FASTA section (or more accurately, the lines following 
the FASTA section).

I've attached a patch to fix this, which just identifies the three-hash directive in GFF3 files that 
indicates records are all done and breaks out of the iteration loop after that is seen. There's 
information about other directives here:

http://www.sequenceontology.org/gff3.shtml

Hope this helps,

   - David Eccles (gringer)