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.