Bacterial RNA-Seq - Problem Counting Reads with htseq

[umnklang]

I'm having problems utilizing htseq.

I have a sam file of paired end illumina reads. These have been sorted by the read name (1st column). Here are the first few lines of data:


HWI-ST261:3960D48ABXX:4:1101:10002:96209 147 gi|237809907|ref|NC_012692.1| 30581 255 79M = 30446 -214 AGGCCCCCAGCGTCATGCAGGCGCTCTCCGCCTACCGCGATGGCGCGCGCATCCATGCCGGCACGCGACCGGGCGCACC #####################################?C<562'93>>6BFHBA<@B=A7(:A@@@IIGHFDHHDDD XA:i:0 MD:Z:1T4G4T3G9G0C0T51 NM:i:7
HWI-ST261:3960D48ABXX:4:1101:10002:96209 99 gi|237809907|ref|NC_012692.1| 30446 255 79M = 30581 214 AGCAGCCGGCGCTTTACTGGCACTTCAGGAACAAGCGGGCGCTGCTCGACGCACTGGCCGAAGCCATGCTGGGGGAGAA FFHGHHHHGHHJJHIIIIIJJJIIIJJIFIIIGIJ@@EHF?8=@85:<0:<BD@59AADB################### XA:i:0 MD:Z:72C6 NM:i:1
HWI-ST261:3960D48ABXX:4:1101:10007:167767 163 gi|237809907|ref|NC_012692.1| 115957 255 79M = 116090 212 ACTGGTCCAGAACCTTGACCGAACGCAGCGGTGGTAACGGCGCAGTGGCGGTTTTCATGGCTTTTTTTTTTCTGTTTTT FFHHHHHJJJGIJJJJIGIIIJJJJJGJJJJFHGHGGHHHFDDDDDCDDDD00;B######################## XA:i:0 MD:Z:63G2A1G0A0C0T0G0T5 NM:i:8
HWI-ST261:3960D48ABXX:4:1101:10007:167767 83 gi|237809907|ref|NC_012692.1| 116090 255 79M = 115957 -212 ATGTTTGATGTTATGGAGCAGCAACGATGTTACGCAGCAGGGCAGTCGCCCTAAAACAAAGTTAGACATCATGAGGGTA DDDDDDDDDDDDDDDDDDDDDEDDFFFFFHHHJJJJJJJJJJJJJJJJJIJJJJJJJJJIHJJJHIGJJJJIHHHHHFF XA:i:0 MD:Z:79 NM:i:0
HWI-ST261:3960D48ABXX:4:1101:10013:104772 163 gi|237809907|ref|NC_012692.1| 119786 255 79M = 119940 233 TCGAAGACGCCCTGCATGCCACCCGTGCGGCGGTCGAGGAAGGCATCGTCCCGGGCGGCGGGGCCGCCCCGGTCCGTGC FFHGHFHGIIIIIGIADFHHJJJJIAGIIIHDD?B6-9:&8<@??9?@DA@D&077;@##################### XA:i:0 MD:Z:61C1T5T1A7 NM:i:4
HWI-ST261:3960D48ABXX:4:1101:10013:104772 83 gi|237809907|ref|NC_012692.1| 119940 255 79M = 119786 -233 GAAAGCCCCGCGCCGCGAGTTCGCCCCCAACGCCGGCGACGAGCCGAGCGTGGTGCTGAACCGCGTCGCCGAAGGCACC ########################################B<@@;>;BCCDDB?@@BEA<IHBEHHHGIGJIFHFADFF XA:i:0 MD:Z:1G9T0G6A3T1A53 NM:i:6

Next I have a GTF file. Here the first few lines:

##gff-version 3
#!gff-spec-version 1.20
#!processor NCBI annotwriter
##sequence-region NC_012692.1 1 166530
##species http://www.ncbi.nlm.nih.gov/Taxonomy...tax.cgi?id=562
NC_012692.1 RefSeq region 1 166530 . + . ID=id0;Dbxref=taxon:562;Is_circular=true;gbkey=Src;genome=plasmid;mol_type=genomic DNA;plasmid-name=pAR060302;strain=AR060302
NC_012692.1 RefSeq gene 323 862 . + . ID=gene0;Name=pAR060302_0001;Dbxref=GeneID:7872685;gbkey=Gene;locus_tag=pAR060302_0001
NC_012692.1 RefSeq CDS 323 862 . + 0 ID=cds0;Name=YP_002894347.1;Parent=gene0;Note=identified by Glimmer 2%3B sequence similarity to pSN254_0002%3B blastp_match %3D ZP_01176338.1;Dbxref=Genbank:YP_002894347.1,GeneID:7872685;gbkey=CDS;product=hypothetical protein;protein_id=YP_002894347.1;transl_table=11
NC_012692.1 RefSeq gene 867 1154 . + . ID=gene1;Name=pAR060302_0002;Dbxref=GeneID:7872573;gbkey=Gene;locus_tag=pAR060302_0002
NC_012692.1 RefSeq CDS 867 1154 . + 0 ID=cds1;Name=YP_002894348.1;Parent=gene1;Note=identified by Glimmer 2%3B sequence similarity to pSN254_0002%3B blastp_match %3D ZP_01176339.1;Dbxref=Genbank:YP_002894348.1,GeneID:7872573;gbkey=CDS;product=hypothetical protein;protein_id=YP_002894348.1;transl_table=11
NC_012692.1 RefSeq gene 1139 2239 . + . ID=gene2;Name=repA;Dbxref=GeneID:7872574;gbkey=Gene;gene=repA;locus_tag=pAR060302_0003
NC_012692.1 RefSeq CDS 1139 2239 . + 0 ID=cds2;Name=YP_002894349.1;Parent=gene2;Note=identified by Glimmer 2%3B sequence similarity to pSN254_0003%3B blastp_match %3D ZP_01176340.1;Dbxref=Genbank:YP_002894349.1,GeneID:7872574;gbkey=CDS;product=plasmid replication protein RepA;protein_id=YP_002894349.1;transl_table=11

What I want to do is count the reads that align to the CDS and return a file with the GeneID and the number of reads that align to it so that I can use these in subsequent edgeR analyses.

The command I run is

-m HTSeq.scripts.count -m intersection-nonempty -t CDS -s no -i Dbxref pAR_sorted.sam NC_012692.gff > test.txt

in my output I get no reads aligning to any CDS

no_feature 49397
ambiguous 0
too_low_aQual 0
not_aligned 0
alignment_not_unique 0


However I know that this cannot be true - I have viewed this data on IGV and there are many many reads aligned to the CDS

Any thoughts?

[kopi-o]

I have run HTSeq successfully on similar GFF files with the options

-i locus_tag -t Gene

I wonder if it is the Dbxref that is messing up things for you, but I am not sure. You could try to use "Gene" and "locus_tag" or "Name" instead.

[umnklang]

I've found the answer.

In my sam file the reference (3 column) the accession was given as gi|237809907|ref|NC_012692.1|.

But my gff file it was simply NC_012692.1. A simple find (NC_012692.1) replace (gi|237809907|ref|NC_012692.1|) in my gff file is all I needed.

[JUANPABLOROSAS5]

Thanks!!!! great detail