SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
htseq-count paolo.kunder Bioinformatics 10 10-22-2014 04:45 AM
Issue with htseq-count cpleis Bioinformatics 8 10-15-2012 09:31 AM
multiBamCov or htseq-count to count read per feature ? NicoBxl Bioinformatics 1 07-03-2012 02:05 AM
htseq-count on UTR emilyjia2000 Bioinformatics 9 04-06-2012 12:13 PM
htseq-count gets more reads? deepsea Bioinformatics 3 03-29-2012 11:27 AM

Reply
 
Thread Tools
Old 10-10-2012, 02:23 AM   #1
EGrassi
Member
 
Location: Turin, Italy

Join Date: Oct 2010
Posts: 66
Default htseq-count eats 42G of memory

Hello,

I've been using DESeq for a while and right now htseq-count (version 0.5.3p9) it's giving me problems with a bam file (smaller than other that I've used with the same options and the same gtf), it runs and keeps occupying more memory, fills it and then start swapping, this is the relevant top output:
Code:
2958 data      20   0 46.5g  46g  676 D    6 97.5 841:10.30 htseq-count
Any idea? Thanks!
EGrassi is offline   Reply With Quote
Old 10-10-2012, 04:41 AM   #2
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

How have you called htseq-count? Please post the precise command line, and any output generated.
Simon Anders is offline   Reply With Quote
Old 10-10-2012, 05:03 AM   #3
EGrassi
Member
 
Location: Turin, Italy

Join Date: Oct 2010
Posts: 66
Default

Thank you for you answer:

Code:
samtools sort -n /rogue/bioinfotree/task/RNAseq/dataset/0.1/GSE27003/alignment_tophat/SRR097789.merged.bam SRR097789_sorted  
samtools view SRR097789_sorted.bam | htseq-count -m intersection-nonempty --stranded=no - annotation.gtf > nonoverlap_nonempty_counts_SRR097789_htseq 2> nonoverlap_nonempty_counts_SRR097789_htseq.err
No output reported, just the processing of the gtf infos.
I've manually inspected the bam file and it seems to have a lot of missing read ids (? I got the data from sra and just dumped it to fastq files, and they seem ok to me. Then standard tophat pipeline...), so now I'm trying to filter them out to see if that's the cause of the problem.
Samtools flagstat on the bam file does not report anything wrong with the bam file itself, apart from a very log percentage of properly paired reads.
EGrassi is offline   Reply With Quote
Old 10-10-2012, 05:09 AM   #4
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Htseq-count starts with reading in the annotation file and reports when its finished with that before looking at the reads. So, if it prints nothing, your problem is with the GTF file. How does it look like?
Simon Anders is offline   Reply With Quote
Old 10-10-2012, 05:13 AM   #5
EGrassi
Member
 
Location: Turin, Italy

Join Date: Oct 2010
Posts: 66
Default

Sorry, I did not point out htat properly maybe: "No output reported, just the processing of the gtf infos."
It finished the processing of the gtf, if you want the number I have to rerun it...but I swear that it printed them and that the same gtf worked with other bam files.

Last edited by EGrassi; 10-10-2012 at 05:31 AM.
EGrassi is offline   Reply With Quote
Old 10-10-2012, 05:36 AM   #6
EGrassi
Member
 
Location: Turin, Italy

Join Date: Oct 2010
Posts: 66
Default

Ok, I tried removing from the .sam files lines with an empty read name (I hope to understand why they're there) and it finished without any problem, the counts seems sensible to me (not from a biological point of view, ok, but they are numbers and sometimes they are different from 0!). If you want a sample of the .bam/.sam files with those strange lines I can give them to you (although I do not know if they are out of the bam format standard and samtools is just being nice to them or if they are ok and this could be a small bug of htseq_count/the library that it uses to scan bam files).
EGrassi is offline   Reply With Quote
Old 10-11-2012, 01:01 AM   #7
EGrassi
Member
 
Location: Turin, Italy

Join Date: Oct 2010
Posts: 66
Default ...because of tophat?

Ok, as long as fastq seems fine to me (they have an id foreach read, for example) I'm starting to think that the strange bam/sam obtained are tophat's responsibility.

I have "normal" lines like these ones:
Code:
SRR097789.29777200      89      chr19   55899346        255     50M     *       0       0       ATGCTCGCGCCNCGNTCAGCAGCATCAGACACATGATCCGCAAGAACAAG      AHFH@44455-!5-!:A:A:DHHHHHHHHHH=HEHDDDHDHHHHFFEGFG      AS:i:-2 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:11A2C35    YT:Z:UU XS:A:+  NH:i:1
And others without the QNAME field:

Code:
       329     chr1    10005   0       50M     *       0       0       CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC
      *       AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU NH:i:20 CC:Z:chr5       CP:i:10131      HI:i:0
Sequences gotten from lines like this one are found in the fastq files, associated with IDs.

Does anyone have an idea? I'm using tophat 2.0 with this command line:
Code:
tophat -p 7 --no-novel-juncs -G /rogue/bioinfotree/prj/ewing-rnaseq/local/share/data/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf --transcriptome-index=transcriptome -o SRR097789.th ../alignment/expanded_genome2 <(zcat /rogue/bioinfotree/task/RNAseq/dataset/0.1/GSE27003/reads//SRR097789_1.fastq.gz) <(zcat /rogue/bioinfotree/task/RNAseq/dataset/0.1/GSE27003/reads//SRR097789_2.fastq.gz)
The --no-novel-juncs option had to be added otherwise it will just freeze and stop after a while (see this other thread: http://seqanswers.com/forums/showthread.php?t=23887)
EGrassi is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:23 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO