R运行非时间数据的平均值 [英] R running average for non-time data

问题描述

这是我现在的情节。

这是从这个代码生成：

pre $ g $ p $ ggplot（data1，aes（x = POS，y = DIFF，color = GT））+
geom_point（）+
facet_grid（〜CHROM，scales =free_x，space =free_x）+
主题（strip.text.x = element_text（size = 40），
strip.background = element_rect（color ='lightblue'，fill ='lightblue'），
legend.position =top，
legend.title = element_text（size = 40，color = ），
legend.text = element_text（size = 40），
legend.key.size = unit（2.5，cm））+
guides（fill = guide_legend（title .position =top，
title =Legend：GT ='REF'+'ALT'），
shape = guide_legend（override.aes = list（size = 10）））+
scale_y_log10（breaks = trans_breaks（log10，function（x）10 ^ x，n = 10））+
scale_x_continuous（breaks = pretty_breaks（n = 3））+
geom_line（stat = hline，
yintercept =mean，
size = 1）

最后一行，geom_line为每个面板创建平均线。

但是现在我想要在每个面板中具有更具体的运行平均值。

即如果panel1（'chr01'）的x轴范围是从0到100,000,000，那么我希望每个1,000,000范围的平均值。

mean1 = mean（x = 0到x = 1,000,000）

mean2 =平均值（x = 1,000,001至x = 2,000,000）

解决方案

提供运行平均值的一种方法是使用 geom_smooth（）使用 loess 局部回归方法。为了展示我提出的解决方案，我使用R函数创建了一个假基因组数据集。您可以调整 span 参数 geom_smooth 以使运行平均值更平滑（更接近1.0）或更粗糙（更接近1 /数据点数）。

 ＃创建示例数据。 
 set.seed（27182）
 
 y1 = rnorm（10000）+ 
c（rep（0，1000），dnorm（seq（-2，5，length.out = 8000））* 3，rep（0,1000））
 y2 = c（rnorm（2000），rnorm（1000，mean = 1.5），rnorm（1000，mean = -1，sd = 2），
 rnorm（2000，sd = 2））
 y3 = rnorm（4000）
 pos = c（sort（runif（10000，min = 0，max = 1e8）），
 （runif（6000，min = 0，max = 6e7）），
 sort（runif（4000，min = 0，max = 4e7）））
 chr = rep（c（chr01， chr02，chr03），c（10000,6000,4000））
 
 data1 = data.frame（CHROM = chr，POS = pos，DIFF = c（y1，y2，y3） ）
 
＃剧情。 
p = ggplot（data1，aes（x = POS，y = DIFF））+ 
 geom_point（alpha = 0.1，size = 1.5）+ 
 geom_smooth（color =darkgoldenrod1，size = 1.5，method =loess，degree = 0，
 span = 0.1，se = FALSE）+ 
 scale_x_continuous（breaks = seq（1e7，3e8，1e7），
 labels = paste seq（10,300,10）），expand = c（0,0））+ 
 xlab（Position，Megabases）+ 
 theme（axis.text.x = element_text（size = 8） ））+ 
 facet_grid（。〜CHROM，scales =free，space =free）
 
 ggsave（filename =plot_1.png，plot = p，width = 10 ，height = 5，dpi = 150）
  

This is the plot I'm having now.

It's generated from this code:

ggplot(data1, aes(x=POS,y=DIFF,colour=GT)) + 
  geom_point() +
  facet_grid(~ CHROM,scales="free_x",space="free_x") + 
  theme(strip.text.x = element_text(size=40),
        strip.background = element_rect(color='lightblue',fill='lightblue'),
        legend.position="top",
        legend.title = element_text(size=40,colour="lightblue"),
        legend.text = element_text(size=40),
        legend.key.size = unit(2.5, "cm")) +
  guides(fill = guide_legend(title.position="top",
                             title = "Legend:GT='REF'+'ALT'"),
         shape = guide_legend(override.aes=list(size=10))) +
  scale_y_log10(breaks=trans_breaks("log10", function(x) 10^x, n=10)) + 
  scale_x_continuous(breaks = pretty_breaks(n=3)) +
  geom_line(stat = "hline",
            yintercept = "mean",
            size = 1)

The last line, geom_line creates the mean line for each panel.

But now I want to have the more specific running average inside each panel.

i.e. If panel1('chr01') has x-axis range from 0 to 100,000,000, I would want to have the mean value for each 1,000,000 range.

mean1 = mean(x=0 to x=1,000,000)

mean2 = mean(x=1,000,001 to x=2,000,000)

like that.

解决方案

One way to provide a running mean is with geom_smooth() using the loess local regression method. In order to demonstrate my proposed solution, I created a fake genomic dataset using R functions. You can adjust the span parameter of geom_smooth to make the running mean smoother (closer to 1.0) or rougher (closer to 1/number of data points).

# Create example data.
set.seed(27182)

y1 = rnorm(10000) + 
     c(rep(0, 1000), dnorm(seq(-2, 5, length.out=8000)) * 3, rep(0, 1000))
y2 = c(rnorm(2000), rnorm(1000, mean=1.5), rnorm(1000, mean=-1, sd=2), 
       rnorm(2000, sd=2))
y3 = rnorm(4000)
pos = c(sort(runif(10000, min=0, max=1e8)),
        sort(runif(6000,  min=0, max=6e7)),
        sort(runif(4000,  min=0, max=4e7)))
chr = rep(c("chr01", "chr02", "chr03"), c(10000, 6000, 4000))

data1 = data.frame(CHROM=chr, POS=pos, DIFF=c(y1, y2, y3))

# Plot.
p = ggplot(data1, aes(x=POS, y=DIFF)) +
    geom_point(alpha=0.1, size=1.5) +
    geom_smooth(colour="darkgoldenrod1", size=1.5, method="loess", degree=0, 
        span=0.1, se=FALSE) +
    scale_x_continuous(breaks=seq(1e7, 3e8, 1e7), 
        labels=paste(seq(10, 300, 10)), expand=c(0, 0)) +
    xlab("Position, Megabases") +
    theme(axis.text.x=element_text(size=8)) +
    facet_grid(. ~ CHROM, scales="free", space="free")

ggsave(filename="plot_1.png", plot=p, width=10, height=5, dpi=150)

这篇关于R运行非时间数据的平均值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！