Alzheimer's disease, GWAS & eQTL 분석 (2) GWAS 정리

우선 받아온 GWAS summary data를 먼저 정리하겠습니다.

 

 

library(tidyverse)
library(data.table)
library(stringr)
library(doParallel)
library(coloc)

setDTthreads(threads = 20)

 

우선 필요한 패키지를 불어오고, 제 컴퓨터 코어는 20개라 thread 개수를 디폴트로 20으로 설정했습니다.

data.table로 fread 등 실행할 때의 설정입니다.

 

 

fread("~/driveD/R_projects/alzheimer/jansen1/AD_sumstats_Jansenetal_2019sept.txt.gz") -> AD_gwas

 

 

data.table이 가장 빠르죠. 명령어가 직관적이지 않지만, 속도가 중요해서 우선 불러올 때 fread로 불러옵니다.

 

AD_gwas %>% 
  rename(MAF = EAF) %>% 
  mutate(VARBETA = SE^2) -> AD_gwas

 

EAF라는 이름 대신 MAF로 사용

 

그리고 eQTL을 분석하는 coloc 패키지에서, SE 대신 VAR를 사용하게 됩니다. 그래서 VAR 값을 구합니다 (SE의 제곱)

 

 

(AD_gwas %>% 
  filter(P < 5 * 10^(-8)))$SNP -> target_SNPs

 

우선 분석의 대상으로 삼을 P < 5 x 10^(-8)인 SNP들만 선택합니다.

 

 

AD_gwas %>% 
  filter(SNP %in% target_SNPs) %>% 
  mutate(
    BP_start = BP - 1000000,
    BP_end = BP + 1000000
  ) -> target_SNPs2

 

보통 coloc으로 분석할 경우, 지역을 정해 분석하게 됩니다.

 

보통 1MB 위아래의 지역을 정하거나, 아니면 chromosome을 1개별로 정하거나 해서 분석할 지역을 설정합니다.

 

보통 1MB 정도 설정하는 것 같습니다. (cis니까...)

 

이 지역 설정하는 정도에 따라서, 보고자 하는 결과가 조금 달라집니다.

 

(왜그런지는 coloc 논문을 찾아보세요. 무슨 말인지는 대충 알겠는데, 설명 드릴 정도로

 

이해하지는 못했습니다. ㅜㅜ)

 

 

자. 이제. multicore를 돌려봅시다.

 

보통 그동안 R로 singlecore만 사용하셨을 겁니다.

 

하지만 GWAS에서 의미있게 나오는 SNP(즉, p < 5 x 10^(-8))이 2000여개 입니다.

 

거기에 위아래로 BP 1MB (즉 100만개)의 SNP들을 1 set으로 만들고, 그 각 1 set에 대해서

 

coloc을 돌려야 합니다. 시간도 매우 오래걸릴 뿐더러, 손으로 할 수 없습니다.

 

for 문을 쓰더라도, 역시 single core로 돌리기에는 시간이 오래 걸리기 때문에 multicore를 다 쓰는 방법을 이용해야 합니다.

 

우선 저는 doParallel 패키지를 사용합니다.

 

(클러스터서도 돌아가고 윈도우에서도 돌아가고 리눅스에서도 돌아간다고 합니다...)

 

 

 

 

n.cores <- 19

우선 코어 19개 쓸 거라고 설정해주고 (1개는 OS돌리거나 Rstudio 돌려야 하니까..)

slave.cluster <- makeCluster(n.cores)

core 갯수만큼의 cluster를 생성하고

registerDoParallel(slave.cluster)

cluster를 등록합니다.

 

clusterExport(slave.cluster,
              varlist = c("AD_gwas", "target_SNPs", "target_SNPs2") )

사용할 데이터 프레임을 미리 클러스터에 등록을 해 두어야 합니다

var list에 우리가 사용할 데이터프레임을 등록합니다.

 

 

 

foreach(i = seq_along(target_SNPs2[[1]]),
        .combine = rbind,
        .packages = c("tidyverse", "stringr", "data.table"),
        .inorder = TRUE) %dopar% {
          
          AD_gwas %>% 
            mutate(target = ifelse(CHR == as.vector(target_SNPs2[i , "CHR"]) & 
                                     BP > as.vector(target_SNPs2[i , "BP_start"]) & BP < as.vector(target_SNPs2[i , "BP_end"]),
                                   as.vector(target_SNPs2$SNP[i]), NA)) %>% filter(!is.na(target)) -> temp
          
          return(temp)
        } -> data_target

 

우선 seq_along으로 for 문에 들어갈 i를 설정해 줍니다

.combine은 쪼개서 만들어진 데이터를 rbind 한다는 뜻입니다. (row bind)

cbind, c가 있는데, c는 concatenate로, 데이터프레임으로 만들어지지 않습니다.

 

필요한 패키지는 미리 .packages로 등록시켜 주어야 합니다.

.inorder는 순서대로 한다는 이야기인데, 속도 생각하면 F로 설정해도 무방하지 싶습니다.

 

%dopar%는 parelle로 실행한다는 이야기고

 

%do%는 single로 실행하게 됩니다.

 

그 이후에 필요한 명령을 입력하시면 됩니다.

 

간단히 위에서 만들어 놓은 target_SNPs2라는 데이터프레임에 BP_start와 BP_end 범위 내에서 (BP 주변 1MB)

 

SNP들로만 coloc을 돌릴 수 있는 데이터 셋을 만드는 것입니다.

 

target은 그 데이터 셋들이 어떤 target SNP을 대상으로 하는지 말해주는 것이구요

 

 

하지만, 멀티코어로 돌려도 시간이 꽤 걸립니다.

 

램도 많이 필요합니다...