What kind of tool would be appropriate do detect "broad peaks" in small RNA-seq sequencing data?

MACS2 appears to be developed for ChIP-seq data, but I see that there is a --nomodel option. Would that make this program usable in my case?

I suppose I should also set the --shift and --extsize options, but I'm not sure how:

  --nomodel             Whether or not to build the shifting model. If True,
                        MACS will not build model. by default it means
                        shifting size = 100, try to set extsize to change it.
                        DEFAULT: False
  --shift SHIFT         (NOT the legacy --shiftsize option!) The arbitrary
                        shift in bp. Use discretion while setting it other
                        than default value. When NOMODEL is set, MACS will use
                        this value to move cutting ends (5') towards 5'->3'
                        direction then apply EXTSIZE to extend them to
                        fragments. When this value is negative, ends will be
                        moved toward 3'->5' direction. Recommended to keep it
                        as default 0 for ChIP-Seq datasets, or -1 * half of
                        EXTSIZE together with EXTSIZE option for detecting
                        enriched cutting loci such as certain DNAseI-Seq
                        datasets. Note, you can't set values other than 0 if
                        format is BAMPE for paired-end data. DEFAULT: 0.
  --extsize EXTSIZE     The arbitrary extension size in bp. When nomodel is
                        true, MACS will use this value as fragment size to
                        extend each read towards 3' end, then pile them up.
                        It's exactly twice the number of obsolete SHIFTSIZE.
                        In previous language, each read is moved 5'->3'
                        direction to middle of fragment by 1/2 d, then
                        extended to both direction with 1/2 d. This is
                        equivalent to say each read is extended towards 5'->3'
                        into a d size fragment. DEFAULT: 200. EXTSIZE and
                        SHIFT can be combined when necessary. Check SHIFT

In my case the tags are the same things as the fragments, so this "extension" thing is probably not appropriate.

Basically, I would like to call peaks based on the plain real coverage of my small reads. Are there other tools that would be more appropriate?

  • $\begingroup$ Do you actually need peaks or is your end goal instead to find novel sRNAs? I would think stringTie followed by a considerable amount of filtering would be better. $\endgroup$
    – Devon Ryan
    Jun 8 '19 at 8:20
  • $\begingroup$ I mapped sRNA-seq reads on the genome, and I want to automatically delineate genomic regions that produce them. By "broad peaks", I mean such regions. Not sure what you mean by "novel sRNAs". $\endgroup$
    – bli
    Jun 10 '19 at 7:50
  • $\begingroup$ sRNAs shouldn’t appear as broad peaks unless they’re coming from repeat clusters. See the answer from Sebastian Müller for good tool options. $\endgroup$
    – Devon Ryan
    Jun 10 '19 at 8:34

The sure why any software should be confined to only broad peaks, but there is two packages that should do what I thing you have in mind. Just to clarify, for smallRNAs peak detection is sometimes also referred to as loci or cluster detection:

  • ShortStack which is developed for alignment, annotation, and quantification of small RNAs. The de-novo cluster finding option as copied from the manual might do the trick:
    Unless options --locus or --locifile are used (see below), ShortStack
    will de-novo identify clusters of small RNA accumulation genome-wide.
    Cluster definition is simple: First, all regions containing at least one
    primary alignment are found where the maximum distance between the ends
    of the alignments is <= option --pad (default: 75). Second, if the
    number of alignments in the cluster is >= option --mincov (default:
    0.5rpm), the cluster is kept. The mincov threshold can also be specified
    in terms of reads per million by using a value such as 3.2rpm (which
    specifies the threshold to be 3.2 reads per million). Using a rpm
    threshold allows the sensitivity of cluster discovery to be normalized
    between libraries of different sizes. Alternatively, reads per million
    mapped (rpmm) can be specified: A --mincov of 1.2rpmm indicates 1.2
    reads per million mapped is the threshold. rpm is a fraction of total
    library size, while rpmm is a fraction of only the aligned & placed
    fraction of the library.
  • segmentSeq is an R-package useing a Bayesian approach for small RNA locus detection. It's however really slow and not developed any further by the main author.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.