Question: Find overlapping coordinates between a bed file and multiple files in order and associate the coordinates with the original coordinate in tsv
0
gravatar for skim
23 months ago by
skim30
skim30 wrote:

Hello, I have posted a question here but I found out that bedtools intersect was not sufficient for me. I need to make a tab separated table in format like below

I have one original bed file and 125 cell line bed files that I have to check overlaps and list the overlapping coordinates as a new file in the following format:

original bed file: an ordinary bed6 file

chromosome start end nametag score strand(+/-)

ex) chr1 1000001 1000095 RBP7 70.69 +

cellline bed file: also an ordinary bed6 or bed8 file

what i have to make:

I have to make a new tsv file that appends values to the right of each bed coordinate.

chromosome start end nametag score strand(+/-) cellline1 cellline2 ....etc..... cellline125 Total_celline_overlap_count

ex) chr1 1000001 1000095 RBP7 70.69 + chr1:1000001-1000004 None ...etc.... chr1:1000003-1000015 60

basically, I have to write down in every line 1) whether this original coordinate has an overlap at the given cell line and 2) if there is, write down the overlaps.

Therefore there are 6+125+1 columns in this file.

I found out that bedtools intersect does not provide this function and even if I could code in python to do this, the time consumption is very inefficient if i use the files generated in bedtools intersect. I used glob package to search the 125 files in my directory and I used csv package to make my custom tsv file but they were insufficient to do this job with bedtools intersect. Can there be any help for me? Thank you very much.

sequence next-gen bed overlap • 1.3k views
ADD COMMENTlink modified 23 months ago by Alex Reynolds28k • written 23 months ago by skim30
1

All you need is original bedtools intersect and linux join command.

Algorithm goes like that:

for i in cellLineBed {
    # you need wao to get all reports in -a 
    # possible to overcome this with linux join
    bedtoolsIntersect -a originalBedFile -b ${i}_cellLineBed -wao > result_${i}
}
join result_* by originalBedFile (chromosome start end nametag)

I haven't used join in a long time and don't remember it's syntax, but it should do what you want

ADD REPLYlink modified 23 months ago • written 23 months ago by PoGibas4.8k
1
gravatar for Alex Reynolds
23 months ago by
Alex Reynolds28k
Seattle, WA USA
Alex Reynolds28k wrote:

BEDOPS does this efficiently.

First, make a file called union.bed that contains a custom ID built from the chromosome, start and stop of each element. Then map that file to original.bed, retrieving IDs of overlapping elements (--echo-map-id-uniq) and the total count of overlaps (--count):

$ bedops --everything cellLine*.bed | awk '{ print $1"\t"$2"\t"$3"\t0\t"$1":"$2"-"$3; }' > union.bed
$ bedmap --echo --echo-map-id-uniq --count --delim '\t' --multidelim '\t' original.bed union.bed > answer.bed

An even more efficient approach would use a process substitution, which avoids creating the intermediate file union.bed:

$ bedmap --echo --echo-map-id-uniq --count --delim '\t' --multidelim '\t' original.bed <(bedops --everything cellLine*.bed | awk '{ print $1"\t"$2"\t"$3"\t0\t"$1":"$2"-"$3; }') > answer.bed

Use sort-bed to sort your input BED files original.bed and cellLine*.bed before using bedops and bedmap, if the inputs are not sorted or if their sort order is unknown, e.g.:

$ sort-bed original.unsorted.bed > original.bed
ADD COMMENTlink modified 23 months ago • written 23 months ago by Alex Reynolds28k

Is this only possible after sorting? It would be better off for my post-processing to just preserve the order and just have more processing time.

ADD REPLYlink written 23 months ago by skim30
1

If you need to sort inputs, yes, you would need to sort inputs. Using sorted input is why BEDOPS tools run efficiently. On the plus side, you only have to sort a file once, as the tools bedops, bedmap, etc. write sorted output.

ADD REPLYlink modified 23 months ago • written 23 months ago by Alex Reynolds28k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1096 users visited in the last hour