안녕하세요, 이번 글에서는 오픈 엑세스 논문을 하나 정리해보고자 합니다. 논문 제목은 'Single-cell eQTL mapping reveals cell-type-specific genes associated with the risk of gastric cancer'로, 저번 달에 cell genomics 저널에 개재된 논문입니다.
Single-cell eQTL mapping reveals cell-type-specific genes associated with the risk of gastric cancer - PMC
Summary Most expression quantitative trait locus (eQTL) analyses have been conducted in heterogeneous gastric tissues, limiting understanding of cell-type-specific regulatory mechanisms. Here, we employed a pooled multiplexing strategy to profile 399,683 g
pmc.ncbi.nlm.nih.gov
먼저 연구 배경부터 설명하도록 하겠습니다. 선행 연구들에서는 GWAS를 통해 위암(GC, gastric cancer)의 유전적 기반(genetic basis)을 조사하여 20개가 넘는 suceptibility loci를 파악하였고, 해당 loci는 non-coding genomic region(특히, regulatory elements)에 위치한 유전 변이(genetic variants)인 SNP들로 대표됨을 밝혀냈습니다. 그렇지만, 어떤 유전 변이가 어떤 생물학적 경로를 통해 위암 발달에 영향을 미치는지, 정확한 메커니즘은 아직 밝혀지지 않은 상황입니다.
무엇보다 위암 발달에 대한 유전변이의 영향을 밝히는 대부분의 연구들에서 bulk transcriptome 데이터(샘플 내 전체 세포들의 유전자 발현 평균 데이터)를 사용했기에 조직 내 개별 세포들 간의 유전자 발현의 이질성을 반영하지 못했습니다. 특히 위의 상피는 표면 점액세포(mucous pit cells), 줄기 세포(stem cells), 벽 세포(parietal cells), 경부 점액세포(mucous neck cells), 주세포(chief cells) 등 다양한 세포 유형들로 구성되어 있기에, 유전자 발현에 있어 위 세포 유형별로 특이적인 유전적 조절 메커니즘을 조사해야 생물학적 메커니즘을 제대로 규명할 수 있을 것입니다.
최근에는 조직 내 개별 세포들 사이에서의 유전자 발현의 이질성 즉, cell-type specific gene expression pattern을 파악할 수 있는 single cell RNA-seq 기술이 등장하였습니다. 따라서 이제 우리는 본 기술을 통해 조직 내 세포 조성을 파악하는 cell type annotation을 한 뒤 질병 관련 유전 변이(disease-related genetic variants) 데이터와 통합함으로써 어떤 세포 유형에서 유전 변이가 특이적으로 영향을 가하는지 이해할 수 있습니다.
유전자 발현에 영향을 미치는 유전 변이가 위치하는 게놈 영역을 eQTL(expression quantitative trait loci)이라고 합니다. 그리고 이러한 eQTL에서 발현이 조절되는 변이는 ‘eSNP’라고 하죠. 참고로, GWAS를 통해 특정 형질(ex. 질병)과 통계적으로 유의하게 연관된 SNP 집합을 찾을 수 있는데요. 만약 동일한 SNP에서 eQTL 신호와 GWAS 신호가 모두 검출된다면(co-localization), 해당 SNP가 GWAS 형질(phenotype)에 영향을 가하는 메커니즘이 유전자 발현을 매개한다고 볼 수 있습니다.
세포 유형 특이적으로 유전자 발현에 영향을 가하는 유전변이가 존재하는 영역을 'cell-type-specific eQTL'이라고 합니다. 본 연구에서는 위(stomach) 조직 샘플들을 기반으로 'cell type-specific eQTLs'를 파악한 것은 물론, 위암(gastric cancer, GC)에 대한 GWAS 신호와 eQTL 신호가 함께 나타나는 특정 세포 유형에서의 공위치화 신호(co-localization signal)를 밝혀냈습니다.
또한, cell-type-specific transcriptome-wide association study(TWAS)를 수행함으로써 특정 세포유형에서 위암 위험(GC risk)과 연관된 유전적 조절을 받는 유전자들을 찾아냈는데요. 이를 통해 본 연구팀은 단일 세포 수준에서의 위암과 연관된 유전자 발현 조절 메커니즘 규명에 기여할 수 있었습니다. 이외에도 연구팀은 서로 다른 위 세포 유형들(gastric cell types)에 대한 eQTLs와 eGenes(eQTLs로부터 발현 조절을 받는 유전자)를 검색하고 확인할 수 있는 'scGaTE'라는 공공 엑세스 온라인 플랫폼을 개발하였습니다.
Single-cell Gastric Transcriptome Explorer
ccra.njmu.edu.cn
결국 본 연구의 목적은 scRNA-seq 데이터를 사용해서 위, 보다 나아가 위암 맥락에서 유전 변이가 유전자 발현에 미치는 영향을 단일세포수준에서 탐구하기 위함이었으며, 이를 통해 위암 risk와 연관된 유전자 발현에 대한 유전적 조절에 있어 개인 간 변동성(inter-individual variation)을 이해하고자 했습니다. 그렇다면, 본 연구팀이 어떻게 연구를 설계했는지 간단하게 설명해보도록 하겠습니다.
먼저 본 연구팀은 chinese ancestry인 203명의 사람들에게서 위점막 조직을 채취한 뒤 8-9명씩 27개 풀로 섞어 scRNA-seq 분석을 수행하였습니다. 이때, 동일 기증자를 대상으로 SNP genotyping도 함께 진행하여 이를 토대로 디멀티플렉싱을 할 수 있었고, 결과적으로 총 399,683개의 고품질 세포(19개 아형)를 기증자별로 재구성하였습니다. 이후 canonical markers를 기반으로 cell type annotation을 거쳐 19개의 특징적인 위 세포 아형들을 파악하였습니다*. 다음으로 본 연구팀은 특정 세포 아형의 조성 비율에 유의하게 기여하는 macroscopic factors를 확인한 것은 물론, cell type을 phenotype으로 두고 이와 연관된 SNP 집합을 GWAS를 통해 파악했습니다. 이후 위 세포 아형 조성 비율과 연관되어 나타난 7가지 macroscopic factors와GWAS를 통해 파악된 세포 아형 조성 비율과 연관된 SNP 집합이 gastric cell composition 분산에 기여하는 정도를 파악하였습니다.
* 참고로 전체 세포들 중에서 mucous neck cells와 mucous pit cells로 annotation된 세포가 각각 약 25%, 약 17.4%로 높은 비율을 차지했다고 합니다.
Figure - PMC
Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites.
pmc.ncbi.nlm.nih.gov
다음으로 본 연구팀은 공변량(성별, 연령, H. pylori 감염, 10 개의 유전적 PC, 2개의 PEER factor)에 의한 교란을 제어하면서 세포 유형 특이적으로 유전자의 발현을 유의하게 조절하는 'cell-type-specific eQTL'을 파악하는 eQTL 분석을 진행했습니다. 이후 위암(GC)에 대한 GWAS 신호와 eQTL 신호가 함께 나타나는 특정 세포 유형에서의 공위치화 신호(co-localization signal)를 밝힌 것은 물론, TWAS를 통해 해당 신호로부터 발현이 조절되는 특정 세포 유형에서의 GC risk와 연관된 유전자들을 파악하였습니다. 마지막으로는 해당 유전자들의 기능을 세포 유형 맥락에서 해석하였고, 이를 통해 위암 위험에 영향을 가할 수 있는 세포 유형 특이적인 분자 기전을 제시하였습니다. 결국, 본 연구는 단일세포 수준에서 위암 위험과 관련된 유전 변이에 의해 발현이 조절되는 질병 감수성 유전자들을 파악하면서, 특정 세포 유형에서 해당 유전자들의 기능을 위암 맥락에서 해석했다고 볼 수 있겠습니다.
이제 구체적인 연구 설명을 해보도록 하겠습니다. 먼저 본 연구팀은 GC risk(위암 위험)와 연관있는 18개의 거시적인 요인들(macroscopic factors)*과 위암 세포 조성 간의 관계를 분석했습니다. 결과적으로 18개 요인 중 Helicobacter pylori infection, gastric lesion, sex, cereal fiber intake, smoking, irregular meals, 그리고 spicy diet가 위 세포 아형 조성 비율과 상당히 연관되어 있음이 나타났는데요. 구체적으로, H.pylori-infected individuals에서 CD4+ T cells의 확장이 나타났고, 상피 아형 집단들(epithelial subpopulations)에서는 표재성 위염(superficial gastritis, SG)에 비해 만성 위축성 위염(chronic atrophic gastritis, AG) 조건에서 해당 아형 집단들의 구성 비율이 감소함을 확인했습니다. 또한, 기존에 보고되지 않았던 결과로 여성에게서 남성에 비해 내분비 세포(endocrine cells)의 확장이 특이적으로 나타났다고 하죠.
* 18 macroscopic factors : age, sex, Body mass index (BMI), H. pylori infection, gastric lesion, smoking, alcohol drinking, tea consumption, consumption of fresh fruits and vegetables, cereal fiber intake, meat intake, fried foods intake, pickled foods intake, spicy foods intake, regulation, temperature, speed of meal
이러한 macroscopic factors와 위암 세포 조성 비율 간의 관계는 single-cell cluster association model인 MASC(v0.1.0) 모델을 통해 파악할 수 있었습니다.
GitHub - immunogenomics/masc: MASC: Mixed-effects Association testing for Single Cells
MASC: Mixed-effects Association testing for Single Cells - immunogenomics/masc
github.com
그렇다면 여기서 MASC 모델이 무엇이며 어떻게 분석에 적용되었는지를 잠깐 설명하고 넘어가도록 하겠습니다. MASC(Mixed-effects Association testing for Single Cells)는 “이 세포가 해당 클러스터(= cell type)에 속하는가(1) / 아닌가(0)” 즉, 개별 세포가 특정 세포 아형 군집에 속할 확률을 종속변수로 두는 로지스틱 혼합모형을 19개의 세포 아형마다 적합합니다. 식은 다음과 같습니다.
여기서 식 절편 β₀은 모든 공변량이 0인 기준 상태에서 세포가 특정 세포 아형 군집에 속할 기본 log-odds를 나타내고, β₁은 관심 공변량(예: H. pylori 감염 여부)이 한 단위 증가하거나 음성에서 양성으로 바뀔 때 log-odds가 얼마나 변하는지를 나타냅니다. 다음으로 Xᵢγ 항은 나이, 성별, 배치 등 다른 교란 요인을 고정효과로 보정해 β₁ 추정을 순화하고, bᵢ는 동일 기증자(혹은 샘플)에서 얻은 세포들이 공유하는 내부 상관을 난수효과로 모델링해 과도한 제1종 오류를 방지합니다*. 이렇게 적합된 모형에서 β₁ = 0 가설을 우도비 검정으로 평가해 'cell-type x 공변량' 별 p-값을 구하고, 이를 통해 단일세포 수준 공변량 연관성을 정량화합니다.
* 참고로 제1종 오류는 실제로 효과가 없는데도 통계 검정 결과 귀무가설을 기각하는 거짓 양성(즉, 존재하지 않는 차이나 연관이 있다고 잘못 판단하는 경우)을 말합니다.
예를 들어 설명해보겠습니다. 위 조직 샘플들 가운데 특정 샘플의 기증자는 H.pylori 양성, 다른 샘플의 기증자는 음성이라고 해봅시다. 이후 이들의 scRNA-seq 데이터를 pit cell, chief cell, parietal cell 등으로 annotate한 뒤, 먼저 pit cell에 대해 '이 세포가 pit cell 클러스터에 속할 확률'을 종속변수로 로지스틱 혼합모형을 적합합니다. 여기서 로지스틱 혼합 모형은 0/1 결과의 로그 오즈를 고정효과(공변량)와 무작위효과(군집 간 변이)로 동시에 설명해주는 이분형 회귀 모델을 말합니다. 이때 한 번은 감염 변수(공변량)를 포함한 full 모델로 H. pylori 감염 여부(0 또는 1)를 교정효과로 고려하고, 다른 한번은 감염 변수를 제거하는 null 모델로 기증자 절편을 무작위효과로 고려해 계산합니다. 이 두 모델의 로그 우도를 비교하면 우도비 통계량 2Δℓ이 나오는데 만약 12.3으로 나왔다고 하면 이는 자유도 1인 χ² 분포 기준 P = 0.00045에 해당합니다. 이 값은 0.05보다 훨씬 작으므로 귀무가설 𝛽₁ = 0 즉, '감염 여부가 Pit 세포 비율에 영향을 주지 않는다'는 가설을 기각하고, 감염된 샘플에서 pit 세포가 통계적으로 유의하게 더 많다는 결론을 내립니다. 같은 절차를 Chief, Parietal 세포 등에서도 적용해 얻은 p값들을 감마 결합하면 전역 p값이 나오고, 그 p값 또한 0.05보다 작으면(ex. 0.0013) 이는 H.pylori 감염이 위 상피 세포 구성을 전반적으로 재편한다고 볼 수 있습니다.
결과적으로 본 연구팀은 MASC를 통해 "H. pylori 감염은 Pit 세포를 늘리고, Parietal 세포에는 거의 영향을 주지 않지만, 19개의 세포형 모두를 합친 전역 검정은 P < 0.001(***)로 유의하다"처럼 세포 특이적, 행 전역(19개의 세포형 전부) 두 수준의 정보를 동시에 제공하였습니다.
추가적으로 본 연구팀은 모든 위 세포 유형 annotation 결과를 phenotype으로 두고 GWAS를 수행하였는데요. 그 결과 서로 다른 세포 아형들의 조성 비율과 연관된 68개의 independent genetic loci를 발견했습니다. 예를 들어, 18번 염색체에 위치한 세 개의 유전 변이는 M1-like macrophage의 비율과 상당히 연관되어 있었습니다(rs8087289의 A allele는 비율 증가, rs150615809의 A allele는 비율 감소와 연관됨). 여기서 연구팀은 먼저 각 기증자에서 특정 세포 유형이 차지하는 비율을 구한 뒤, 그 값이 0-1 구간 사이에 치우쳐 있는 문제를 완화하기 위해 arcsine-squeare-root 변환을 적용했습니다(빈도 자료에 대해 분산을 일정하게 만들고 정규분포(특히 중간구간에서)에 가까워지도록 설계한 변환). 변환 후에도 여전히 0값이 많아 정규성을 만족하지 못한 희귀 세포 군집(Enterocytes, CD4, Naive T, M1·M2 대식세포, Arterial ECs)은 분포를 5분위로 나눠 범주형(0~4) 변수로 치환하였습니다. 그리고나서 PLINK v1.9를 이용해 GWAS 분석을 진행했는데요. 구체적으로, 변환된 세포 비율을 종속변수로, SNP 대립유전자 수를 주효과(β₁)로 두고 연령, 성별, H.Pylori 감염 여부, 그리고 집단 구조를 보정하기 위한 유전적 주성분 10개(PC1 ~ PC10)을 공변량으로 포함한 선형 회귀를 세포 유형마다 수행했습니다. 다시 말해, 변환된 세포 비율에 대한 SNP 연관성을 연령, 성별, H.pylory 상태 등의 요인을 공변량으로 두고 테스트한 것입니다.
이때, 모든 SNP에 대해 얻은 p값 가운데 상호 연관성이 낮도록 LD-clumping(r² < 0.5, 250 kb 창)으로 추린 8498개의 독립적인 대표 SNP를 대상으로 연구-전 보정 α = 0.05의 Bonferroni 기준을 적용해 유의수준 6.0x10^-6으로 설정했고, 이를 넘는 SNP를 해당 세포 유형의 조성 비율에 유의하게 영향을 주는 변이로 최종 보고하였습니다. 이러한 방법을 통해 연구팀은 M1 macrophage 비율이 18번 염색체의 특정 SNP인 rs9963169와 강하게 연관되어 있음(β = 0.94, P = 1.67 x 10^-8)을 밝혀낸 것이죠. 종합적인 결과는 아래 링크 Figure 2D에서 확인할 수 있습니다.
Figure - PMC
Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites.
pmc.ncbi.nlm.nih.gov
Figure2D는 특정 염색체 위치에 존재하는 유전 변이와 세포 비율 사이의 연관성을 보여주는 Manhattan plot인데요. 같은 염색체 안에서 -log10P 값(y축)이 특히 높게 솟은 구간(피크)은 그 위치에 해당 세포형 비율에 영향을 주는 유전 변이가 존재한다는 뜻입니다. 예를 들어, 노란색(M1 macrophage) 점이 18번 염색체 근처에 군집해 있다는 것은 그 염색체 영역에 존재하는 변이들이 위장 조직 내 M1 대식세포 비율을 변화시킬 가능성이 높다는 것을 시사합니다. 또한, 15번 염색체 영역처럼 서로 다른 색이 같은 피크에 섞여 있다면, 하나의 변이가 여러 세포형의 비율에 동시 영향을 미치는 다면 발현(pleiotropy)하고 있음을 시사합니다. 종합적으로, 본 Manhattan plot은 어떤 게놈 영역에 위치한 유전 변이가 어떤 위장 세포 아형의 상대적 풍부도에 영향을 주는지, 그리고 본 영향이 세포 유형 특이적인지를 직관적으로 보여준다고 할 수 있겠습니다.
다음으로 본 연구팀은 위 세포 아형 조성과 연관이 있다고 나타난 7가지 macroscopic factors*와 genetic factors(GWAS를 통해 파악한 gSNPs)의 위 세포 조성(gastric cell composition)의 분산에 대한 기여를 평가했습니다. 여기서 각 예측변수가 세포 조성 분산(=변동성)에 얼마나 기여하는지를 평가할 땐, relimpo R 패키지의 lmg 지표(Lindeman–Merenda–Gold 방법)를 사용하였다고 하는데요. 이는 모든 가능한 변수 순서를 고려해 각 변수의 제곱합(SS)을 계산하고 이를 평균함으로써 순서 의존성을 제거한 공평한 분산 분해 값을 산출하는 절차입니다.
* H. pylori infection, gastric lesion, sex, cereal fiber intake, smoking, irregular meals and spicy diet
여기서 순서 의존성을 제거하는 이유는 분산을 계산할 때 먼저 변수 순서를 고려해 단계적 제곱합 SS(추가되는 변수의 R^2 증가분 x 총제곱합 SST)을 구하기 때문입니다. 즉, 각 순열(만약 설명변수가 3개면 3x2x1=6가지 순열)에서 '앞에 이미 들어간 변수들을 모두 포함한 모델'과 '거기에 방금 추가된 변수까지 포함한 모델'의 R^2 차이를 해당 변수의 기여도로 기록하기에, 각 변수별로 모든 순열에서 자신이 마지막으로 들어가면서 확보한 ΔR²을 평균해 변수 입력 순서가 설명력 분배에 미치는 편향을 없애는 것입니다. 이렇게 얻은 lmg 값은 각 요인이 해당 세포 아형 비율의 총 분산 중 얼마만큼을 설명하는지를 나타냅니다. 연구팀은 이를 통해 7가지 거시적인 요인과 유전적 요인 각각이 세포 조성 변동성에 미치는 상대적 영향력을 정량화했습니다. 결과적으로 해당 요인들의 기여는 세포 유형들마다 다르게 나타났으며(15.4% ~ 67.0%), genetic factors(9.5%–37.6%)와 H. pylori infection(0.1%–21.8%)이 세포 유형 조성 분산 비율에 대한 기여가 높았다고 합니다.
다음으로 본 연구팀은 앞서 말했듯 eQTL 분석을 수행하여 5,443개 유전자의 발현을 조절하는 8,498개의 독립적인 eQTL을 밝혔습니다. 여기서 eQTL로부터 조절 받는지 확인할 후보 유전자를 선택할 때 다음과 같은 기준을 고려했습니다. 먼저 각 세포 유형마다 Seurat으로 정규화-계수(UMI) 행렬을 만들고, 1% 이상의 세포에서 UMI > 0으로 검출된 유전자만 cis-eQTL 분석에 투입했습니다. 또한, 특정 세포 유형 중에서 해당 기증자에게 최소 5개의 세포가 확보된 세포 유형이 유전자-세포형 테스트에 포함됐습니다. 그렇게 유전자 전사 시작점(TSS) ± 1 Mb 범위의 모든 SNP에 대해 FastQTL 선형 회귀를 돌리고, 공변량(성별·연령·H. pylori 감염, 10 개의 유전적 PC, PEER factor 2개)을 조정했습니다. 본 식은 다음과 같습니다.
Gx : 기증자 별 gene X의 평균 발현으로 구성된 matrix
β0 : 모든 설명변수(유전자형, 공변량)가 0일 때, 기증자의 예상 geneX의 평균 발현 값.
geno : 각 기증자별 해당 구간 SNP들의 유전자형(0/1/2 대립유전자 수)
β1 : SNP의 대립유전자 수(0, 1, 2)가 한 개 증가할 때 geneX 발현이 평균적으로 변화하는 정도로, β1 > 0이면 발현이 증가하고 β1 < 0이면 발현이 감소함. 이 값의 유의성(β1 = 0 vs ≠ 0)이 바로 eQTL 판별 기준.
β2 : 성별, 연령, H.pylori 감염, 10개 유전적 PC, 2개 PEER factor 등 각 공변량이 발현에 미치는 효과를 나타내는 계수들의 집합으로, 예로, β2.age는 나이가 1살 늘 때 평균 발현 변화량을 뜻함.
C : sex, age, H. pylori infection, PCs, 그리고 PEER factors를 포함하는 공변량 매트릭스
ε : 모델이 설명하지 못한 오차, 개체 특이적 변동. 평균 0, 분산 𝜎2를 가정
여기서 β1이 통계적으로 유의하면 해당 SNP가 특정 세포형에서 특정 유전자 발현을 독립적으로 조절하는 cis-eQTL(eSNP)로 간주되고, β2 항은 인구 구조, 배치 효과 등 교란 요인을 보정해주어 β1 추정치를 편향 없이 해석할 수 있게 해줍니다. 여기서 공변량으로 선택한 '10개의 유전적 PC'와 'PEER factor'에 대해 잠깐 설명하고 넘어가도록 하겠습니다. 먼저 기증자별 SNP 행렬에 대해 주성분분석(PCA)을 수행하면, 대립유전자 빈도의 공분산 구조를 요약하는 선형 축(PC)이 얻어집니다. 이때 PC1, PC2는 흔히 '유럽 <-> 동아시아'처럼 계통 및 인구 집단(ancestry) 차이를 설명하는 축이고, 뒤쪽 PC들은 보다 미세한 혈연 또는 혼합(admixture) 정보를 담는다고 하는데요. eQTL 회귀식에 이 축들을 공변량으로 넣으면, 유전자형과 유전자 발현이 같은 조상의 사람끼리 비슷하다는 구조적 상관(집단 분할, stratification) 때문에 생길 수 있는 가짜 상관을 효과적으로 제거할 수 있습니다. 예로, 어떤 SNP와 MUC20 유전자 발현이 양(+)의 상관을 보였지만, 실제로는 모두 '동아시아 계열 기증자'에게서 높은 값이었다고 가정해봅시다. 이때 PC1, PC2를 공변량에 넣으면 '동아시아 vs 기타' 축을 회귀식이 설명하므로, SNP-유전자 상관은 사라져 가짜 eQTL을 방지할 수 있습니다. 참고로 본 연구에서 공변량으로 10개의 PC를 선택한 것은 GTEx 등 대규모 eQTL 연구에서 PC 5~10 정도까지 넣으면 추가 PC의 편익까지 크게 줄어든다는 경험적 기준을 따른 것이라고 하네요.
다음으로 2개 요인을 사용한 PEER(Probabilistic Estimation of Expression Residuals)는 발현 행렬을 베이지안 요인분석으로 분해해 실험배치, 시퀀싱 깊이, 세포 품질, 세포주기 상태처럼 측정되지 않은(숨은) 변동원을 잠재 요인(latent factor)으로 추정합니다. 연구팀은 먼저 mucous-neck cell에서 1~10개의 PEER factor를 하나씩 추가해보며 eGene 수 변화를 확인했고, 두번째 factor까지 넣었을 때 잡음이 최대한 제거되면서도 검정력이 크게 감소하지 않았다고 보고하였습니다. 따라서 최종 모델에는 PEER1, PEER2 두 개만 포함해 과보정(over-fitting)을 피했다고 하죠. 이렇게 적절한 요인 개수를 설정해 과보정을 피하면서 배치 차이, 노이즈 등을 잡아내면 실제로 보고자 하는 변화와 무관한 신호를 제거해 실제 생물학적 eQTL을 더 잘 검출할 수 있습니다. 요약하면 유전적 PC 10개는 '인구 구조로 인한 교란변수'를, PEER factor 2개는 '숨겨진 기술적/생물학적 변동'을 보정해 eQTL 분석의 거짓 양성률(FDR)을 낮추고 검출력을 높이기 위한 필수 보조 변수로 역할을 한다고 볼 수 있겠습니다.
이후 연구팀은 SNP-유전자 쌍의 nominal p-값을 전 유전체 False Discovery Rate(FDR) 절차로 교정한 뒤 FDR ≤ 0.05를 만족하면 해당 SNP를 cis-eQTL(eSNP)로, 그로부터 조절받는 유전자는 eGene으로 정의했습니다. 구체적으로 연구팀은 각 유전자-SNP 조합마다 선형 회귀를 돌려 얻은 nominal p-값(β₁ = 0 가설의 유의성에 대한 p-값)을 하나의 긴 목록으로 모으고, FastQTL의 기본 절차를 따라 두 단계로 다중검정을 제어했습니다. 첫번째 단계는 유전자 내부(per-gene) 보정으로, 유전자별로 최대 10,000회(FastQTL 기본값) 유전자형 행을 무작위 재배열해(적어도 인구 구조, 공변량은 그대로 유지) 회귀를 반복 수행했습니다. 이렇게 얻은 permutation 분포로부터, 각 SNP-유전자 쌍의 nominal p가 '무작위일 때 이 값보다 작을 확률'에 해당하는 empirical p-값으로 변환되었습니다. 이 단계는 유전자마다 시험 횟수가 다르므로, 같은 nominal p라도 쉽게 나왔는지, 어렵게 나왔는지를 균형 있게 비교하려는 목적입니다.
두번째는 전 유전체(global) 보정 단계로, 모든 유전자-SNP 쌍의 emprical p를 일렬로 정렬해 Benjamini–Hochberg(BH) 알고리즘을 적용했습니다. BH 알고리즘은 정렬된 순위를 이용해 각 시험의 q-값(FDR에 대한 최소 α)을 계산하고, q ≤ 0.05(즉, 5 % FDR 이하)인 쌍을 significant cis-eQTL로 선언했습니다. 이때 유전자-SNP쌍이 기준을 통과하면 해당 SNP는 eSNP, 그 유전자는 최소한 하나의 eSNP를 가지므로 eGene으로 규정됩니다. 예를 들어 MUC20 주변 1 Mb 창에 3,200 개의 SNP가 있다면, permutation 단계에서 이들 3,200 개의 nominal p가 모두 empirical p로 바뀝니다. 이후 모든 유전자 합계 약 20 만 개(모든 세포형 합쳐 수억 개) 시험 중, BH 절차를 거친 뒤 MUC20-rs98765의 q-값이 0.012라면 q ≤ 0.05이므로 rs98765는 eSNP, MUC20은 eGene이 됩니다. 만일 MUC20 주변 다른 SNP들이 q > 0.05이면, eSNP는 rs98765 하나뿐이지만 eGene 지위는 유지됩니다.
이후에는 조건부 분석(conditional analysis)이 이어졌는데요. 가장 작은 p값을 가진 SNP(rsX)를 eSNP1으로 지정한 뒤, 다시 같은 모델을 돌리되 rsX의 유전자형을 공변량에 추가해 그 효과를 제거하고 남은 SNP들을 재검정할 수 있습니다. 이렇게 하여 독립된 두번째 신호(rsY•eSNP2)가 발견되면 또 다시 rsX/rsY를 함께 보정해 세번째 신호를 찾는 방식으로, FDR ≤ 0.05를 만족하는 독립 신호가 더 이상 없을 때까지 반복했습니다.
결과적으로 연구팀은 19개의 위장 세포아형에서 5,443개의 eGene과 8,498개의 독립 eSNP를 확보했으며, 이는 5% FDR 기준(거짓 양성률 5 % 이하의 신뢰 수준)을 만족하는 cis-조절 신호 집합이라고 할 수 있겠습니다. 여기서 어떤 유전자(eGene)는 하나의 강력한 eSNP만 가졌고(단일 신호), 다른 유전자는 순차 조건부 분석을 통해 두 개 이상의 독립적인 변이가 존재함이 밝혀져(다중 신호) eGene 개수를 넘어서는 cis-eQTL 개수가 나타나게 되었습니다. 이와 더불어 본 연구팀은 앞서 말했듯 'scGATE'라는 플랫폼을 개발하여 서로 다른 위 세포 아형들에서 유의하게 나타난 eQTLs와 eGenes를 간편하게 브라우징할 수 있도록 했습니다.
Single-cell Gastric Transcriptome Explorer
ccra.njmu.edu.cn
다음으로 본 연구팀은 위 세포유형들 사이에서 공유되는 cis-eQTL signal의 수를 비교했습니다. 결과적으로 eQTL의 수는 세포 유형마다 다르게 나타났는데요. mucous neck cells에서의 1,915개부터 myfibroblasts에서의 24개까지 그 스펙트럼이 다양했습니다. 이는 cell-type specifi eQTL의 수가 많다는 결과로도 볼 수 있는데, 실제로 약 80% 이상의 eQTLs가 위 세포 아형 특이적으로 나타났다고 합니다. 본 결과는 아래 Fig3B에서 확인하실 수 있습니다.
Figure - PMC
Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites.
pmc.ncbi.nlm.nih.gov
참고로 Fig3B에서 막대 하나(upset plot의 각 열)는 어떤 세포 유형 조합에서 동일한 eQTL 신호가 관찰됐는지를 나타내는 집합(intersection) 단위입니다. 즉, 막대는 어떤 세포형 조합이 유전자 발현 조절 신호를 공유하는지를 패턴별로 구분한 것이며, 왼쪽의 색 막대들은 각 세포형이 자기만의 eQTL을 갖고 있는 경우이고, 오른쪽의 검은 막대들은 여러 세포형이 동일 변이를 함께 사용하는 경우를 집계한 결과입니다.
각 세포 유형에서 파악된 eQTLs의 개수는 세포 유형의 풍부도(cell type abundance)는 물론, 샘플( individuals)의 수 모두와 정적인 상관관계가 나타났다고 하는데요. 이는 세포 수나 샘플의 수가 늘어날 수록 발견되는 eQTL의 수가 늘어난다는 기존에 알려진 사실과 부합하는 결과로 볼 수 있습니다. 무엇보다 eQTLs의 약 80%가 cell-type-specific함에도 불구하고 오직 92개의 eGenes(1.7%)만 개별 세포 유형에서 독점적으로 발현되고 대다수 eGene이 여러 세포 유형에서 유의하게 발현됨이 확인되었습니다(FigS4A 참조). 이는 cell-type-specific eQTL이 단순한 유전자 발현량 차이만으로는 설명되지 않음을 시사합니다. 무엇보다 UMI count가 0보다 크다는 정보만으로는 발현 유무를 넘어선 조절 민감도를 포착할 수 없죠. 따라서 본 연구에서는 eSNP를 대상으로 LD 구조와 조건부 회귀 분석을 수행함으로써 효과 크기가 세포 유형별로 상이함을 확인하였으며, 이는 크로마틴 접근성, 전사인자 결합 조합, 후성유전학적 상태 등 세포 고유의 전사 환경이 eQTL 특이성을 규정하는 주요 요인임을 뒷받침합니다.
다음으로 본 연구팀은 19개의 위 세포 아형들마다 공유되는 혹은 고유한 조절 메커니즘이 있는지 조사하기 위해 세포 유형 쌍(19 x 18 = 342 cell type pairs)마다 eGene의 발현이 동일한 변이에 의해 조절되는지 혹은 서로 독립적인 변이에 의해 조절되는지를 조건부 회귀분석을 수행하여 파악했습니다. 구체적으로, 세포 유형 A에서 lead eSNP를 이용해 eGene 발현량(y)과 유전형(geno) 간 단일 선형 회귀를 수행하여 원래 회귀 계수 borig를 추정한 뒤 이 모델에 세포 유형 B의 lead eSNP를 공변량으로 추가하고, 다시 회귀를 수행해 새로운 회귀 계수 bcond와 p 값을 계산하였습니다. 여기서 효과 크기 변화는 FC = bcond/borig로 정의하며, 유의성: p ≤ 0.05, 크기 차이: −log₂ FC ≥ 2를 동시에 만족하면 '공유된(shared) eQTL'로, 그렇지 않으면 세포 유형별 '독립(independent) eQTL'로 분류했습니다. 본 연구팀은 이러한 절차를 가능한 342개 세포 유형 쌍에 모두 적용하여, 각 eGene이 서로 다른 세포 유형의 서로 다른 변이에 의해 태그되는지를 일괄적으로 판별했습니다.
결과적으로, 대부분의 epithelium eQTLs(enterocytes 제외)는 adjustment 후 회귀 계수(regression coefficients)에서의 상당한 변화를 보였으며, 이는 epithelial cells 사이에서의 공유된 eQTL 패턴을 암시한다고 볼 수 있습니다. eGene 발현에 대한 A 세포의 리드 eSNP 신호가 B 세포의 리드 eSNP를 통계모델에 추가(조건부)했을 때 유의하게 변화한다면(약해진다면), 두 SNP가 서로 독립적인 조절 요소가 아니라 같은 표적 유전자 발현 변화를 매개하는 동일 신호(tagged signal)을 포착하고 있다는 뜻이기 때문입니다. 그렇게 epithelium 사이에서 공유된 eQTL 신호를 보인 eGenes(GSTM1, GSTM3, and GPX)들은 글루타치온 대사 및 유기산 생합성 경로에 enrich하게 나타났다고 합니다.
이때 세포 유형 간 공유 eGenes의 eSNP들은 LD r² > 0.5로 강한 연관성을 나타내는 반면, 세포별로 독립적인 eGenes의 eSNP들은 LD r² ≤ 0.5로 낮은 연관성을 보였습니다(Figure S4F 참조). 공유 Genes를 태그하는 서로 다른 리드 eSNP 쌍이 높은 LD를 보인다는 것은 두 세포 유형에서 각각 검출된 eQTL 신호가 서로 독립적인 변이가 아니라 동일한 인과 변이를 반영하고 있음을 의미합니다.
구체적으로 보면, LD(Linkage Disequilibriun)란 한 위치의 대립형질(allele)이 다른 위치의 대립형질과 얼마나 함께 유전되는지를 나타내는 지표이므로, r^2이 높다는 것은 두 SNP가 같은 하플로타입 블록(haplotype block) 내에 있어 한 변이가 곧 다른 변이를 강하게 예측한다는 뜻입니다. 즉, A 세포에서 lead eSNP1로 검출된 신호와 B세포에서 lead eSNP2로 검출된 신호가 실질적으로는 동일한 인과 조절 변이(causal regulatory variant)를 태그하기 때문에 두 세포에서 같은 유전적 효과가 공유된 것으로 해석할 수 있습니다 .반면, LD가 낮은 독립 신호들은 서로 다른 하플로타입 상의 별개 변이가 각 세포에서 eQTL 효과를 내고 있음을 보여주어, 세포 특이적 조절 메커니즘이 다름을 시사합니다. 이는 크로마틴 접근성, 전사인자 결합, 후성유전학적 표지 등 세포마다 다른 전사 환경이 유전변이의 조절 효과를 증폭 혹은 완화시키기 때문으로 추측할 수 있겠습니다.
종합적으로 연구진은 우선 각 eGene에 대해 리드 eSNP를 조건부 회귀 및 LD 분석으로 교차 검정하여 조건화 이후에도 효과 크기가 유지되면 독립 신호, 감소하면 동일 신호로 간주하는 방식으로 342개의 세포쌍에서 실제로 공유되는 변이를 리스트화했습니다. 이후 본 연구팀은 어떤 조절 변이(regulatory variants)가 각 세포 유형에서의 eQTL profiles 중에서 공유되는지 multivariate adaptive shrinkage(mashR)를 통해 평가했습니다. 즉, FastQTL 결과를 mashR로 변환해 모든 유의 SNP-gene 쌍의 효과량 벡터를 동시에 모델링한 뒤, 방향이 같고 크기가 0.5배 이내로 유사한 경우를 공유 효과로 정의하면서 세포 간 eQTL 프로필의 전반적 중첩도를 추정했습니다. 그 결과 Figure S4G에서 보이듯 epithelium과 stroma cells의 세부집단이 많은 수의 eQTLs를 공유한 반면, 면역계 세포는 공유도가 낮게 나타났습니다. 이는 상피와 기질 계열이 위 조직 맥락에서 유사한 전사 및 후성유전 환경을 지닌다고 추론할 수 있겠습니다. 결국 세포 고유의 크로마틴 접근성, 전사인자 결합, 후성유전 상태가 eQTL 특이성을 좌우한다는 점이 변이 리스트와 공유 지표 두 수준에서 모두 확인된 셈이라고 할 수 있겠습니다.
글이 길어지니 다음 글에서 이어 정리해보도록 하겠습니다. 본 글에는 오류가 많을 수 있다는 점 양해부탁드리겠습니다. 감사합니다!
참고자료
1. Bian, Lijun et al. Single-cell eQTL mapping reveals cell-type specific genes associated with the risk of gastric cancer, Cell Genomics, Volume 5, Issue 4, 100812
2. ChatGPT
