[R] sliding window analysis with rollapply
R. Michael Weylandt
michael.weylandt at gmail.com
Thu Oct 18 23:42:21 CEST 2012
Hi Li,
I'm afraid I don't have much time to look at this right now -- I'm
sort of swamped for the next 5 days -- so I'm forwarding back to the
list so more eyeballs pass over it. In general, you should always cc
the list so the archives remain complete and so that you get a quicker
response.
On Thu, Oct 18, 2012 at 5:12 PM, Wang, Li <li.wang at ttu.edu> wrote:
> Hi, Michael
>
> Thanks very much for your reply!
>
> My main concern is to calculate average Fst within a window. The window size, say 1000,000, means from position 1-1000,000. In this following case, no position within this window. But if the window is 1000,000 -- 2000,000, then the former 9 positions were within the window, and we can calculated the average Fst value. That is what I mean by rolling Fst again pos. I donot know if my intention could be achieved in R.
>
> pos Fst
> 264643 0.0845204537
> 264673 0.0060222786
> 264735 0.0865209872
> 265951 -0.0038305848
> 267707 0.0865209872
> 268666 0.0185608883
> 268741 0.0493999072
> 269643 0.0424427481
> 269771 0.1294509289
> 643220 -0.0160247705
> 643956 0.0106426936
> 644001 0.0411787797
> 644170 0.0302527151
> 644252 0.0067342624
> 644730 0.0073714424
> 645180 0.0331635458
> 648070 0.0919606246
> 648193 -0.0077790431
> 648477 0.0331635458
> 649139 0.0302527151
> 649536 0.0240243946
> 649537 0.0060054306
> 649810 0.0246680573
> 649860 0.0302527151
> 649861 0.0194870768
> 650399 0.0331635458
> 650427 0.0380601183
> 650602 0.0331635458
> 652468 0.0246680573
> 653322 0.039720035
> 653477 0.0246680573
> 653539 0.0380601183
> 653808 -0.0362652864
> 653839 0.0589780782
> 653863 0.0589780782
> 653979 0.1080478911
> 653998 0.1080478911
> 654016 0.0163679102
> 654104 0.0167725359
> 654156 0.1436247159
> 715529 -0.0750541068
> 715621 -0.1207825145
> 716629 -0.0696369342
> 719352 -0.0188278314
> 719359 -0.0188278314
> 719606 0.0355536571
> 720851 -0.0188278314
> 720907 0.0073714424
> 720948 -0.0661275377
> 721436 -0.1073576949
> 723537 0.0110375276
> 724027 0.0134550135
> 724119 -0.030982485
> 724725 0.1468391016
> 728718 -0.1216726698
> 728727 -0.043005885
> 730714 0.1111296269
> 731410 0.0239519768
> 732418 0.0676868746
>
> Sorry for my stupid question. I am very willing to be educated.
>
> Best regards
> Li
>
>
> ________________________________________
> From: R. Michael Weylandt [michael.weylandt at gmail.com]
> Sent: Wednesday, October 17, 2012 7:51 AM
> To: Wang, Li
> Cc: r-help at r-project.org
> Subject: Re: [R] sliding window analysis with rollapply
>
> On Tue, Oct 16, 2012 at 11:17 PM, Wang, Li <li.wang at ttu.edu> wrote:
>> Dear List members
>>
>> I want to do the sliding window analysis of some specific values. Here is my code:
>>
>> require(zoo)
>> dat <- read.table("chr1.txt", header = TRUE, sep="\t")
>> dat2 <- cbind(dat[1,3]) #The first column is also important. It represents the position of the site on the chromosome.
>> TS <- zoo(c(dat2))
>
> Two Style comments:
>
> 1) Perhaps you just want read.zoo()
>
> 2) No need to have the c() when you aren't combining things.
>
>> a <- rollapply(TS, width=1000000, by=200000, FUN=mean, align="left") #The third column values should roll against the first column values. I might be wrong here with the code.
>
> 3) rollmean() directly could also work here, but I'm not at all sure
> what you mean by "the third column values should roll against the
> first column values"
>
>> plot1 <- plot.default(a, col="red", type="l") #It returns me some errors for this command.
>
> 4) Generally better to call the generic plot() instead of the
> specialized plot.default()
>
> Michael
>
>>
>> Error in plot.window(...) : need finite 'xlim' values
>> In addition: Warning messages:
>> 1: In min(x) : no non-missing arguments to min; returning Inf
>> 2: In max(x) : no non-missing arguments to max; returning -Inf
>> 3: In min(x) : no non-missing arguments to min; returning Inf
>> 4: In max(x) : no non-missing arguments to max; returning -Inf
>>
>> Part of the original data looks like as follows:
>> pos Fit Fst Fis
>> 12794 0.928380160665041 -0.0877263098996843 0.934156378600823
>> 12816 0.901040947621283 0.0382417096943425 0.897106109324759
>> 12821 0.901040947621283 0.0382417096943425 0.897106109324759
>> 12827 0.901040947621283 0.0382417096943425 0.897106109324759
>> 12855 0.909446752933768 0.00347893012209298 0.909130624726955
>> 13324 0.498125727088629 0.0975914515920549 0.443850267379679
>> 13338 0.53827858671811 0.103194947279406 0.485148514851485
>> 13339 -0.059306330783135 -0.0198369843236668 -0.038701622971286
>> 13379 0.610641507226514 0.175889055559966 0.527541169789892
>> 13381 0.144593428517224 -0.0343016958072985 0.172962226640159
>> 13390 0.216526396327467 0.179801071155318 0.0447761194029849
>> 13454 0.498125727088629 0.0975914515920549 0.443850267379679
>> 13457 0.309860615307135 -0.0501399771889025 0.342812006319115
>> 13462 0.536517915763086 -0.0302077737766018 0.550108147080029
>>
>> That is, Fst against pos.
>>
>> I am very willing to be educated concerning the mistakes I made for the code.
>>
>> All the best
>> Li
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> 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.
More information about the R-help
mailing list