htseq-count not working with miRbase gff?
I have tried to count mapped reads on the human miRNA regions from the miRBAse human GFF
the output was just:
__no_feature 70216353 __ambiguous 0 __too_low_aQual 0 __not_aligned 1383340 __alignment_not_unique 0
So there should be something I do wrong?
I have run it - htseq-count SAMfile Gfffile (SAMfile is from bowtie with option -S)
The GFF is from MirBase:
chr1 . miRNA_primary_transcript 17369 17436 . - . ID=MI0022705;Alias=MI0022705;Name=hsa-mir-6859-1 chr1 . miRNA 17409 17431 . - . ID=MIMAT0027618;Alias=MIMAT0027618;Name=hsa-miR-6859-5p;Derives_from=MI0022705 chr1 . miRNA 17369 17391 . - . ID=MIMAT0027619;Alias=MIMAT0027619;Name=hsa-miR-6859-3p;Derives_from=MI0022705 chr1 . miRNA_primary_transcript 30366 30503 . + . ID=MI0006363;Alias=MI0006363;Name=hsa-mir-1302-2 chr1 . miRNA 30438 30458 . + . ID=MIMAT0005890;Alias=MIMAT0005890;Name=hsa-miR-1302;Derives_from=MI0006363 chr1 . miRNA_primary_transcript 187891 187958 . - . ID=MI0026420;Alias=MI0026420;Name=hsa-mir-6859-2 chr1 . miRNA 187931 187953 . - . ID=MIMAT0027618_1;Alias=MIMAT0027618;Name=hsa-miR-6859-5p;Derives_from=MI0026420 chr1 . miRNA 187891 187913 . - . ID=MIMAT0027619_1;Alias=MIMAT0027619;Name=hsa-miR-6859-3p;Derives_from=MI0026420 chr1 . miRNA_primary_transcript 632325 632413 . - . ID=MI0022558;Alias=MI0022558;Name=hsa-mir-6723 chr1 . miRNA 632382 632403 . - . ID=MIMAT0025855;Alias=MIMAT0025855;Name=hsa-miR-6723-5p;Derives_from=MI0022558
• 5,652 views
•
link
updated
by
Ram
•
written
by
manekineko
•
written
by
manekineko
1 answer
Did you use the -t flag to change the feature type to miRNA_primary_transcript or miRNA? As far as HTSeq would be concerned (without using this flag) you have given it no features to look for.
The default looks for exon in the third column
You should also use the -I flag to use 'ID' as the id attribute (default looks in the last column for 'gene_id')
Ensure you read the manuals before using these tools.
Log in to answer this question.
More posts like this
-
featurecounts not working on mirbase annotation file
written by demoraesdiogo2017 •
-
Remove multiple gene id from bed file
written by anii •
-
miRNA-seq count by "ID" or "Alias" using mirBase mmu.gff3
written by Rafael Soler •
-
Htseq is giving me 0 counts using the GFF3 of miRBase
written by Rafael Soler •
-
mirBase data to write a script to understand which is the dominant arm of a mirna duplex
written by bio_elle •
-
Counting reads with feature counts.
written by aranyak111 •
-
Micro analysis pipeline
written by aranyak111 •
-
convert hsa miRBase gff3 into gtf
written by anna.coenenstass •
-
How to find nearest CpG islands around miRNA genes using granges?
written by Adam •
-
Find Mirna Sequences
written by esmaeel •
Oh I forgot to use
-t, is it also necessary to use-i ID?And should the SAM file be sorted?
I will read the manual.
Sorry, I just edited my answer to say this as you submitted this.
Yes you should use
-i ID.Your Sam file is likely already sorted (many aligners do this for you) or HTSeq would have thrown an error.
Thanks!
since my SAM is quite big (15G) I wait quite a lot, befor I see results. Is there a way to see if it is counting properly with this flags in the middle of the run?
You can use samtools to subsection a random percentage of your input file, for example 5% would use
Sample Sam File
Is there a way to count reads but with option if they are longer than XXnt?
Not as far as I am aware.
You could control this during the preprocessing or alignment stages.
It is also possible to extract only reads of a specified length from a sam file Bam : Extract Only Alignment With A Specific Length
Also, please mark my answer if accepted if it has solved your original issue with read counting.
I found the solution - BBmap have tool to reformat SAM
Hi manekineko,
What program exactly did you download to reformat your sam file? Also, what do the options minlength and maxlength do?
Thanks!