Question: Htseq Count
0
gravatar for sethugunja
5.1 years ago by
sethugunja60
Sweden
sethugunja60 wrote:

Hi guys,

I want to use htseq-edgeR to validate my differentially expression cuffdiff results.

For this, first htseq-count have to be used to count the number of reads for the features i have but htseq-count says:

If you have paired-end data, you have to sort the SAM file by read name first.

I converted the tophat bam file into sam file and sorted by readname by using samtools, and given the sorted sam file as an input for htseq-count.... but it says in the slurm file: Warning: Read HISEQ:64:D25AKACXX:1:2211:9696:78749 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?) Warning: Read HISEQ:64:D25AKACXX:1:2309:5639:89052 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?) Warning: Skipping read 'HISEQ:64:D25AKACXX:1:2309:5639:89052', because chromosome 'GL000191.1', to which it has been aligned, did not appear in the GFF file. Warning: Skipping read 'HISEQ:64:D25AKACXX:1:1112:7444:68675', because chromosome 'GL000191.1', to which it has been aligned, did not appear in the GFF file.

I didnt understand why i got this warning even I have given the sorted file and the output file looks like this:

HISEQ:64:D25AKACXX:1:1101:10000:13676 83 chr1 179120449 50 100M = 179120415 -134 AGGTATAGTTTAAGTTCCTAGATTCTCCTTGAACCTGACCTCTGGTCATTTCAGATTAAGTAATCTCTTCCTTTGCCTAAGCTCAAAGCACCCATTGCCT DDDDDDDDDDDDDECCEEEFFFFFHHHHHJIE@JJIHFJJIJJJJJJJJJJJJJJJJJJJJIJJJJJJJJJIJJJJJJJJJJJJJJJHHHHHFFFFFCCC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:100 YT:Z:UU XS:A:- NH:i:1 XF:Z:no_feature HISEQ:64:D25AKACXX:1:1101:10000:13676 163 chr1 179120415 50 100M = 179120449 134 CCCGTCTGTGTATTCAAATCCTATATTTCTTTCAAGGTATAGTTTAAGTTCCTAGATTCTCCTTGAACCTGACCTCTGGTCATTTCAGATTAAGTAATCT BCCFDDFFHHFHHJJJJJJJJJJJJJJJJJJJJJJJJCHIIJHIJJJJJJIIJJJJJJJJJJJJJJJJJJJJJJJJJJIHIJJJHHGHHHHFFFFFFFED AS:i:-5 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:3C96 YT:Z:UU XS:A:- NH:i:1 XF:Z:no_feature HISEQ:64:D25AKACXX:1:1101:10000:17369 83 chr11 73120439 50 100M = 73120346 -193 GCTGTGCTGCAGCCAATTCTATGACTGTGTTCTCTGCTCCGTTTTGAGTTCCTAGAGGAAAGGGATGGACACTGATGCAGATGCTGCCCTCACCTGCTGA DDDDDCDDDDDDDDDCCDCDEECEEDFFFFFHHHHIJJJJJJIIIJJJIJIJJJJJJJJJJJJJJIJIIIGJJJJJJJJJJJIJJIIHHHHHFFFFFCCC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:100 YT:Z:UU XS:A:+ NH:i:1 XF:Z:ENSG00000054965.6 HISEQ:64:D25AKACXX:1:1101:10000:17369 163 chr11 73120346 50 100M = 73120439 193 CGTCCTTCTTGGTGAAATCTTCCTCATCCCCATAGACAGATTTAGAGGCTCCTTCCTCTGGGTTCCCAAGCCCCCTGGCACAGGATCTCTTGTGCTGTGC ?@BDFFFFHHHHHJJJIJJJJJJJJJJJJJJJJJJJJJJIJJJJJIIJJIHIIJJJJGHEII8CHGIJJJHEHHFFCDEEEDDDDDDDDECDDCDDCCDC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:100 YT:Z:UU XS:A:+ NH:i:1 XF:Z:ENSG00000054965.6

Do you think I m going in a correct way?

Is their any proper way of doing the htseq count? Please suggest me .

Thanks, Sethu

• 9.2k views
ADD COMMENTlink modified 5.1 years ago by Devon Ryan86k • written 5.1 years ago by sethugunja60
3
gravatar for Devon Ryan
5.1 years ago by
Devon Ryan86k
Freiburg, Germany
Devon Ryan86k wrote:

The warnings are fine and you're not doing anything wrong. The "claims to have an aligned mate which could not be found" warnings are due to the fact that tophat won't put in both reads of a pair in the alignment file if only one of them maps (the better solution would be to put both in there). HTSeq handles this properly but issues a warning to alert you. You can ignore this warning in your case. The "because chromosome 'GL000191.1', to which it has been aligned, did not appear in the GFF file" warnings are due to there being contigs without annotated genes. That's also fine in your case, the warning is just to alert you in case that would indicate an error.

ADD COMMENTlink written 5.1 years ago by Devon Ryan86k

Thankyou very much for your clarification!

ADD REPLYlink written 5.1 years ago by sethugunja60

Hi, The output file which I got from htseq-count is a very large file of size around 10 gb. I didnt get any count number but looks like :

HISEQ:64:D25AKACXX:1:1101:10000:13676 83 chr1 179120449 50 100M = 179120415 -134 AGGTATAGTTTAAGTTCCTAGATTCTCCTTGAACCTGACCTCTGGTCATTTCAGATTAAGTAATCTCTTCCTTTGCCTAAGCTCAAAGCACCCATTGCCT DDDDDDDDDDDDDECCEEEFFFFFHHHHHJIE@JJIHFJJIJJJJJJJJJJJJJJJJJJJJIJJJJJJJJJIJJJJJJJJJJJJJJJHHHHHFFFFFCCC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:100 YT:Z:UU XS:A:- NH:i:1 XF:Z:no_feature

HISEQ:64:D25AKACXX:1:1101:10000:13676 163 chr1 179120415 50 100M = 179120449 134 CCCGTCTGTGTATTCAAATCCTATATTTCTTTCAAGGTATAGTTTAAGTTCCTAGATTCTCCTTGAACCTGACCTCTGGTCATTTCAGATTAAGTAATCT BCCFDDFFHHFHHJJJJJJJJJJJJJJJJJJJJJJJJCHIIJHIJJJJJJIIJJJJJJJJJJJJJJJJJJJJJJJJJJIHIJJJHHGHHHHFFFFFFFED AS:i:-5 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:3C96 YT:Z:UU XS:A:- NH:i:1 XF:Z:no_feature

HISEQ:64:D25AKACXX:1:1101:10000:17369 83 chr11 73120439 50 100M = 73120346 -193 GCTGTGCTGCAGCCAATTCTATGACTGTGTTCTCTGCTCCGTTTTGAGTTCCTAGAGGAAAGGGATGGACACTGATGCAGATGCTGCCCTCACCTGCTGA DDDDDCDDDDDDDDDCCDCDEECEEDFFFFFHHHHIJJJJJJIIIJJJIJIJJJJJJJJJJJJJJIJIIIGJJJJJJJJJJJIJJIIHHHHHFFFFFCCC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:100 YT:Z:UU XS:A:+ NH:i:1 XF:Z:ENSG00000054965.6

HISEQ:64:D25AKACXX:1:1101:10000:17369 163 chr11 73120346 50 100M = 73120439 193 CGTCCTTCTTGGTGAAATCTTCCTCATCCCCATAGACAGATTTAGAGGCTCCTTCCTCTGGGTTCCCAAGCCCCCTGGCACAGGATCTCTTGTGCTGTGC ?@BDFFFFHHHHHJJJIJJJJJJJJJJJJJJJJJJJJJJIJJJJJIIJJIHIIJJJJGHEII8CHGIJJJHEHHFFCDEEEDDDDDDDDECDDCDDCCDC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:100 YT:Z:UU XS:A:+ NH:i:1 XF:Z:ENSG00000054965.6

I m gng to use EdgeR to identify the differentially expressed genes which needs count number and geneIds. I heard that htseq count output file contains gene Ids and count numbers , but my output file doesnt have these variables.

Where did I go wrong?

Thanks Sethu

ADD REPLYlink modified 5.0 years ago • written 5.0 years ago by sethugunja60

That's not the output file you need. Don't use the -o option unless you need to debug something or have other special needs. Instead, use htseq-count [options] file.sam annotation.gtf > output.txt. The output file will then be small and what you actually need.

ADD REPLYlink written 5.0 years ago by Devon Ryan86k

yes, I did the same way as you suggested, but the text file itself contains the above data with 10 gb size

ADD REPLYlink written 5.0 years ago by sethugunja60

No, you didn't do what I suggested, or else the file snippet you posted would have never been produced. Post the exact command you used with nothing removed.

ADD REPLYlink written 5.0 years ago by Devon Ryan86k

Hi, following is the pipeline i have used:

samtools sort -n accepted_hits.bam accepted_hits_bam_sorted

samtools view -h accepted_hits_bam_sorted.bam >accepted_hits_sorted.sam

htseq-count -s no -i gene_id -o htseq_P-6_PA_results.txt accepted_hits_sorted.sam gencode.v17.annotation.gtf

ADD REPLYlink modified 5.0 years ago • written 5.0 years ago by sethugunja60

Yes, exactly, so you didn't do what I suggested ("Don't use the -o option..."). Instead you need htseq-count -s no -i gene_id accepted_hits_sorted.sam gencode.v17.annotation.gtf > htseq_P-6_PA_results.txt though you could also just skip a step in your pipeline and samtools view accepted_hits_bam_sorted.bam | htseq-count -s no -i gene_id - gencode.v17.annotation.gtf > htseq_P-6_PA_results.txt. Technically you could pipe in the sort step too, but that won't really save you much time.

ADD REPLYlink written 5.0 years ago by Devon Ryan86k

Thank you for your prompt reply.

ADD REPLYlink written 5.0 years ago by sethugunja60
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1305 users visited in the last hour