Alzheimer's disease, GWAS & eQTL 분석 (4) coloc 분석

coloc은 colocalization을 분석하는 패키지 입니다.

 

자세한 사항은 관련 논문 참조해 주세요

 

쉽게 설명하면 GWAS의서의 summary statistics, eQTL에서의 summary statistics를 이용해서

 

abf를 계산하는 것이라고 합니다.

 

 

 

자 이제 coloc 패키지의 coloc.abf를 이용해서 자동으로 분석해 보도록 합시다.

 

자동화를 위하여 함수를 만들어야 합니다.

 

 

automated_coloc <- function(target_gwas, eqtl, target_SNP, eqtl_sample_no) {
  
  inner_join(
    eqtl,
    target_gwas %>% filter(target == target_SNP), 
    by = c("CHR","BP"), suffix = c("_eqtl", "_gwas")
    ) -> merged_data
  
  merged_data %>% 
    filter(!is.na(MAF_gwas)) -> merged_data
  
  
  if (NROW(merged_data) == 0) {
    return(NA)
  } else{
  
  coloc.abf(dataset1 = list(beta = merged_data$BETA_gwas, varbeta = merged_data$VARBETA_gwas,
                            N = merged_data$Nsum, MAF = merged_data$MAF_gwas,type = "quant"), 
            dataset2 = list(beta = merged_data$BETA_eqtl, varbeta = merged_data$VARBETA_eqtl,
                            N = eqtl_sample_no, MAF = merged_data$MAF_eqtl, type = "quant")) -> result_coloc
  
  result_coloc$results -> result_coloc
  
  result_coloc %>% 
    mutate(
      no = as.numeric(str_sub(snp, 5, nchar(snp)))
    ) %>% select(no, everything()) %>% 
    arrange(no) -> result_coloc
  
  cbind(result_coloc, merged_data) %>% 
    select(SNP_eqtl, SNP.PP.H4, everything(), -no, -snp) %>% data.table() -> final_result_coloc
    #filter(SNP ==target) %>% data.table() -> final_result_coloc
    
  return(final_result_coloc)
  }
}

 

 

함수의 변수는 gwas summary statistics, eqtl summary, target으로 하는 SNP, eqtl의 샘플 갯수 입니다.

 

우선 맨 처음 만들었던 gwas의 target SNP (p < 5 x 10^(-8) 이었던 2000여개의 SNP)마다 위아래로 100MB의 SNP

 

들을 추려서 데이터셋을 만들었습니다.

 

이제 그 각 target SNP 별로 분석을 시행해야 합니다. 전체 gene set으로 분석하면 안된다고 합니다. cis eQTL이니

 

어느 정도 지역을 정해두어야 하겠죠. 그게 보통 1MB라고 합니다.

 

우선 eQTL과 1MB 위아래 정도의 GWAS만 머지합니다.

 

분석에 MAF가 필요한데 (없으면 에러), MAF가 없는 데이터들을 지워줍니다.

 

sdy라는 값을 구하기 위해서 MAF와 n수가 필요합니다. (coloc 패키지 참조)

 

우선 coloc.abf를 이용해서 gwas 데이터와 eqtl 데이터를 list 형식으로 넘깁니다.

 

coloc.abf 위해 2개 summary data set이 필요하고, beta, var, N, MAF, type이 최소로 필요하다고 합니다.

 

N, MAF는 sdy를 계산하기 위함이고, type은 cc (case control)인지 quant인지 (연속 형질) 구분하기 위해 필요합니다.

 

원본 GWAS 파일에 beta로 되어있는 것 보니 (OR가 아니라), 연속 형질인가 봅니다.

 

 

data_target %>%
  mutate(target_id_b37 = paste(CHR, target_BP, sep = "_") ) %>%
  filter(
    target_id_b37 %in%  (eQTL_Amygdala %>%
                                mutate(variant_id_b37_2 = paste(CHR, BP, sep = "_") ))$variant_id_b37_2
  ) -> data_target_final

 

그리고 우선 실행하기 전에 마지막으로 eQTL데이터가 있는 것들만 추려냅니다.

 

 

for (file in merge_list) {
  
  get(file) -> temp1
  
  paste("data_target",  str_split(file, "_", n=2, simplify = T)[[2]], sep = "_") -> new_file_name
  
  data_target %>%
    mutate(target_id_b37 = paste(CHR, target_BP, sep = "_") ) %>%
    filter(
      target_id_b37 %in%  (temp1 %>%
                             mutate(variant_id_b37_2 = paste(CHR, BP, sep = "_") ))$variant_id_b37_2
    ) -> temp2
  assign(new_file_name, temp2)

}

 

역시 for 문을 이용해서 자동화...

 

간단하게 data_target이라고 했던 위에서 만든 큰 파일 중에서

target SNP이 각 eQTL데이터 셋에 있는지 검색해서,

eQTL에 있는 것만 하나씩 저장하는 것입니다.

 

 

이제 마지막으로 실제 분석을 하도록 하겠습니다.

 

여기서도 굉장히 오래 걸리기 때문에,

 

multicore를 이용하도록 하겠습니다.

 

 

n.cores <- 19

코어 갯수를 세팅하고

slave.cluster <- makeCluster(n.cores)

클러스터로 등록하고

registerDoParallel(slave.cluster)

클러스터를 역시 패러럴로 하기 위해 등록

clusterExport(slave.cluster,
               varlist = c("data_target_Amygdala", "target_SNPs", "eQTL_Amygdala") )

변수를 등록합니다.

 

아무래도 여기에서 부터는 하나씩 시행하는 것이 좋습니다.

용량과 연산이 너무 커집니다 -_-

변수명만 신경쓰세요

 

data_target_Amygdala 였으니, eQTL_Amygdala를 써야합니다.

그리고 Amygdala의 eQTL sample 수는 129개라고 합니다.

이 역시 GTEx portal에 있습니다.

(genotype과 같이 있는 샘플 수를 사용해야할 것 같습니다.

GTEx map 만들 때 당연히 genotype도 같이 사용하기 때문에..)

 

foreach(i = seq_along(target_SNPs),
        .combine = rbind,
        .packages = c("tidyverse", "stringr", "data.table", "coloc"),
        .inorder = TRUE) %dopar% {

          automated_coloc(data_target_Amygdala, eQTL_Amygdala, target_SNPs[[i]], 129) -> temp
          
          if( all(is.na(temp)) ) return(NULL)
          
          return(temp)
          
        } -> result_eQTL_Amygdala

stopCluster(slave.cluster)

cluster를 멈춰주고

 

 

 

 

 

자 이제 automated_coloc이라는 함수를 연속적으로 시행하게 됩니다

19개의 코어(노예)가 나누어 시행을 하게 되고

결측값이 있으면 NULL를 되돌려 줍니다 (return(NULL))

(그래야 에러가 안나더군요).

 

 

얻어진 결과에서

 

 

 

그리고 결과를 fwrite로 저장해 줍니다.