Htseq Count

Question

[sethugunja]

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

Answer 1

[Devon Ryan]

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.

[sethugunja]

Thankyou very much for your clarification!

[sethugunja]

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

[Devon Ryan]

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.

[sethugunja]

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

[Devon Ryan]

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.

[sethugunja]

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

[Devon Ryan]

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.

[sethugunja]

Thank you for your prompt reply.