유전체 연구와 자료의 대표성

이 포스트는 제 연구를 소개하는 글입니다. 이 포스트를 시작으로 진화생물학 및 통계유전학이 사회와 어떤 관계를 맺고 있는지 소개하는 글을 몇 편 쓰려고 합니다. 대중적으로 유명한 진화심리학 등이 과학적 발견의 사회적 의의에 초점을 맞췄다면 저는 반대로 사회적 구조가 유전학 연구에 어떤 영향을 미치는지 다루려고 합니다.

저는 서울대학교 의과대학을 다녔습니다. 서울대학교 의과대학에는 크게 3개의 수련 병원이 있습니다. 연건동에 위치한 서울대학교병원 (소위 본원이라고 불립니다), 분당서울대학교병원 그리고 동작구에 위치한 시립보라매병원입니다.

서울대학교병원 전경

세 병원은 각각 성격이 다릅니다. 본원이 상급종합병원으로서 암이나 희귀질환 같은 난치성 질환 치료와 연구에 집중한다면 보라매병원은 시립병원으로서 취약계층 진료에 그 목적을 두고 있습니다. 이런 차이는 학생으로 실습기간 동안 충분히 느낄 수 있었습니다. 점심시간이 넉넉한 날이면 현관 근처에서 커피를 마시며 멍때리곤 했는데 과장 조금 보태면 10분에 한 번씩 양복을 입은 기사님이 운전하는 S클래스에서 하차하는 어르신을 봅니다. 온 동네 고급차량은 병원에서 다 본 것 같습니다. 한편, 보라매병원은 완전히 다른 세상입니다. 건강보험료를 지급하지 못하여 공공부조의 일종인 의료급여를 받는 환자가 굉장히 많습니다. 외래에서 뵙는 환자 중에 병이 상당히 진행된 상태에서 내원하는 비수도권 환자도 적지 않게 보입니다.

동시에 요즘은 빅데이터와 인공지능이 각광을 받는 시대입니다. 이런 흐름에 발맞춰 학교와 병원에서도 십수년 간 축적된 환자 데이터를 바탕으로 새로운 지식을 창출하려는 연구를 활발하게 하고 있습니다. 그런데 저는 그런 생각이 들었습니다. 50억 정도 하는 압구정 현대 아파트에 사는 어르신과 경남 양산 시골에 사는 어르신이 서울대병원에 내원할 기회가 동등하지 않은데 이런 편향된 환자의 모임에서 생산된 데이터로 좋은 연구를 할 수 있을까?

인과추론 (Causal inference)는 관찰된 상관관계에서 인과성을 창출하는 방법을 연구하는 분야로 제 의문 역시 인과추론에서 오랫동안 연구된 문제입니다. 이처럼 대표성 없는 자료에서 인과성에서 비롯되지 않는 무의미한 상관관계가 생기는 현상을 선택편향 (selection bias)라고 부릅니다. 선택편향은 우리 생활 곳곳에 퍼져있습니다. 흔히, 많은 기혼 여성이 만났던 사람 중에 현재 배우자가 성격은 좋지만 외모는 제일 못났다고 한탄하곤 합니다. 외모와 인품의 반비례 관계는 어떤 인과성에 의한 것일까요? 선택편향에 대한 이론은 여기에 대한 확고한 답을 줍니다. 외모와 인품 모두 교제의 가능성을 높이기 때문에 (교제가 외모와 인품의 결과이기 때문에) 교제했던 사람만 모아서 보면 둘이 반비례하는 것처럼 보이는 것입니다. 교제가 외모와 인품 모두와 완전히 독립적으로 일어났다면 이런 상관관계는 생기지 않았을 것입니다.

위 같은 이유로 서울대병원의 빅데이터는 여러 요인의 결과로 서울대병원에 내원할 수 있었던 환자만 분석 대상으로 삼기 때문에 앞선 예시와 같은 문제를 겪을 수 밖에 없습니다. 이런 의료 빅데이터의 한 종류로 대규모 유전체 연구가 있습니다. 저는 유전학을 공부하고 이런 자료를 분석하고 있습니다. 그런데 유전체 연구도 연구에 자발적으로 참여하는 사람만 포함하고 있으므로 선택편향으로부터 자유로울 수 없습니다. 그렇다면 그 선택편향의 크기는 얼마고 어떤 요인의 영향을 받는지 밝힐 필요가 있고 제 연구는 정확히 이 문제를 다룹니다.

2021년 네이쳐 유전학 (Nature Genetics)에 출판된 연구는 생물학적 성별과 상염색체 유전 변이 사이의 상관관계를 분석했습니다. 생물학적 성별은 성염색체가 결정하고 상염색체와 무관하므로 아무런 상관관계가 나타나지 않을 것으로 기대할 수 있으나 결과는 충격적이게도 그렇지 않았습니다. 놀랍게도 수백 개의 상염색체 변이가 성별과 상관관계가 있는 것으로 드러났습니다.

Pirastu 외가 밝혀낸 성별-상염색체 변이 상관관계. 빨간 점선 위에 있는 점이 모두 성별과 관련성이 있는 상염색체 변이이다.

이 연구는 선택편향를 이러한 상관관계의 원인으로 지목했으나 정확히 그 수학적인 기전을 밝히진 못했습니다. 제 연구는 그들이 밝히지 못한 수학적인 기전과 일반적인 이론을 제안하는 것입니다.

이한빈 외.는 유전체 연구에서 발생하는 선택편향에 대한 일반적 이론을 제안합니다

왜 이런 연구를 했는지 생각해보면 의과대학에서의 경험이 결정적이었습니다. 저는 대학에 오기까지 20년을 부산에서 살았는데 부산은 제2 도시로서의 위상에도 불구하고 전라권이나 경북권에 비해 의료 환경이 매우 낙후됐습니다. 때문에 서울에서 경험한 의료환경은 그 자체로 놀라운 경험이었고 그 안에서도 본원과 보라매 병원의 차이와 같은 불평등의 결과가 항상 마음 한 켠에 있었습니다. 그런 와중에 네이쳐에 출판된 연구를 접하면서 제 문제의식을 구체화할 수 있는 기회라고 생각하여 금번의 연구를 시작할 수 있었습니다.

비록 논문 본문에 적지 않았지만 제 연구는 몇 가지 중요한 점을 시사하고 있습니다. 하나, 대규모 유전체 연구에서 선택편향을 제거하기 위해서는 각 피연구자가 연구에 참여할 확률을 알아야 합니다. 그리고 그 확률을 알기 위해서는 어떤 경위로 그 사람이 연구에 참여하게 되는지 반드시 이해할 필요가 있습니다. 둘, 이런 사실을 한국의 맥락으로 가져오면 한국은 지역 간 의료불평등과 계층에 따른 의료접근성의 차이가 대단히 큰 나라입니다. 그런 나라에서 아무리 많은 데이터를 바탕으로 연구를 한다고 한들 선택편향으로부터 자유로울 수 없습니다. 때문에 엄밀한 연구를 위해서라도 건강불평등의 원인과 해결책을 고민할 수 밖에 없는 것입니다. 마지막으로, 사회적으로 의료접근성이 균등한 곳일수록 최첨단의 연구가 더 쉽고 간편하다는 것입니다. 그런 사회에서는 의료접근성으로 인한 선택편향이 없을테니까요.

저는 집단유전학과 인과추론의 교집합에 위치한 사람입니다. 그 자체로 의료접근성이나 불평등을 연구하진 않는다는 뜻입니다. 그러나 제가 다루는 자료와 데이터는 사회에서 만들어지고 바로 그렇기 때문에 의료접근성과 불평등을 무시할 수가 없습니다. 금번의 연구는 그 점을 수학적인 관점에서 간접적으로 보여주고 있습니다.

Marker SNP as a noisy measurement

Genome-wide association study (GWAS) runs millions of regressions across the genome to identify genomic regions that are responsible for the variatiblity in a trait. The idea stems from the concept of linkage disequilibrium (LD), a phenomenon where two genomic regions are statistically correlated due to physical linkage. Technically, the correlation may come from population phenomenons but the post will consider only physical linkage for simplicity. There are approximately 3 billion base pairs in the human genome (even larger for plant genomes) but the actual number of regressions are usually around several million. As long as the genotyped markers are located dense enough, they will cover the whole genome through LD. Even if the true causal variant is not genotyped, it is likely that it’s effect will be detected in a regression of a marker nearby. Although REAL whole genome regressions based on WGS data are released nowadays, this rationale makes genotype array based GWAS a cost-effective method for identifying trait associated loci.

More …

Instrumental variables in econometrics and epidemiology

Mendelian randomization (MR) became very famous in the recent years. Even clinical journals that are quite conservative in using causal language outside randomized trials also embrace MR studies nowadays. In the methodological field, many new MR methods have been proposed. However, I’m very skeptical about them as they impose extremely strict parametric assumptions. For example, with very few execptions (see Shi 2021 and 2022), they assume strict homogeneity (see Morrison 2020). Furthermore, the biggest problem in current literature is that we are not very transparent about the parameter we are estimating. The target parameter is absent!

In this post, I review the two approaches in instrumental variable. The first one is the local average treatment effect (LATE) framework famous in econometrics (Imbens and Angrist, 1994). The second ons is the structural mean model (SMM) framework proposed by the Harvard CI group (Robins, 1994).

We start from the switching formula.

\[Y_i = D_i \cdot Y_i(1) + (1-D_i) \cdot Y_i(0)\]

which gives

\[Y_i = \tau_i \cdot D_i + Y_i(0)\]

after reordering. $\tau_i = Y_i(1) - Y_i(0)$ is the individual treatment effect.

Learning the importance of heterogeneous TE in IV was a very illuminating experience. At the same time, the plethora of many assumptions was very confusing. How does these different assumptions relate to each other? I don’t have a full-blown answer to this question. However, I think that the two approaches by econometricians and epidemiologists can be classified in terms of the stage in which the assumption is applied. The former, which is usually referred to monotonicity is a matter of how the instrument $Z$ affects the instrument $D$. On the other hand, the latter (e.g. no treatment effect modification (NEM) assumption) imposes restriction on the mode of treatment $D$ influencing the outcome $Y$. The following derivations will make this point clear.

Monotonicity

The proof is essentially the same as in Mostly Harmless Econometrics (Theorem 4.4.1, p.155). Apply $E[\cdot \vert Z]$ to equation (2).

\[E[Y_i \vert Z_i] = E[\tau_i \cdot D_i \vert Z_i] + E[Y_i(0) \vert Z_i]\]

The last term is just a constant $E[Y_i(0)]$ by the exclusion criteria $Y(d) \perp\kern-5pt\perp Z$.

\[E[Y_i \vert Z_i=1] - E[Y_i \vert Z_i=0] \\ = E[\tau_i \cdot D_i \vert Z_i=1] - E[\tau_i \cdot D_i \vert Z_i=0]\]

To simplify this equation, we have two choices. One is to impose some condition on $\tau_i$ and the other is to impose some condition on $D_i$. Monotonicity does the latter by assuming no defiers. Let $C_i = D_i(1) - D_i(0)$ be the compliance status. Then

\[E[\tau_i \cdot D_i \vert Z_i] \\= E[ E[ \tau_i \cdot D_i \vert Z_i, C_i ] \vert Z_i] \\= \sum_c E[\tau_i \cdot D_i \vert Z_i, C_i=c] \cdot P(C_i =c \vert Z_i)\]

By exogeneity of the instrument, $P(C_i \vert Z_i) = P(C_i)$. For $Z=1$,

\[= \sum_c E[\tau_i \cdot D_i \vert Z_i=1, C_i=c] \cdot P(C_i =c) \\ = E[\tau_i \cdot 1 \vert Z_i=1, C_i=1] \cdot P(C_i =1) \\ + E[\tau_i \cdot D_i \vert Z_i=1, C_i=0] \cdot P(C_i=0)\]

and for $Z=0$,

\[= E[\tau_i \cdot 0 \vert Z_i=0, C_i=1] \cdot P(C_i =1) \\ + E[\tau_i \cdot D_i \vert Z_i=0, C_i=0] \cdot P(C_i=0)\]

Finally, subsituting (6) and (7) into (4) and the exclusion criteria guarantees

\[E[\tau_i \cdot D_i \vert Z_i=1] - E[\tau_i \cdot D_i \vert Z_i=0] \\ = E[\tau_i \cdot 1 \vert C_i =1] \cdot P(C_i=1) \\ = E[\tau_i \vert C_i=1] \cdot E[D_i(1) - D_i(0)] \\ = E[\tau_i \vert C_i=1] \cdot (E[D_i \vert Z=1] - E[D_i \vert Z=0])\]

where the last line came from exogeneity. The proof shows that there is no restriction on $\tau_i$ and only the assumptions involving $D_i$ was used to derive the Wald ratio.

Structural Mean Model (SMM)

The additive SMM is

\[E[Y-Y(0) \vert D, Z] = (\psi_0 + \psi_1 Z) \cdot D\]

When I saw this for the first time, I was confused. Most importantly, the interpretation of $\psi_0$ and $\psi_1$ wasn’t very transparent to me. Therefore, I thought it was better to start from equation (2) to make it more sensible.

Applying $E[\cdot \vert D,Z]$ to equation (2) gives

\[E[Y_i - Y_i(0) \vert D_i,Z_i ] = E[\tau_i \cdot D_i \vert D_i, Z_i] \\ = E[\tau_i \vert D_i, Z_i] \cdot D_i \\ = E[\tau_i \vert Z_i] \cdot D_i\]

So, what was happening in the SMM was the specification of $E[\tau_i \vert Z_i]$ which is

\[E[\tau_i \vert Z_i] = \psi_0 + \psi_1 Z_i\]

NEM imposes $\psi_1 = 0$ so that $\psi_0$ eventually becomes the PATE. Embracing the soul of Wooldridge’s idea that I wrote, using $Z - \mu_Z$ instead of $Z$ would have given $\psi_0$ the same interpretation without the NEM although the problem of identification would have remained (the number of moment condition is smaller than the number of estimands).

To obtain the Wald ratio, apply the law of iterated expectation on (assuming NEM),

\[E[Y_i - Y_i(0) \vert Z] = \\ E[ E[Y_i - Y_i(0) \vert D,Z] \vert Z ] \\ = E[ \psi_0 \cdot D \vert Z ] \\ = \psi_0 \cdot E[ D \vert Z ]\]

and substituting it into equation (4) gives the desired result. Examining the consequence of NEM on equation (11) directly shows that SMM-based approaches acheives identification by restricting the mode of $D \rightarrow Y$.

Concluding remark

My impression of these result is that monotonicity and SMM-based assumptions are not necessarily stronger/weaker than the other. The choice ultimately depends on which stage, $Z \rightarrow D$ or $D \rightarrow Y$, will be restricted through imposing additional assumptions. May be some person can come up with an alternative identification strategy by imposing restrictions on both but each of them being weaker for each stage.

Regress-out in single-cell biology

I first intended this content as a research paper but I’m currently not planning any additional single-cell projects so I write this as a blog post. In many, I mean really many single-cell analysis notes, people suggest regressing out certain variables prior to analysis. For example, in the cell-cycle vignette in Seurat manpage, they recommend regressing out cell-cycle scores prior to PCA. I think such practice is not legitimate because it wipes out not only the variation due to cell-cycle but also other biological variations that are correlated to cell-cycle. I show, in this post, that this is ture in general.

More …