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

이제 다운로드 받은 eQTL 데이터를 불러옵니다

 

AD gwas 데이터를 사용했으니

 

eQTL 데이터는 brain을 사용해 보도록 하겠습니다.

 

 

dir <- ("~/driveD/ref_gene/GTEx_v8/GTEx_Analysis_v8_eQTL/")

우선 이렇게 하시면 eQTL 데이터가 있는 폴더를 정해주시고

file_list <- list.files(dir)

이러면 이 안에 그 폴더 내 파일명들이 모두 저장이 됩니다.

file_list[str_detect(file_list, "Brain") & str_detect(file_list, "signif") ]  -> file_list

그 파일명 중 Brain이라는 단어와 signi라는 단어가 있는 것만 선택합니다.

왜냐하면 egene이 있고 significant 가 있는데, significant가 필요합니다.

(훨씬 양이 많습니다. egene은, 말 그대로 gene만 저장해놓은 파일 같습니다.)

for (file in file_list) {
  
  assign(file, fread(paste(dir, file, sep = "/"), nThread = 20) )
  
}

자 이제, for 문을 이용해서 파일들을 연속해서 불러옵니다

 

paste를 이용해서 파일 주소를 만들고

fread에 인자로 넣어서 불러오고

그것을 assign으로 file_list의 파일 명으로 저장을 하도록 만듭니다.

 

자동화 만세..

 

 

근데 AD gwas 데이터는 hg19 입니다.

eQTL은 hg38이구요.

 

SNP (즉, rs ID)만으로는 loss되는 경우가 많아서, (같은 SNP인데 이름 다른 경우가 종종 있지요..)

 

그래서 저는 BP (base position)으로 매칭을 시킵니다.

우선 GTEx 포털에, 아래쪽으로 내리다보면 파일을 찾을 수 있습니다.

(Reference에 존재합니다.)

hg38을 19로 변경할 수 있도록 만든 lookup table 입니다.

 

wget https://storage.googleapis.com/gtex_analysis_v8/reference/GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.lookup_table.txt.gz

역시 리눅스에서 내려받습니다.

 

 

fread("~/driveD/ref_gene/GTEx_v8/GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.lookup_table.txt") -> hg38_to_hg19

그 테이블을 불러옵니다.

setkey(hg38_to_hg19, variant_id)

그리고 dplyr도 느리기 때문에, data.table을 활용할 계획입니다.

data.table 패키지의 merge 기능을 사용하기 위해서는 미리

setkey 설정을 해놓는 것이 필요합니다.

우리는 variant_id (1_123456_T_A 식으로 되어있습니다.)를 이용해서 merge할 것이기 때문에

setkey를 variand_id로 할 겁니다.

 

 

이제 eQTL 파일 불러온 것들을 hg19로 변경할 겁니다

 

ls()하면 모든 데이터들의 이름이 나타나고

ls()[str_detect(ls(), "Brain_")]

 

이렇게 하면 그 중 Brain_ 이라는 단어가 들어있는 것을 추려내지요.

 

ls()[str_detect(ls(), "Brain_")] -> merge_list

 

merge_list를 만들어서 역시 for 문을 이용해서 merge를 할 겁니다

 

for (file in merge_list) {
  
  str_sub(file, 7, nchar(file)) -> temp
  (str_split(temp, ".v8", 2) %>% unlist())[[1]] -> temp
  paste("eQTL",temp,sep="_") -> temp
  print(temp)
  assign(temp, merge(hg38_to_hg19[ , list(variant_id, ref, alt, rs_id_dbSNP151_GRCh38p7, variant_id_b37)], get(file), by = "variant_id", all=F)  )
  
}

 

이제 eQT_Amygdala에서부터 eQTL_Substantia_nigra 파일들이 생성되었습니다.

 

모두 hg38 -> hg19로 변경된 파일들입니다.

다만 약간의 최종 정리는 좀 필요해서 수정을 가합니다.

왜냐하면 GWAS 파일과 merge가 필요하기 때문입니다.

 

ls()[str_detect(ls(), "eQTL_")] -> merge_list

이제는 eQTL_로 시작하는 파일 명단을 만들고

 

 

for (file in merge_list) {
  
  get(file) %>% 
    mutate(
      BP = (str_split(variant_id_b37, "_", n=3, simplify = T))[ ,2],
      CHR = (str_split(variant_id_b37, "_", n=3, simplify = T))[ ,1]
    ) %>% 
    filter(
      CHR %in% as.character(seq(1,22,1))
    ) %>% 
    mutate(
      CHR = str_sub(CHR, 1, nchar(CHR)),
      VARBETA = (slope_se)^2
    ) %>% 
    mutate(
      BP = as.integer(BP),
      CHR = as.integer(CHR)
    ) %>% 
    rename(
      SNP = rs_id_dbSNP151_GRCh38p7,
      P = pval_nominal,
      A1 = alt,
      A2 = ref,
      BETA = slope,
      SE = slope_se,
      MAF = maf
    ) %>% 
    select(CHR, BP, SNP, BETA, SE, VARBETA, P, MAF, A1, A2, gene_id, variant_id_b37, variant_id, gene_id, everything()) %>% 
    arrange(CHR, BP) -> temp
  assign(file, temp)

}

 

향후 과정을 위한 간단한 설정을 해줍니다.

 

 

이제 한 번 저장합시다.

 

아마 컴퓨터가 돌아간 시간만 1~2시간이 넘게 걸렸을 겁니다.

 

rm(hg38_to_hg19)

우선 파일 용량이 매우 큰 lookup table을 지웁시다.

 

save.image("AD_gwas_eqtl.RData")

그리고 저장하였습니다.