4
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ Is your result C1_2 value for L01 supposed to be 005000 or 205000? $\endgroup$
    – Ram RS
    Feb 6, 2023 at 21:22
  • 1
    $\begingroup$ I think that is a typo @RamRS L01 should be 205000 (lower table), either that or its not possible to solve. $\endgroup$
    – M__
    Feb 6, 2023 at 21:24
  • $\begingroup$ If @RamRS answers solves points 1 to 3, I'd recommend accepting this as the answer and then asking a further question on the 3000 threshold. In Python you'd need to construct a new column based the differences between C1_2 .. well thats how I'd do it anyway. $\endgroup$
    – M__
    Feb 6, 2023 at 21:34

2 Answers 2

4
$\begingroup$

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

$\endgroup$
3
  • 2
    $\begingroup$ Wonderful use of cumsum + lag. Very nifty solution! $\endgroup$
    – Ram RS
    Feb 7, 2023 at 18:34
  • 1
    $\begingroup$ Thank you @jared_mamrot! I tried using your code. It gives the correct result with my data. Except one thing: if there are two SNPs with different positions in C1_2 but same C3 values within the 3000 window, it lists all those associations. I have not thought about this situation earlier. $\endgroup$ Feb 7, 2023 at 19:43
  • 1
    $\begingroup$ Sounds like an easy fix @user1567654 - instead of slice_min(C3) try slice_min(C3, with_ties = FALSE) - that should do the trick, but if not, leave a comment and I'll look into it :) $\endgroup$ Feb 7, 2023 at 22:27
1
$\begingroup$

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.
$\endgroup$
2
  • 1
    $\begingroup$ Oh wait you have that 3000 condition. That's not easy to implement. You'll need to create a new column with a unique value per 3000-bin that can be used to group bins together, which is not simple as the 3000's bounds are not fixed. $\endgroup$
    – Ram RS
    Feb 6, 2023 at 21:29
  • $\begingroup$ I think this is basically the answer and the 3000 threshold is a separate question. $\endgroup$
    – M__
    Feb 6, 2023 at 21:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.