Tool to retrieve regions with reads?

Question

[manekineko]

Hi,

Is there a tool for retrieving genomic regions occupied with reads as in the picture (reads are ~17-25nt)? I was thinking to use assembly on the reads to reconstruct the regions covered with reads but maybe they are not so dense/overlapped, so using velvet with kmer 15 I got some regions reconstructed but very very small part? So any idea will help, probably there is more smart method to do what I want?

[manekineko]

No, it is a plant and yes they can be TE regions expressing small RNAs, but I'm interested in them too. So want to collect such regions and have the sequence that was occupied by such cloud of reads.

Answer 1

[James Ashmore]

If I understand your question correctly, you can retrieve the regions using bedtools. First make a coverage track of reads along the genome (only recording regions with coverage > 0). Then merge any regions which are book-ended (You can change this limit, to merge regions which are within N bp of each other). Then get the FASTA sequence of the merged regions.

bedtools genomcov -ibam sample.bam -bg > sample.bedgraph
sort -k1,1 -k2,2n sample.bedgraph > sample.sorted.bedgraph
bedtools merge -i sample.sorted.bedgraph > sample.sorted.merged.bedgraph
bedtools getfasta -fi genome.fasta -bed sample.sorted.merged.bedgraph -fo regions.fasta

[manekineko]

Thanks , seems what I need probably, will try it!

[manekineko]

If I want to make GFF of the regions and not get their sequence how should I proceed?

Particulary I want to see in sample.sorted.merged.bedgraph where are intersections with another GFF?

Can I do it just by merging all my mapped reads with close to each other within Nbp? so skipping genomecov?

[James Ashmore]

A simple approach would be to change the GFF file into a BED file using awk:

awk -v OFS="\t" '{print $1, $4 - 1, $5}' input.gff > output.bed

Subtract 1 from the start position to make it 0-based (thanks to dariober for pointing this out.)

Then use bedtools to look for intersections between the two BED files:

bedtools intersect -a first.bed -b second.bed

[dariober]

Shouldn't it be

awk -v OFS="\t" '{print $1, $4 - 1, $5}' ...

Since start position in gff is 1-based while in bed is 0-based?

[James Ashmore]

Ah yes you are correct, I'll edit my comment, thank you

Answer 2

[Devon Ryan]

It sounds like you're after stringTie, or cufflinks, or similar. Note that if you have a well annotated genome (e.g., human or mouse), then these regions are likely just repeat regions (e.g., LTRs).

Answer 3

[dariober]

As I understand it, you have reads (RNAseq?) mapped to a reference genome and you want to find clusters of reads, right? If so, you could try to scan your bam file with a peak caller devised for ChIP-Seq data, e.g. macs. Alternatively, you could divide your genome in overlapping windows of desired length, say 100bp windows, count reads in each window and consider as clusters windows with read count above a threshold (in the end that's what peak callers do). for this see the bedtools suite.

Once you have the positions of the clusters you can obtain the sequence with a tool like fastaFromBed (in bedtools).

[manekineko]

Not RNA-seq, it is smallRNA-seq, but probably finding the regions of clusters will be better with some peak caller rather than trying to assemble reads...

Answer 4

[Malcolm.Cook]

There are so many variants of your question, it is hard to know exactly what you desire to do.

Trying to literally understand your word "genomic regions occupied with reads", in your coverage map picture thre appear to be six such regions, separated from each other by short regions of zero coverage ("regions NOT occupied by reads").

It that is what you want, then bedtools merge (which works with bam and bed files) will directly answer this question, and identify the coordinates of the six regions.

Docs: http://bedtools.readthedocs.org/en/latest/content/tools/merge.html

Answer 5

[A. Domingues]

If you think these might be piRNA clusters, running proTrac on your data might help you:

proTRAC predicts and analyzes genomic piRNA clusters based on mapped piRNA sequence reads. proTRAC 2.0 and later versions apply a sliding window approach to detect loci that exhibit high sequence read coverage. Subsequently, sequences mapped to these loci are analyzed with respect to typical piRNA and piRNA cluster characteristics to ensure high specificity.

You can find all the scripts needed here.