htseq-count : 0 counts per miRNA

[JPerales]

Hi everyone,

I have fastq files from Illumina small RNA-Seq. My first goal is to analyze differential expression of known miRNAs.

I trimmed the adapter 3' with cutadapt and mapped the reads with Bowtie. I'm running htseq-count with the hsa.gff3 file from miRBase, but I get 0 counts per miRNA.

I don't understand what I am doing wrong. I need some help.

Thanks a lot!

I used these commands:

bowtie --sam -n 0 -l 15 -e 99999 -k 100 --best --chunkmbs 128 Reference_genome/h_sapiens_37_asm.ebwt/h_sapiens_37_asm reads_trimmed.fasq > aligned.sam

htseq-count -m union -t miRNA -s yes -i Name aligned.sam miRBase/genomes/hsa.gff3

[Simon Anders]

Are you sure the coordinates match, i.e., does the GFF file from miRNA also use assembly 37? And do the chromosome names match? (i.e., either "chr3" in both files, or just "3" in both files, but no mix.)

[JPerales]

Hi Simon,

Thanks for your fast answer. I have little experience in bioinformatics approach.

I used the indexed reference genome from Bowtie official web-page (Assembly 37).
I think that you are in the correct way about chromosome names match.

I show you the head of my SAM file:

egrep -ve "^@" aligned.sam | head -n 10

Head:

HWUSI:3:1:7622:1103 1:N:0:GTTTCG	4	*	0	0	*	*	0	0	AAACCGTTACCATNACTGAGTT	HHHHGHHHHHEGG#G;@=@>HH	XM:i:0
HWUSI:3:1:4646:1103 1:N:0:GTTTCG	4	*	0	0	*	*	0	0	TTGGGATGTATTCNTACTGTCTGATGTGGAATTCTTGGGT	HBHHHGGDBG>BC#BABBCAGGBGGGGGGDGGG<@*??##	XM:i:0
HWUSI:3:1:11655:1104 1:N:0:GTTTCG	4	*	0	0	*	*	0	0	GCATTGGTGGTTCNGTGGTAGAATTCTCGC	HHHHHHHHHEDDE#EAA?CBHDGHHHHHHH	XM:i:0
HWUSI:3:1:3653:1104 1:N:0:GTTTCG	4	*	0	0	*	*	0	0	TGGTTTTCGGAACNGAGGC	HHDHHGDGG<:?:#=977>	XM:i:0
HWUSI:3:1:10014:1105 1:N:0:GTTTCG	4	*	0	0	*	*	0	0	ACGAGAACTTTGANGGCCGAAGTGGAGAAGGGTGGAATTC	IIIIIIIIIIFFF#EBEEDBIIHIIFHIIIII?IIDGEGG	XM:i:0
HWUSI:3:1:7574:1106	16	gi|224589820|ref|NC_000008.10|	141742723	255	20M	*	0	0	CTAGACTGTGAGCTCCTCGA	BGGGEGHGHHDDGGDEGGGG	XA:i:0	MD:Z:20	NM:i:0
HWUSI:3:1:7574:1106	16	gi|224589816|ref|NC_000004.11|	17448685	255	20M	*	0	0	CTAGACTGTGAGCTCCTCGA	BGGGEGHGHHDDGGDEGGGG	XA:i:0	MD:Z:20	NM:i:0
HWUSI:3:1:7574:1106	16	gi|224589805|ref|NC_000014.8|	100575773	255	20M	*	0	0	CTAGACTGTGAGCTCCTCGA	BGGGEGHGHHDDGGDEGGGG	XA:i:0	MD:Z:20	NM:i:0
HWUSI:3:1:7574:1106	16	gi|224589822|ref|NC_000023.10|	53408605	255	20M	*	0	0	CTAGACTGTGAGCTCCTCGA	BGGGEGHGHHDDGGDEGGGG	XA:i:0	MD:Z:20	NM:i:0
HWUSI:3:1:7574:1106	0	gi|224589800|ref|NC_000001.10|	99358453	255	20M	*	0	0	TCGAGGAGCTCACAGTCTAG	GGGGEDGGDDHHGHGEGGGB	XA:i:0	MD:Z:20	NM:i:0

The 3th field from my SAM file (RNAME ~ reference sequence name) doesn't seem like "chr3" for example (first field in my GFF file). i.e, RNAME = gi|224589820|ref|NC_000008.10|, this is chr8.

Can you give me some advice?

Maybe I should use other reference genome, or do a conversion in this field with "sed -e".

Thanks so much for your help.

[Simon Anders]

Maybe I should use other reference genome, or do a conversion in this field with "sed -e".

The reference genome should match (both your file refer to assembly 37), so, yes, you should simply fix the chromosome names with 'sed'.

[JPerales]

SOLVED!
Thanks you so much Simon.

I fixed RNAME field in the SAM file to his corresponding chromosome name and it works very well. It's my solution : fix it with "sed", "awk", or scripting.