LD score regression 

LD score regression 은 Genome-wide association study(GWAS) 에서 특정 trait 의 polegenicity 를 추정하기 위해 사용하는 방법이다. LD score regression 은 GWAS summary statatistics 를 기반으로 SNP-heritability 추정, SNP-heritability 기반 genetic correlation 의 추정 등 다양한 measure 들을 계산하는데에 활용되고 있다. 본 문서에서는 LD score regression 의 등장 배경과 의미에 대해 알아보고자 한다. 

polygenic trait 에 대해 GWAS 를 수행한 후에, 각 SNP 들의 p-value 의 분포를 살펴보면 null distribution 과 비교하여 값들이 낮게 나타나는 것을 확인할 수 있다. 이렇게 높게 나타나는 검정 통계량은 형질이 polygenic 함을 의미할 수도 있지만, confounding bias 나 population stratification 가 영향을 주었을 수도 있다. polygenicity 로부터 위와 같은 bias 를 분리해내는 방법이 LD score regression 이며, LD score regression 은 이 과정에서 Linkage Disequillibrium (LD) 와 검정 통계량 (test statistics) 의 관계를 이용한다. 

"Both polygenicity (i.e. many small genetic effects) and confounding biases, such as cryptic relatedness and population stratification, can yield inflated distributions of test statistics in genome-wide association studies (GWAS). "

LD score regression 의 아이디어

어떤 SNP j 에 대해서 이 LD 관계에 있는 SNP 들이 많을 수록, polygenic 한 trait 에 대해서는 test statistics 이 높게 나올 가능성이 높다.  LD 는 유전자 변이간의 연관성을 의미한다. 만약, 어떤 SNP 이 LD 관계에 있는 SNP 이 많다고 하면, causal variant 와 LD 관계일 가능성이 높고, causal variant 와 LD 관계라면, test statistics 가 높게 나온다. 따라서 LD 관계에 있는 SNP 이 많으면, test statistics 가 높게 나올 가능성이 높다. LD score regression 은 이렇게 LD 관계에 있는 SNP 이 많은 SNP 일 수록 test stat 이 높게 나올 가능성이 높다라는 관계를 이용하는 방법이다. 그리고, 이러한 경향성이 강한 trait 일 수록 polygenicity 가 강하다고 말할 수 있다. (즉, LD score 와 test stat 의 연관성이 강할 수록 polygenic effect 로 phenotype 을 설명할 수 있는 비중이 높다.)

GWAS test statistics

LD score regression 에서는 chi-square value 를 regression 의 종속변수로 선정한다. chi-square value 는 무엇일까? 일반적으로 GWAS 결과로 effect size (beta) 와 standard deviation (sd) 값이 나오게 된다. 이 때, beta/sd 를 z-value 라고 한다. (이는 beta = 0 이라는 귀무가설 하에 구한 z-score 이다.) z-value 관측된 beta 값이 0으로부터 몇 standard deviation 떨어져 있는지를 의미한다. 이 때 chi-square value 는 z value 의 제곱으로 계산된다. z-value 는 표준정규분포를 따르며, 표준 정규 분포의 확률변수 z 의 제곱은 자유도가 1인 chi-square 분포를 따르기 때문이다. 

LD 와 LD score 는 무엇일까?

먼저 LD 와 LD score 를 계산하는 방법을 간단히 알아보자. 일반적으로 두개의 SNP A,B 의 LD 와 관련된 지표 D 와 r^2 은 아래와 같이 계산된다. D 값이 높을 수록 A,B 변이는 함께 나타날 가능성이 높음을 의미 한다. 만약, P(A)=0.3, P(B)=0.4, P(A,B) = 0.15 라고 하면, D는 0.03 으로 계산되며, A,B 는 LD 관계가 아닐 것으로 판단된다. 

$$D_{AB} = P(A \cap B) - P(A) P(B)$$  

일반적으로 많이 사용되는 지표인 r^2 은 아래와 같이 계산된다. 

$$r_{AB}^2 = \frac{D^2}{P(A)(1-P(A))P(B)(1-P(B))}$$  

특정 SNP j 에 대한 LD score 는 아래와 같이 계산된다. LD score 는 각각의 SNP 에 대해 다른 모든 SNP 들과의 LD 값 (r^2) 들을 더한 값으로 볼 수 있다. 

$$l_j = 1 + \sum_{k \neq j}r_{jk}^2$$  

LD score 와 test statistics (chi square value) 의 관계를 아래와 같이 시각화해볼 수 있다. 아래 차트는 LD score 의 bin 과 평균 chi-square value 의 관계를 보여준다. 직선은 아래 점들을 대상으로 단순 선형 회귀 분석을 한 결과를 표현한다 (아래 차트에서 각 점들에 해당 하는 SNP 의 갯수에 가중치를 두어 regression 을 돌리면 결국 전체 snp 과 chi-square 를 대상으로 regression 을 돌린 것과 같은 값이 나오게 될 것이다). 

이를 LD score regression 이라고 하며, 이 선형 회귀 분석에서 기울기는 polygenicity 를 반영하고, 절편은 bias 를 반영한다. LD 와 chi-square 의 연관성 (기울기) 이 polygenicity 이며, 전반적으로 chi-square 가 inflation 이 된 정도 (절편) 가 bias 라는 것이다. 이러한 방법을 통해 polygenicity 와 bias 를 분해할 수 있게 된다. 

또한, LD score regression 에서 기울기는 heritability 를 반영한다. 아래 그림과 같이 heritability 가 높을 수록 LD score regression 의 기울기가 커지게 된다. 

또한, 기울기는 샘플 사이즈와 사용한 SNP의 전체 개수에도 영향을 받는다. 샘플 사이즈 N 이 커질 수록 chi-square 값이 커지고, 사용한 SNP 의 개수가 많아질 수록, LD score 의 값이 기본적으로 높아진다. 이를 고려하여, 특정 SNP j 의 test stat 을 설명하는 아래와 같은 regression 모델을 고려할 수 있다.

$$E[\chi^2 | l_j] = Nh^2/M l_j + Na + 1$$

기울기는 heritability 와 N, M 으로 분해하여 나타낸다. 또한 절편은 Na + 1 로 표현되는데, 이 때, a 가 population structure 또는 confounding bias 와 같은 요인으로 인해 test stat 이 inflation 된 정도를 의미한다. 1이 더해진 이유는, LD score 가 0 인 SNP (즉, 그 어떤 SNP 과도 LD 관계에 있지 않은 SNP) 의 경우, chi-square value 는 causal variant 가 아닌 이상 기댓값은 1일 것이다 (자유도가 1인 chi-square distribution 의 평균값). 따라서, 절편은 1에 가까울 것이며, 1에서 벗어난 만큼을 bias 로 판단하겠다는 의미를 가진다. 

Bivariate LD score regression

두개의 trait X, Y 에 대한 test statistics 를 이용해 LD score regression 을 하는 것을 Bivariate LD score regression 이라고 한다. 구체적으로, 두개의 trait X, Y 에 대한 각각의 z-value 의 곱에 대하여 LD score regression 을 한다. 앞선 LD score regression 에 대해서는 chi-square value 를 사용했는데, chi-square value 는 z-value 의 곱이다. 즉, z^2 대신에 z_x * z_y 를 넣어서 regression 을 한다는 것이다.

참고) X_n 이 표준정규분포로부터 추출된 random variable 일 때, X_n 의 제곱의 합은 자유도가 n인 chi-square 분포를 따른다. 

위에서는 자유도가 1인 chi-square value 이기 때문에 z^2 = chi-square 가 된다. 

$$X^2_1 + X^2_2 + ... X^2_n \sim \chi^2(n)$$  

이 때의 기울기는 무슨 의미를 가질까? 만약 두개의 z-value 간에 아무런 연관성이 없다면, 1을 중심으로 퍼져있는 분포를 나타내게 된다. (z^2 의 분포는 자유도가 1인 카이제곱 분포를 따르며, 자유도가 1인 카이제곱 분포의 기댓값은 1이기 때문) 하지만, z-value 간에 연관성이 있다면, 기울기를 갖게 되며, 이 때의 기울기는 두 trait 간의 유전적 연관성을 의미한다. 두 trait 간에 유전적 연관성을 나타내는 지표로 co-heritability 가 있다. 기울기는 co-heritability 를 반영한다고 볼 수 있다. 이를 모델링하면 아래와 같이 표현할 수 있다. 

$$E[z_{xj}z_{yj}] = \frac{\sqrt{N_xN_y}h_{xy}^2}{M}l_j + \frac{\rho N_s}{\sqrt{N_xN_y}}$$  

왼쪽 그림에서 검은색 선은 유전적 연관성이 없는 trait 에 대한 기울기를 보여주며, 오른쪽 그림에서의 검은색 선은 연관성이 있는 trait 에 대한 기울기를 보여준다. 기울기가 가파를 수록 두 trait 간에 유전적 연관성이 높다고 볼 수 있다. 

참고자료

- LD Score Regression Distinguishes Confounding from Polygenicity in Genome-Wide Association Studies (2015)

- https://cnsgenomics.com/data/teaching/GNGWS22/module4/Lecture11_from_pgc_stat_bulik_2015.pdf

- An atlas of genetic correlations across human diseases and traits (2015)

- https://annahutch.github.io/PhD/LD-score-regression.html

genetic score, heritability, co-heritability, genetic correlation 관련 개념 정리

phenotype Y 의 분산은 genetics 로 설명되는 분산과 environment 로 설명되는 분산으로 나누어진다.  유전율 (heritability) 는 phenotype(또는 trait) 의 분산에서 genetics 를 통해 설명되는 분산을 의미한다. heritability 는 0~1 사이의 값을 가진다. 

$$Var(Y) = Var(G) + Var(E)$$

$$h^2 = Var(G)/Var(Y)$$  

genetics 로 설명되는 분산이란 무엇일까? phenotype 을 예측하기 위해 유전정보를 이용해 어떠한 score 를 만들고, 이를 genetic score 라고 하자. genetic score 는 phenotype 에 대한 예측 값이며, 이 값이 높을 수록 phenotype 의 값이 높을 가능성이 높음을 의미한다 (만약 질병과 같은 binary trait 인 경우, 질병의 걸릴 가능성이 높음을 의미한다.)

genetic score 는 유전체 정보를 이용해 구한 phenotype Y 에 대한 예측값이다. 따라서 아래와 같이 쓸 수 있는데 heritability 의 식이 결정계수의 식과 닮아 있음을 알 수 있다. 결정계수는 전체 분산중 어떠한 모델을 통하 예측값의 분산이 차지하는 비율이며, 이것이 곧, 모델을 통해 설명되는 분산을 의미한다. 

$$h^2 = Var(\hat{Y}) / Var(Y) = r^2$$  

genetics 로 설명되는 분산은 genetic score 의 분산으로 정의할 수 있다. 만약 genetic score 를 구할 때, additive genetic effect 만 고려하여, additive genetic score 를 구해 heritability 를 구한 것을 narrow-sense heritability 라고 한다.  

$$h_n^2 = Var(AG)/Var(Y)$$

만약 , Y 가 standardization 이 되어 있다고하면, Y의 평균은 0이고, Y의 분산은 1이다. 그러면, 간단히, additive genetic score 의 분산이 바로 narrow sense heritability 가 된다. 

"If the traits are standardized (that is, phenotypic variance = 1) and the genetic values consider only the additive genetic effects, then the genetic variances are narrow-sense heritabilities."

$$h_n^2 = Var(AG)$$

두 가지 trait 의 유전적인 연관성을 정량적으로 표현하는 지표로 coheritability 라는 개념이 있다.

"Co-heritability is an important concept that characterizes the genetic associations within pairs of quantitative traits."

co-heritability 는 아래와 같이 정의되며, -1~1 사이의 값을 가진다. 

$$h_{x,y} = \frac{Cov(g_x,g_y)}{Var(X)Var(Y)}$$

이 식의 의미를 살펴보면 분자의 covariance 에 Cov(X,Y) 가 오게 된다면, pearson 상관계수와 같음을 알 수 있다. 이 식은 Cov(X,Y) 대신에 X,Y 에 대한 genetic score 를 대입시킴으로써, 두 trait 의 유전적 상관성을 표현했다고 볼 수 있다. 여기서도 마찬가지로 trait X,Y 를 평균이 0이고 분산이 1인 표준화된 trait 을 사용했다면, Var(X) = Var(Y) = 1 이기 때문에 아래와 같다. 

$$h_{x,y} = Cov(g_x,g_y)$$

두 가지 trait 의 유전적인 연관성을 정량적으로 표현하는 지표로 genetic correlation 이라는 개념도 있다. genetic correlation 은 아래와 같이 정의된다. 

"The genetic correlation is a quantitative genetic parameter that describes the genetic relationship between two traits"

$$\rho_g = \frac{Cov(g_x,g_y)}{\sqrt{Var(g_x)Var(g_y)}}$$  

위 식은 pearson 상관계수의 식과 같으며, genetic correlation 의 통계적인 의미는 X,Y 의 genetic score 상관성 (pearson 상관계수) 라고도 할 수 있다. 만약 두가지 trait, 예를 들어 키와 발가락 길이의 유전적 연관성이 높다라고 한다면, 유전자를 통해 예측한 키 (키에 대한 genetic score) 와, 예측된 발가락 길이 (발가락 길이의 genetic score) 의 연관성이 높을 것이다. 이를 수치화한 것이 genetic correlation 이라고 볼 수 있다. genetic correlation 도 마찬가지로 -1~1사이의 값을 가진다. 

genetic correlation 과 co-heritability 모두, 두가지 trait 의 유전적 연관성을 표현한다. 둘의 차이점은 무엇일까?  trait X,Y 가 표준화 되어있다고 하면 genetic correlation 은 아래와 같이 정의된다. 아래 식을 보면, genetic correlation 은 co-heritability 가 X,Y 각각의 trait 의 heritability 로 보정된 식임을 알 수 있다. 

$$\rho_g = \frac{h_{x,y}}{\sqrt{h^2_x h^2_y}}$$  

따라서, 두 trait 의 heritabilty 값이 작더라도, genetic correlation 은 높을 수 있다. 예를 들어, 발가락 길이와 키의 heritability 가 10% 라고 하자 (실제로는 더 높을 것이나 예시임). 즉, 전체 분산에서 genetic score 의 분산이 차지하는 부분이 10% 이다. 하지만, 두개의 genetic score 의 연관성이 높다라고 하면, genetic correlation 은 높게 추정될 수 있다. 따라서, genetic correlation 을 해석할 때, trait 을 genetics 가 설명하는 비중 (heritability) 도 함께 고려해야할 필요가 있다. 

참고자료

- Genetic correlations of polygenic disease traits: from theory to practice, Nature review genetics, 2020

- Optimal Estimation of Co-heritability in High-dimensional Linear Models

- Statistical methods for SNP heritability estimation and partition: A review

List 를 이용한 방법

# dictionary 생성
dict <- list(name = "John", age = 30, city = "New York")

# dictionary 사용
dict$name
# [1] "John"

dict$age
# [1] 30

dict$city
# [1] "New York"

vector 를 활용한 방법

ㄴ setNames 함수를 활용

# dictionary 생성
dict <- setNames(c("John", 30, "New York"), c("name", "age", "city"))

# dictionary 사용
dict$name
# [1] "John"

dict$age
# [1] 30

dict$city
# [1] "New York"

hash 함수를 활용한 방법

library(hash)
h <- hash() 
h[['a']] <- 'a'
h[['b']] <- 'b'
h[['c']] <- 'c'
h[['d']] <- 'd'

h[['a']]

R 에서 특정 변수를 카테고리화 하고 싶을 때가 많다. 

다양한 방법이 있지만, 

아래 cut 함수를 사용하는 코드로 0~5, 6~10, 11~15, ... >100 으로 카테고리화가 가능하다. 

cat <- seq(0,100,5)
df$cat <- cut(df$x, breaks = c(cat, Inf), labels = cat)
df$cat <- factor(df$cat, levels=cat)

-> breaks 의 element 보다 labels 의 elements 의 갯수가 1개 적다. 

좀 더 일반적으로는 다음과 같다.

# 예시 데이터 생성
set.seed(123)
data <- data.frame(id = 1:10, value = rnorm(10, mean = 50, sd = 10))

# 카테고리화
data$cat <- cut(data$value, breaks = c(0, 25, 50, 75, 100), labels = c("low", "medium-low", "medium-high", "high"))

mutate_at 을 이용해서 아래와 같이 일괄적으로 cut 함수를 여러 컬럼에 적용할 수 있다. 

data <- data %>%
  mutate_at(vars(a, b, c), 
            list(cat = ~cut(., 
                            breaks = c(0, 25, 50, 75, 100),
                            labels = c("0", "1", "2", "3"),
                            include.lowest = TRUE)))

Hmisc 의 Lag 변수를 통해 timeseries 데이터의 lag 변수를 만들 수 있다. 

만약, 그룹별 Lag 변수를 만들고 싶으면 dplyr group_by 를 통해 만들 수 있다. 

library(Hmisc)
data <- data %>% group_by(gender, age) %>% mutate(lag = Lag(variable, 1))

반복문을 돌면서 여러개의 dataframe 을 만들고, 

이것들을 합쳐서 최종적인 결과 dataframe 을 만드는 경우 아래와 같이 함

brand_names <- c("BBQ", "BHC")
tmp_dfs <- lapply(brand_names, function(x){
  p <- paste0(x, "_payment")
  b <- paste0(x, "_buzz")
  tmp_df[1, 'brand'] <- x
  tmp_df[1, 'cor'] <- cor.test(merged[,p], merged[,b])$estimate
  tmp_df
})

df_result <- bind_rows(tmp_dfs)
df_result

출력
> df_result
        cor brand
1 0.6971354   BBQ
2 0.4675438   BHC

모든 na값을 0로 바꾸기

df[is.na(df)] <- 0

특정 컬럼의 na 값을 0으로 바꾸기

df[is.na(df$col_name), 'col_name'] <- 0

여러 컬럼의 na 값을 0으로 바꾸기

impute_var <- c('a','b','c')
df[impute_var][is.na(df[impute_var])] <- 0

Pyspark export to delimited file

def myConcat(*cols):
    concat_columns = []
    for c in cols[:-1]:
        concat_columns.append(F.coalesce(c, F.lit("*")))
        concat_columns.append(F.lit("\t"))  
    concat_columns.append(F.coalesce(cols[-1], F.lit("*")))
    return F.concat(*concat_columns)

# combined column 에 모든 변수를 \t 로 concat 한 값 저장 
data_text = data.withColumn("combined", myConcat(*data.columns)).select("combined")
data_text.coalesce(1).write.format("text").option("header", "false").mode("overwrite").save(path)

출처 : https://stackoverflow.com/questions/17837871/how-to-copy-file-from-hdfs-to-the-local-file-system

가끔 .gitignore 파일에 올리고 싶지 않은 파일이나 폴더를 지정후, remote repo 에 push 를 해도 적용되지 않을 때가 있다. 

이 때는 캐시를 지우고 시도한다. 

git rm --cached -r .
git add .

Jupyter notebook에서 python 파일을 자동으로 갱신하는 방법

예를 들어 jupyter notebook 파일에서 model.py 라는 파일을 import 해서 사용하고 있다고 하자. 

만약 model.py 파일을 수정하면, 커널을 재시작한 후, 다시 import 를 해야한다. 

아래 구문을 통해 model.py 파일이 수정되는 즉시 적용되도록 만들 수 있다. 

%load_ext autoreload
%autoreload 2