[R] Identifying a change in events between bins

Mon Mar 19 16:07:16 CET 2012

Your description was too general for me to know exactly what you want but 
perhaps this will help you solve your own problem
> set.seed(123)
> evtlist = sample(c('fwd','rev'),100,replace=TRUE)
> evtlist
  [1] "fwd" "rev" "fwd" "rev" "rev" "fwd" "rev" "rev" "rev" "fwd" "rev" 
"fwd" "rev"
[14] "rev" "fwd" "rev" "fwd" "fwd" "fwd" "rev" "rev" "rev" "rev" "rev" "rev" 
"rev"
[27] "rev" "rev" "fwd" "fwd" "rev" "rev" "rev" "rev" "fwd" "fwd" "rev" "fwd" 
"fwd"
[40] "fwd" "fwd" "fwd" "fwd" "fwd" "fwd" "fwd" "fwd" "fwd" "fwd" "rev" "fwd" 
"fwd"
[53] "rev" "fwd" "rev" "fwd" "fwd" "rev" "rev" "fwd" "rev" "fwd" "fwd" "fwd" 
"rev"
[66] "fwd" "rev" "rev" "rev" "fwd" "rev" "rev" "rev" "fwd" "fwd" "fwd" "fwd" 
"rev"
[79] "fwd" "fwd" "fwd" "rev" "fwd" "rev" "fwd" "fwd" "rev" "rev" "rev" "fwd" 
"fwd"
[92] "rev" "fwd" "rev" "fwd" "fwd" "rev" "fwd" "fwd" "rev"
> rle(evtlist)
Run Length Encoding
  lengths: int [1:50] 1 1 1 2 1 3 1 1 1 2 ...
  values : chr [1:50] "fwd" "rev" "fwd" "rev" "fwd" "rev" "fwd" "rev" "fwd" 
...
> rle(evtlist)$lengths
[1]  1  1  1  2  1  3  1  1  1  2  1  1  3  9  2  4  2  1 12  1  2  1  1  1 
2  2  1  1
[29]  3  1  1  3  1  3  4  1  3  1  1  1  2  3  2  1  1  1  2  1  2  1
>
see ?rle

Rob
-----Original Message----- 
From: Mark Hills
Sent: Friday, March 16, 2012 3:55 PM
To: r-help at r-project.org
Subject: [R] Identifying a change in events between bins

Hi there,

First off, despite this being my first post here, I have scanned the R help 
forums a lot in the past few months to help with some questions, so a big 
thank you to the community as a whole for being so helpful!

I'm somewhat of an R newbie, and have run up against a problem that I can't 
seem to solve.  If anyone is able to help I would really appreciate it!

I'm looking at a number of events across a chromosome, and have written a 
program that collects them into different bins, based on a specified 
binsize.  The events are directional, either forward or reverse, and a 
chromosome can either be fwd/fwd (all the events fall into the fwd bins), 
rev/rev (all the events fall into the rev bins) or fwd/rev (events are 
evenly split).  In some cases, chromosomes switch from one state to another 
(eg fwd/fwd to fwd/rev).  There are a number of rules that dictate my data. 
First, while there is stochastic variation, the sum of fwd and rev in each 
bin should have approximately the same value. If I were to take the total 
number of events and divide them by the number of bins to get an average 
count per bin, I would expect approximately that value in each bin; in the 
case of fwd/fwd it would be about average number in the fwd column and close 
to zero in the rev column, in rev/rev it would be about the average number 
in the rev column and close to zero in the fwd column, and in fwd/rev it 
would be about half the average number in both.

Hopefully my png attachment worked and you can see an example.  The top plot 
shows fwd reads, the 2nd shows rev reads and the third shows fwd minus rev 
reads.

What I would like to be able to do is to automatically assign regions in 
which the chromosome switches from one state to another. From the graphs 
(and from the read.table output below) you can see that this particular 
chromosome is fwd/fwd from bin 1 to 59, fwd/rev from bin 61 to 73, and 
rev/rev for the remainder of the chromosomes.

These are generated from a read.table that looks like this:

bin  fwd  rev
50  484   2
51  366   4
52  527   6
53  635   2
54  573   6
55  506   4
56  600   6
57  560   2
58  504   2
59  545   0
60  501  68
61  419 223
62  252 109
63  259 138
64  355 189
65  218 125
66  140  57
67   45  31
68  276 144
69  263 152
70  330 193
71  439 204
72  347 207
73   10 611
74    6 619
75    2 578
76    7 372
77    6 436
78    4 373
79    8 417
80    2 276

My question is this:

1. Is there an obvious way to automatically identify these regions?

I am not sure how I can go about scanning previous lines within a read.table 
to find a point at which the values change.  In the above example, I would 
like the program to identify that the fwd graph shifts from ~1x the average 
to ~0.5x the average between bin 61 and 62, and from ~0.5x the average to 
~0x the average between bin 72 and 73. Conversely I'd like to identify the 
rev graph shifting from ~0x average to ~0.5x average between bins 59 and 60, 
and from 0.5x average to 1x average from bin 72 to 73.  Finally, I'd like to 
cross-reference the output from fwd and rev to  only pull out reciprocal 
switches (ie those that occur within 3 bins of each other in both fwd and 
rev data sets).

What I've been trying to gt to work is to generate values based on 0, 0.5 
and 1x the average events, and trying to pull out the range of bins that 
fall into each of those categories (possibly 1 SD higher or lower to account 
for the stochastic variation), but I'm not really sure how to go about that.

2. If I can find a way to identify a shift between bins, is there any way to 
then look in smaller bin sizes across those regions.  The bins shown above 
are for 200,000 bases of DNA.  If my program automatically found an event 
between bin 72 and 73 (14,400,000 bases to 14,600,000), is it possible to 
feed that region into say 50,000 base bins to narrow down the location of 
the event?

I'm not sure if either of these questions is a simple task or incredibly 
time consuming and inappropriate to ask a message board.  If you can offer 
any advice though, I would really appreciate it!

Thanks,

Mark

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

------------------------------------------
Robert W. Baer, Ph.D.
Professor of Physiology
Kirksville College of Osteopathic Medicine
A. T. Still University of Health Sciences
800 W. Jefferson St.
Kirksville, MO 63501
660-626-2322
FAX 660-626-2965 