Splitting genomic intervals based on chromosome arm

Question

[jtull]

I have a list of genomic intervals of interest and I am interested in calculating fraction of chromosome arm they make up.

For example,

chr     start      stop
1         50,000,000     100,000,000
1         120,000,000     150,000,000

And using chromosome arm size information (e.g. chr1p spans region 0-123,400,000, chr1q spans between 123,400,000 and 248,956,422), For this, first I need to split the intervals using chromosome arm boundaries, like:

chr  start   stop arm
1     50,000,000     100,000,000   p
1     120,000,000     123,400,000    p
1      123,400,000      150,000,000   q

Then I will merge the ones on the same chromosome and arm and calculate the fraction of the chromosome arm they make up. Do you have any suggestions on how to split the intervals? Is there an easy function/way or I need to write a script? Thanks.

Answer 1

[Jianyu]

I think you can use bedtools intersect to compare your intervals with each chromosome arm, then annotate each output with the arm used.

A simple bash example:

old_IFS=$IFS
IFS=$'\n'
for i in `cat arm.bed`; do
        name=`echo $i | cut -f 1`
        echo $i | sed 's/\([0-9]*\)[a-z]/\1/g' | bedtools intersect -a interval.bed -b - > interval.$name
        sed -i 's/$/\t'$name'/g' interval.$name
        cat interval.$name >> output.interval
done
IFS=${old_IFS}

interval.bed:

1   50000000    100000000
1   120000000   150000000

arm.bed:

1p  0   123400000
1q  123400000   248956422

output.interval:

1   50000000    100000000   1p
1   120000000   123400000   1p
1   123400000   150000000   1q

Note: I assume the chromosome arm name = chromosome name (which consists of number) + a character

[jtull]

Thank you, it works just as I wanted! Just a note, maybe add rm interval.*[pq] in the end of the bash script to clean up the temporary interval files.