Find site feature according to GFF3 file

Question

[biolab]

Dear all,

I have a list of site information exemplified as below. Also I have a GFF3 file. Now I need to find where those sites locate, for example, my ideal output is shown at bottom. Could anyone help to give me suggestions or solutions? I much appreciated your helps. THANKS A LOT!

Site file:

Site No.   Chr  Start   End
1          1    3000    3100
2          1    3200    3280

GFF3 file:

Chr1 MSU_osa1r7 gene 2903 10817 . + . ID=LOC_Os01g01010;Name=LOC_Os01g01010
Chr1 MSU_osa1r7 mRNA 2903 10817 . + . ID=LOC_Os01g01010.1;Name=LOC_Os01g01010.1
Chr1 MSU_osa1r7 exon 2903 3268 . + . ID=LOC_Os01g01010.1:exon_1
Chr1 MSU_osa1r7 intron 3269 3300 . + . ID=LOC_Os01g01010.1:intron_1

Ideal output:

Site No.  Chr  Start   End    Pos
1         1    3000    3100   ID=LOC_Os01g01010.1:exon_1
2         1    3200    3280   ID=LOC_Os01g01010.1:exon_1;ID=LOC_Os01g01010.1:intron_1

[PoGibas]

Does your site file contain any chromosome information?

[biolab]

Hi Pgibas, yes, I forget to add it.

[Devon Ryan]

A bit of coding with R and the GenomicRanges packages should make this doable.

Answer 1

[PoGibas]

Modify your Site and GFF3 files using awk (change order of columns - make bed files).

Use bedtools intersect to get wanted output.

[biolab]

Thanks a lot, Pgibas and Devon, I need a bit further guidance. Regarding to the bedtools method, I follow the below steps, unfortunately got an error : Error: unable to open file or unable to determine types for file f1.bed

Step1 - Prepare site file and gff3 file, both of which should be in the same format.

Site file: f1.bed

Chr1    3000    3100   +
Chr1    3200    3280   +

GFF3 file: f2.bed

Chr1 ID=LOC_Os01g01010;Name=LOC_Os01g01010  2903 10817  +
Chr1 ID=LOC_Os01g01010.1;Name=LOC_Os01g01010.1  2903 10817  +
Chr1 ID=LOC_Os01g01010.1:exon_1  2903 3268  +
Chr1 ID=LOC_Os01g01010.1:intron_1  3269 3300  +

Step2 -- Use bedtools intersect to find their overlapped region.

I use this command: bedtools intersect -a f1.bed -b f2.bed

[PoGibas]

f2.bed isn't a bed file (take a look at this explanation). Fields in the bed file should be: chrom start end.

Change your f2 file into something like this:

Chr1  2903 10817 + ID=LOC_Os01g01010;Name=LOC_Os01g01010

[biolab]

Hi Pgibas, Thank you very much.