Как перегруппировать строки на основе значений корреляции

У меня есть один фрейм данных под названием snps генетических вариантов:

             ID Group
 1: 1:12345:A:G     1
 2: 1:12346:T:C     1
 3: 1:23457:A:G     1
 4:  3:1234:A:G     2
 5: 3:12345:A:G     2
 6: 1:99991:A:T     3
 7: 1:99992:C:T     3
 8: 1:99993:A:G     3
 9: 1:99994:T:C     3
10: 1:99995:A:G     3
11:   4:777:A:G     4 
12:  4:7778:T:C     4 
13:  4:7774:A:T     4 
14:  4:7771:C:G     4

Затем у меня есть еще один файл данных под названием ld, который измеряет корреляцию некоторых из этих вариантов друг с другом:

           SNP    lead_snp  R2
1: 1:12346:T:C 1:12345:A:G 0.6
2: 1:23457:A:G 1:12346:T:C 0.1
3: 3:12345:A:G  3:1234:A:G 0.5
4:  3:1234:A:G 3:12345:A:G 0.5
5: 1:99991:A:T 1:99992:C:T 0.2
6: 1:99991:A:T 1:99993:A:G 0.7
7: 1:99994:T:C 1:99991:A:T 0.1
8: 1:99992:C:T 1:99994:T:C 0.6
 9:   4:777:A:G  4:7778:T:C 0.7
10:  4:7774:A:T  4:7771:C:G 0.8

Я собираюсь переназначить уже существующие группы в snps$Group в зависимости от того, коррелируют ли какие-либо snps/варианты <0,4.

По сути, я хочу выполнить:

В каждой уже существующей группе найдите любые snps, которые имеют корреляцию только <0,4 с любым другим snp в этой группе.
Назначьте этот snp новому номеру группы.
Если остальные snps либо отсутствуют в фрейме данных ld (у них нет показателей R2), либо имеют R2 >0,4 с любым другим snps, оставьте их в этой группе.

Прямо сейчас мой код для этого:

reassign_groups <- function(df, ld, threshold = 0.4) {
  df <- df %>% arrange(Group)
  new_group_id <- max(df$Group, na.rm = TRUE) + 1
  low_r2_snps <- data.table()  # Data table to store SNPs with <0.4 R2
  
  for (group_id in unique(df$Group)) {
    snps_in_group <- df[df$Group == group_id, ]
    n <- nrow(snps_in_group)
    
    for (i in 1:n) {
      for (j in (i+1):n) {
        snp1 <- snps_in_group$ID[i]
        snp2 <- snps_in_group$ID[j]
        
        ld_check <- ld[(lead_snp == snp1 & SNP == snp2) | (lead_snp == snp2 & SNP == snp1)]
        if (nrow(ld_check) > 0 && any(ld_check$R2 < threshold)) {
          df$Group[df$ID == snp2] <- new_group_id
          
          # Store SNPs with <0.4 R2 in low_r2_snps
          low_r2_snps <- rbind(low_r2_snps, data.table(snp1 = snp1, snp2 = snp2, R2 = ld_check$R2))
          
          new_group_id <- new_group_id + 1
        }
      }
    }
  }
  return(list(updated_df = df, low_r2_snps = low_r2_snps))
}

# Reassign groups based on LD criteria and get SNPs with <0.4 R2
result <- reassign_groups(snps, ld)

new_groups <- result$updated_df
low_r2_snps <- result$low_r2_snps

Однако это не совсем правильно. Для моего примера выводится:

             ID Group
 1: 1:12345:A:G     1
 2: 1:12346:T:C     1
 3: 1:23457:A:G     4
 4:  3:1234:A:G     2
 5: 3:12345:A:G     2
 6: 1:99991:A:T     3
 7: 1:99992:C:T     5
 8: 1:99993:A:G     3
 9: 1:99994:T:C     6
10: 1:99995:A:G     3

Ожидаемый порядок этих данных в целом (редактирование: включая новый крайний случай для группы 4):

            ID Old_group New_group
 1: 1:12345:A:G         1         1
 2: 1:12346:T:C         1         1
 3: 1:23457:A:G         1         5
 4:  3:1234:A:G         2         2
 5: 3:12345:A:G         2         2
 6: 1:99991:A:T         3         3
 7: 1:99992:C:T         3         3
 8: 1:99993:A:G         3         3
 9: 1:99994:T:C         3         3
10: 1:99995:A:G         3         3
11:   4:777:A:G         4         4
12:  4:7778:T:C         4         4
13:  4:7774:A:T         4         4
14:  4:7771:C:G         4         4

Только 1:23457:A:G был присвоен новый номер группы, так как он имеет <0,4 со всеми остальными snp в старой группе. Некоторые другие snps имеют <0,4 для некоторых snps, но >0,4 для других в их группе (поэтому они остаются в группе).

Как я могу исправить свой код, чтобы достичь этого?

Обновлено: обновленные данные примера:


ld <- structure(list(SNP = c("1:12346:T:C", "1:23457:A:G", "3:12345:A:G", 
                             "3:1234:A:G", "1:99991:A:T", "1:99991:A:T", "1:99994:T:C", "1:99992:C:T", 
                             "4:777:A:G", "4:7774:A:T"), lead_snp = c("1:12345:A:G", "1:12346:T:C", 
                                                                      "3:1234:A:G", "3:12345:A:G", "1:99992:C:T", "1:99993:A:G", "1:99991:A:T", 
                                                                      "1:99994:T:C", "4:7778:T:C", "4:7771:C:G"), R2 = c(0.6, 0.1, 
                                                                                                                         0.5, 0.5, 0.2, 0.7, 0.1, 0.6, 0.7, 0.8)), row.names = c(NA, -10L
                                                                                                                         ), class = c("data.table", "data.frame"))

snps <-structure(list(ID = c("1:12345:A:G", "1:12346:T:C", "1:23457:A:G", 
                             "3:1234:A:G", "3:12345:A:G", "1:99991:A:T", "1:99992:C:T", "1:99993:A:G", 
                             "1:99994:T:C", "1:99995:A:G", "4:777:A:G", "4:7778:T:C", "4:7774:A:T", 
                             "4:7771:C:G"), Group = c(1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 
                                                      3L, 4L, 4L, 4L, 4L)), row.names = c(NA, -14L), class = c("data.table", 
                                                                                                               "data.frame"))

Вам следует проверить stats::kmeans для кластеризации групп. Это проще реализовать, чем вашу методологию группировки.

— 21.05.2024 17:11

Спасибо за ваш ответ. Метод, который я собираюсь здесь применить, необходим для группировки вариантов по локусам. Я назвал это здесь корреляцией для доступного понимания, но рассчитанный здесь R2 измеряет неравновесие по сцеплению между генетическими вариантами - и это затем используется для группировки вариантов в генетических исследованиях. Обычно люди используют инструменты командной строки для геномики, предназначенные для выполнения этой задачи, но у меня нет подходящих типов данных, поэтому я пытаюсь написать версию на R.

— 21.05.2024 17:20

r data.table tidyverse

21.05.2024 17:01

Стоит ли изучать PHP в 2026-2027 годах?

Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...

Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией

В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.

Приемы CSS-макетирования - floats и Flexbox

Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...

Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest

В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...

Концепция локализации и ее применение в приложениях React ⚡️

Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...

Пользовательский скаляр GraphQL

Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...

151

Перейти к ответу Данный вопрос помечен как решенный

Ответы 3

Итак, я до сих пор не совсем понимаю, что вы делаете, но я опробовал этот способ создания групп с нуля, начиная с наиболее коррелирующих генов.

# stack all the correlations we saw, so just the ID column has all options
correlations <- rbind(ld[,.(ID = SNP, ID_other = lead_snp, R2)], 
                      ld[,.(ID = lead_snp, ID_other = SNP, R2)])
# snap them to your list and order starting with the highest correlations
all_correlations <- merge(snps, correlations)[order(R2, decreasing = T)]
all_correlations[,new_Group:=0] # initialize a new group column
group_counter = 0 
seen <- list()
seen_counter <- 0
for(i in 1:nrow(all_correlations)){
  if (!(all_correlations$ID[i] %in% unlist(seen)) & 
     !(all_correlations$ID_other[i] %in% unlist(seen)) |
     all_correlations$R2[i] < 0.4){
    group_counter = group_counter + 1
    all_correlations$new_Group[i] = group_counter 
    seen_counter = seen_counter + 1
    seen[[i]] <- all_correlations$ID[i]
    seen_counter = seen_counter + 1
    seen[[i]] <- all_correlations$ID_other[i]
  } else {
    all_correlations$new_Group[i] = group_counter 
  }
}
output <- all_correlations[!duplicated(ID), 
                           .(ID, Group, 
                             new_Group = as.integer(factor(new_Group)))][order(Group)]

Я уверен, что есть лучший способ сделать то, что вы хотите, но, похоже, это работает.

output
#            ID Group new_Group
#1: 1:12345:A:G     1         2
#2: 1:12346:T:C     1         2
#3: 1:23457:A:G     1         5
#4: 3:12345:A:G     2         4
#5:  3:1234:A:G     2         4
#6: 1:99991:A:T     3         1
#7: 1:99993:A:G     3         1
#8: 1:99992:C:T     3         3
#9: 1:99994:T:C     3         3

22.05.2024 07:03

Ответ принят как подходящий

Новый ответ после изменения требований

В этом случае нет необходимости в igraph, нам просто нужно получить максимальный R2, который имеет любой SNP в группе с любым SNP, если меньше 0,4, затем установить новую группу.

ldg <- merge(ld, snps[, .(SNP = ID, Group)], all.x = TRUE)
ldg <- unique(rbind(ldg[, .(ID = SNP, R2)],
                    ldg[, .(ID = lead_snp, R2)]))
res <- merge(snps, ldg, all.x = TRUE)
setorder(res, rn)
res <- res[, .(R2max = max(R2)), by = .(ID, Group, rn)]

#if R2 missing set to 1
res[ is.na(R2max), R2max := 1 ]
# if R2 < 4 set group to NA
res[ R2max < 0.4, Group := NA ]
# get max group ID
gmax <- max(res$Group, na.rm = TRUE)
# set new Group ID for NAs
res[ is.na(Group), Group := gmax + .I ]
res
#             ID Group rn R2max
# 1: 1:12345:A:G     1  1   0.6
# 2: 1:12346:T:C     1  2   0.6
# 3: 1:23457:A:G     5  3   0.1 ### <- new group assigned as 5
# 4:  3:1234:A:G     2  4   0.5
# 5: 3:12345:A:G     2  5   0.5
# 6: 1:99991:A:T     3  6   0.7
# 7: 1:99992:C:T     3  7   0.6
# 8: 1:99993:A:G     3  8   0.7
# 9: 1:99994:T:C     3  9   0.6
#10: 1:99995:A:G     3 10   1.0
#11:   4:777:A:G     4 11   0.7
#12:  4:7778:T:C     4 12   0.7
#13:  4:7774:A:T     4 13   0.8
#14:  4:7771:C:G     4 14   0.8

Старый ответ

Использование членства в igraph на основе высокого ld:

#Example data:
ld <- structure(list(SNP = c("1:12346:T:C", "1:23457:A:G", "3:12345:A:G", 
                             "3:1234:A:G", "1:99991:A:T", "1:99991:A:T", "1:99994:T:C", "1:99992:C:T"
), lead_snp = c("1:12345:A:G", "1:12346:T:C", "3:1234:A:G", "3:12345:A:G", 
                "1:99992:C:T", "1:99993:A:G", "1:99991:A:T", "1:99994:T:C"), 
R2 = c(0.6, 0.1, 0.5, 0.5, 0.2, 0.7, 0.1, 0.6)), row.names = c(NA, 
                                                               -8L), class = c("data.table", "data.frame"))

snps <-structure(list(ID = c("1:12345:A:G", "1:12346:T:C", "1:23457:A:G", 
                             "3:1234:A:G", "3:12345:A:G", "1:99991:A:T", "1:99992:C:T", "1:99993:A:G", 
                             "1:99994:T:C", "1:99995:A:G"), Group = c(1L, 1L, 1L, 2L, 2L, 
                                                                      3L, 3L, 3L, 3L, 3L)), row.names = c(NA, -10L), class = c("data.table", 
                                                                                                                               "data.frame"))


library(igraph)
# get clusters where R2 is more than 0.4
g <- graph_from_data_frame(ld[ R2 > 0.4, ])
plot(g)

gGroup <- setNames(stack(components(g)$membership), c("membership", "ID"))

# add rownumber to reorder after merge
snps[, rn := .I ]
res <- merge(snps, gGroup, all.x = TRUE)
setorder(res, rn)

# if there is no ld, then keep same group
res[ !ID %in% unlist(ld[, 1:2 ]), newGroup := paste0("a_", Group) ]
# no membership means low ld, make new group
res[ is.na(newGroup) & is.na(membership), newGroup := paste0("b_", .I) ]
# update group based on igraph membership
res[ is.na(newGroup) & !is.na(membership),
     newGroup := ifelse(Group == membership, paste0("a_", membership),
                        paste0("c_", membership)) ]

#convert to integer groups
res[, newGroup := as.integer(as.factor(newGroup)) ]
res
#              ID Group rn membership newGroup
#  1: 1:12345:A:G     1  1          1        1
#  2: 1:12346:T:C     1  2          1        1
#  3: 1:23457:A:G     1  3         NA        4
#  4:  3:1234:A:G     2  4          2        2
#  5: 3:12345:A:G     2  5          2        2
#  6: 1:99991:A:T     3  6          3        3
#  7: 1:99992:C:T     3  7          4        5
#  8: 1:99993:A:G     3  8          3        3
#  9: 1:99994:T:C     3  9          4        5
# 10: 1:99995:A:G     3 10         NA        3

Большое вам спасибо за это. Это очень близко к тому, что мне нужно, и очень быстро работает со всеми моими данными. Однако мне нужно было немного обновить свой вопрос и ожидаемый результат. Я добавил награду к вопросу, можно ли перенацелить этот подход на новые требования (включение snps в новые группы только тогда, когда они составляют <0,4 R2 для каждого отдельного snp в этой группе, и если это будет только тогда, когда создается новая группа).

— 23.05.2024 20:14

22.05.2024 14:03

Наверное, ты можешь попробовать это

lvls <- union(snps$ID, unlist(ld[, .(SNP, lead_snp)]))
m <- matrix(NA,
    length(lvls),
    length(lvls),
    dimnames = list(lvls, lvls)
)
m <- with(ld, {
    m[SNP, lead_snp] <- R2
    m[lead_snp, SNP] <- R2
    m
})
setcolorder(
    snps[
        ,
        c(
            .SD,
            .(NewGroup = nafill(+(rowMeans(m[ID, ID] >= 0.4, TRUE) == 0),
                fill = 0
            ))
        ),
        Group
    ][
        ,
        NewGroup := pmax(Group, NewGroup * (max(Group) + cumsum(NewGroup)))
    ], "ID",
    before = 1
)[]

который дает

            ID Group NewGroup
         <char> <int>    <int>
 1: 1:12345:A:G     1        1
 2: 1:12346:T:C     1        1
 3: 1:23457:A:G     1        5
 4:  3:1234:A:G     2        2
 5: 3:12345:A:G     2        2
 6: 1:99991:A:T     3        3
 7: 1:99992:C:T     3        3
 8: 1:99993:A:G     3        3
 9: 1:99994:T:C     3        3
10: 1:99995:A:G     3        3
11:   4:777:A:G     4        4
12:  4:7778:T:C     4        4
13:  4:7774:A:T     4        4
14:  4:7771:C:G     4        4

23.05.2024 23:34