question to pick values based on conditions using a loop

Question

I have a dataset with following information. Column C1_1 is chromosome number, C1_2 is SNP position and C3 is the p-value. I want to pick the most significant association within a genomic region of 3000 bp:

C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_148000  L002g001    1.1e-3

I want to split first column into two as:

C1_1    C1_2    C2  C3
L01 005000  L002g034    1.5e-12
L01 003000  L001g045    2.3e-6
L02 145000  L004g034    1.1e-2
L02 155000  L002g050    1.1e-4
L02 148000  L002g001    1.1e-3

Then I want to sort subsets such as L01 and L02 separately based on C1_2 column:

C1_1    C1_2    C2  C3
L01 003000  L001g045    2.3e-6
L01 005000  L002g034    1.5e-12
L02 145000  L004g034    1.1e-2
L02 148000  L002g001    1.1e-3
L02 155000  L002g050    1.1e-4

Then I want to pick one value from each subset such as for L02 – I want to pick the most significant value based on C3 column under the difference of 3000. For example, L02-148000 is chosen as shown below as its C3 value is lower than L02-145000, and then the next value L02-155000 is chosen as the difference is more than 3000 than the previous value of L02-148000:

C1_1    C1_2    C2  C3
L01 005000  L002g034    1.5e-12
L02 148000  L002g001    1.1e-3
L02 155000  L002g050    1.1e-4

I have a long dataset with 500,000 rows of values and the subset goes from L01-L24. I think it will need a loop to process these many values to get the final values.

This is just an example, I have shown one for L01 because the difference in L01 203000 and L01 205000 is less than 3000. I have shown two values for L02 because the difference between L02 145000 and L02 148000 is not more than 3000. Whereas the difference between L02 155000 and L02 148000 is more than 3000.

Answer 1

This question is deceptively challenging, but here is a potential solution:

library(tidyverse)

df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_148000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 2.3e-06
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_148000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 3 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     5000 L002g034 1.5e-12
#> 2 L02   148000 L002g001 1.1e- 3
#> 3 L02   155000 L002g050 1.1e- 4

## additional example data
df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_158000  L004g034    1.1e-2
L02_159000  L002g050    1.1e-4
L03_165000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 2.3e-06
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_158000 L004g034 1.1e-02
#> 6 L02_159000 L002g050 1.1e-04
#> 7 L03_165000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 5 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     5000 L002g034 1.5e-12
#> 2 L02   145000 L004g034 1.1e- 2
#> 3 L02   155000 L002g050 1.1e- 4
#> 4 L02   159000 L002g050 1.1e- 4
#> 5 L03   165000 L002g001 1.1e- 3

^{Created on 2023-02-07 with reprex v2.0.2}

Edit

If you have two SNPs with the same C3 value (i.e. the top two rows in the example dataset below) you can select one of them using slice_min(C3, with_ties = FALSE), e.g. the code below will select the 'first' entry when two C3 values are the same within the same group/subgroup:

library(tidyverse)

df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    1.5e-12
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_158000  L004g034    1.1e-2
L02_159000  L002g050    1.1e-4
L03_165000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 1.5e-12
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_158000 L004g034 1.1e-02
#> 6 L02_159000 L002g050 1.1e-04
#> 7 L03_165000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3, with_ties = FALSE) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 4 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     3000 L001g045 1.5e-12
#> 2 L02   145000 L004g034 1.1e- 2
#> 3 L02   155000 L002g050 1.1e- 4
#> 4 L03   165000 L002g001 1.1e- 3

^{Created on 2023-02-08 with reprex v2.0.2}

Answer 2

EDIT: Ignore this answer, it does not address one of your requirements - that of the 3000-long bin. You're going to need some partitioning operations which are not straightforward.

@jared_marmot has cracked it in their solution.

This should be easy to do with dplyr/tidyr. Use

## Code not tested, might need tweaking
library(tidyr)
library(dplyr)
df_split <- df %>% separate(C1, sep = "_", into = c("C1_1","C1_2"))
df_aggr <- df_split %>% select(C1_1, C3) %>% group_by(C1_1) %>% summarise(C3 = min(C3)) %>% inner_join(df_split)

There are two problems that make your 3000 window complicated:

1. The start and end are not fixed: it's not 1 to n in steps of 3000; that'd be a straightforward bin.
2. We cannot even go with the smallest coordinate as the start, as you may have data points at 1500, 2500, 7200, 7600. If those data points exist, a set of 3000-long windows starting at 1500 would be 1500-4500, 4500-7500, etc. making your data go in 3 bins whereas your data points belong in 2 bins. You need windows with independent start coordinates. I've never encountered this requirement in my career even as a puzzle so I'm starting to think this might be an X-Y problem that needs to be addressed at a different level.