SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
HTSeq dealing with "*" qualities kamsen Bioinformatics 39 10-15-2018 05:21 AM
"allele balance ratio" and "quality by depth" in VCF files efoss Bioinformatics 2 10-25-2011 11:13 AM
The position file formats ".clocs" and "_pos.txt"? Ist there any difference? elgor Illumina/Solexa 0 06-27-2011 07:55 AM
"Systems biology and administration" & "Genome generation: no engineering allowed" seb567 Bioinformatics 0 05-25-2010 12:19 PM
SEQanswers second "publication": "How to map billions of short reads onto genomes" ECO Literature Watch 0 06-29-2009 11:49 PM

Reply
 
Thread Tools
Old 07-11-2012, 07:39 AM   #1
dmsoanes
Junior Member
 
Location: Exeter

Join Date: Jul 2012
Posts: 2
Default HTseq count - all reads are "alignment_not_unique"

I have been having problems using HTseq-count.
I’ll go through my protocol
1: Take tophat output (accepted_hits.bam), sort them by read and convert to bam file.

samtools sort -n mst12_t14_2_accepted_hits.bam mst12_t14_2_accepted_hits_read_sorted

samtools view mst12_t14_2_accepted_hits_read_sorted.bam > mst12_t14_2_accepted_hits_read_sorted.sam


2: Use HTseq-count to count reads for each gene using GFF file used for tophat alignment

htseq-count -m union -s no -t exon -i gene_id mst12_t14_2_accepted_hits_read_sorted.sam magnaporthe_oryzae_70-15_8_transcripts.gff > test.txt

There were no error messages but when I looked at the output, the count for every gene is zero and all the reads are in “alignment_not_unique” category. Reads in this category supposedly have NH > 1 in SAM file, but the reads I can see all have NH=1.
Below is a sample of the SAM file and the GFF.

D3P26HQ1:129:D0TA0ACXX:2:1101:1108:73648 73 supercont8.4 4401713 50 35M * 0 0
CCTTGATGTGGACGCGCCCNGCGACCTTGCTGGCG ?=@BDDDDCF>D>EGG6FF#1:?DHDB@DH9D>;F AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:19G15 YT:Z:UU NH:i:1
D3P26HQ1:129:D0TA0ACXX:2:1101:1109:25458 89 supercont8.4 1740219 50 35M * 0 0
AGCGAATGACGAGCGNCGCTATCATGATAACGTCA ?091D?)1:?<<333#@2EDE<DC<>C;?2=41?: AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:15A19 YT:Z:UU NH:i:1
D3P26HQ1:129:D0TA0ACXX:2:1101:1112:73199 99 supercont8.1 3567120 50 36M = 3567285 210
TGCACAGAGAGATGGTTGCNTATATATATGCAAAAA ?<;ADDDD+<C;DE@2A@@#33A1CBEEC999CDDD AS:i:-1 XN:i:0 XM:i:1 XO:i:0
XG:i:0 NM:i:1 MD:Z:19A16 YT:Z:UU NH:i:1
D3P26HQ1:129:D0TA0ACXX:2:1101:1112:73199 147 supercont8.1 3567285 50 45M = 3567120 -210 CACAGACAAAGCGGCTCGTCTTGGCGCGGGCACGCATTGGCTCGG @CBHE?C?EGA@A;-(.9DDED0BD@@;<)B2DA<3HD80+BB?@ AS:i:-5 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:15G6T22 YT:Z:UU NH:i:1
D3P26HQ1:129:D0TA0ACXX:2:1101:1113:7390 89 supercont8.6 1611568 50 34M * 0 0 AGTACGTCTCGGCTNTCCAGAAGGCCAAGATCAC ?:)1CC?)1C:+33#3BC:EECA++A?BBB+??? AS:i:-3 XM:i:2 XO:i:0 XG:i:0 MD:Z:2G11T19 NM:i:2 NH:i:1

supercont8.5 MG8_CALLGENES_FINAL_1 start_codon 4470284 4470286 . - 0 gene_id "MGG_00005"; transcript_id "MGG_00005T0";
supercont8.5 MG8_CALLGENES_FINAL_1 stop_codon 4469688 4469690 . - 0 gene_id "MGG_00005"; transcript_id "MGG_00005T0";
supercont8.5 MG8_CALLGENES_FINAL_1 exon 4470173 4470286 . - . gene_id "MGG_00005"; transcript_id "MGG_00005T0";
supercont8.5 MG8_CALLGENES_FINAL_1 CDS 4470173 4470286 . - 0 gene_id "MGG_00005"; transcript_id "MGG_00005T0";
supercont8.5 MG8_CALLGENES_FINAL_1 exon 4470034 4470120 . - . gene_id "MGG_00005"; transcript_id "MGG_00005T0";
supercont8.5 MG8_CALLGENES_FINAL_1 CDS 4470034 4470120 . - 0 gene_id "MGG_00005"; transcript_id "MGG_00005T0";
supercont8.5 MG8_CALLGENES_FINAL_1 exon 4469273 4469969 . - . gene_id "MGG_00005"; transcript_id "MGG_00005T0";
supercont8.5 MG8_CALLGENES_FINAL_1 CDS 4469691 4469969 . - 0 gene_id "MGG_00005"; transcript_id "MGG_00005T0";


Any ideas?

Cheers, Darren
dmsoanes is offline   Reply With Quote
Old 07-11-2012, 09:55 AM   #2
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

I'd guess you are using the version of two weeks ago, which contained a bug misreading the "NH" tag. Try the newest version.
Simon Anders is offline   Reply With Quote
Old 07-12-2012, 05:07 AM   #3
dmsoanes
Junior Member
 
Location: Exeter

Join Date: Jul 2012
Posts: 2
Default

Thanks for your swift reponse Simon, I have installed the latest version (HTSeq-0.5.3p9) and everything works fine.
dmsoanes 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 06:00 PM.


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