How to best detect the "peaks" in RNA-seq data that are not assigned to any gene?

Question

I encountered that many reads from single-cell RNA seq data were lost in the analysis because not assigned to any gene (genome: galgal6). I am trying to find an approach than could give me all the "peaks" in this data, that is to say regions where there is a high density of reads. Once I get there, I would extract "peaks" that are not assigned to any gene, and get the distance to the closest genes.

I am thinking of using first bedtools genomecov to get the "peaks", and then bedtools intersect to get all reads assigned to a gene. Finally, I extract reads that are not in the output of bedtools intersect to get a new file with all my reads from the peak, that are not assigned to any gene.

Does it make sense this way ? Can you see any other way of doing ?

Here are some tiny data:

> cat reads.bed
chr9    505479  505498
chr9    508014  508037
chr9    514603  514633
chr9    529519  529540
chr9    529519  529540
chr9    529519  529540

> cat tiny.galGal6.chrom.sizes
chr9    24153086

> bedtools genomecov -bg -split -i reads.bed -g tiny.galGal6.chrom.sizes
chr9    505479  505498  1
chr9    508014  508037  1
chr9    514603  514633  1
chr9    529519  529540  3

The last line would correspond to a peak. However, another difficulty is how should I define the thresholds ?

Any help would be more than welcome.

Answer 1

Basically what you're discovering is that there are unannotated expressed features, so your task isn't really finding peaks, but rather finding novel expressed transcripts. For that, you can use stringTie or similar programs. Ensure that the BAM file you give to stringTie have all of the cells together.

Answer 2

Many pipelines have been set up for this such as software developed for the FANTOM Project at RIKEN. The software for this is available from our website. These were set up for analysis of CAGE sequencing data for exactly this purpose of calling peaks of previously unannotated gene expression. You can also find data for many species (including chicken) from FANTOM5 where this has already been done.