Peak calling for microRNA

Question

[JulianC]

Hi everyone!

I'm working with microRNA, I have miRNA-Seq from different tissue and I want to perform DEA among them, but only considering a small genomic region. I aligned my reads against my reference, and I visualized them on IGV. The peaks I see do not correspond to any annotated miRNA for my reference genome, so I have not an annotation. In order to perform DEA, I want to create a Count Table containing the number of reads for each peak I see in my region of interest (I want to use it for DESeq2) . How can I do this peak-calling operation for mapped miRNAs? I have only the BAM and SAM files derived from the alignment (performed with Bowtie2) and the reference genome (FASTA). Thanks in advance!

[Devon Ryan]

Have you looked at the repeatMasker track for that area on the UCSC genome browser? It's not unheard of to get a bunch of other ncRNAs (fragmented rRNAs, etc.) in smallRNAseq datasets? Are these the appropriate length for miRNA in your species?

[JulianC]

In UCSC genome browser there's not the new assembly for my species, not yet. However from NCBI report I know that in my species there are 728 annotated miRNA (but I think I found different ones) with a mean lenght of 23 bp. My problem is essentially to get the read counting for the peaks I see on IGV (few peaks actually)

Answer 1

[Devon Ryan]

There are probably some miRNA specific tools for this, but worst case scenario you could use a peak caller like MACS2. If you use that, you'll need to use something like --keep-dup 1000000 to ensure it keep the big stacks of reads you're likely seeing intact. You might have to filter the results a bit to remove false negative, or possible prefilter your BAM files if you have a bunch of alignments of the wrong length (you might be able to use alignmentSieve from deepTools for this, though I don't think I've ever tested it with single-end reads), but in general I would think that that should work.