07. 한 샘플에서 아웃라이어 유전자 찾기
RNA-seq로 발현량 행렬을 얻은 뒤, 환자 한 명의 발현량 벡터 하나만 있을 때, 공개 레퍼런스를 이용해 비정상적으로 튀는 유전자를 어떻게 고르는지에 집중합니다.
이 노트의 목표: N=1 분석에서 비교 기준이 왜 필요한지, 백분위·z-score·fold change·max 비교가 각각 무엇을 말하는지, 그리고 “발현이 튄다”를 어디까지 해석할 수 있는지 이해한다.
운영 중인 서버 한 대의 지표만 보고 이상 여부를 판단할 수 없는 상황과 같습니다. CPU 80%가 높은지 알려면 같은 역할의 서버가 보통 어느 범위에서 움직이는지 알아야 합니다. 여기서는 서버 지표 벡터 대신 샘플의 유전자 발현 벡터, 과거 관측 분포 대신 정상조직·암 레퍼런스 코호트를 사용합니다.
for gene in common_genes: baseline = reference_samples[gene] score[gene] = compare(patient_sample[gene], baseline)
candidates = rank(score)이 의사코드의 출력은 원인 판정이 아니라 추가로 조사할 후보 순위입니다.
1. 한 샘플로 무엇을 알 수 있나: 샘플이 딱 하나일 때
섹션 제목: “1. 한 샘플로 무엇을 알 수 있나: 샘플이 딱 하나일 때”지금까지 만든 것은 환자 한 명의 종양 조직 RNA-seq 프로파일 하나 입니다. 유전자 수만 개(주석 기준 약 6만 개) 각각의 발현량(TPM 또는 raw count)이 담긴 벡터 하나.
DataFrame으로 보면 행은 유전자, 열은 샘플, 값은 발현량입니다. 환자 한 명은 열 하나이고, 같은 유전자의 레퍼런스 값들은 가로로 긴 분포를 만듭니다.
| gene | patient | normal_001 | normal_002 | cancer_001 | … |
|---|---|---|---|---|---|
MDM2 | 1235.72 | 4.10 | 5.32 | 12.80 | … |
ACTB | 850.31 | 790.20 | 910.44 | 822.17 | … |
그래서… 이 중 뭐가 중요한가? 발현량 6만 줄을 눈으로 다 볼 수는 없습니다. 필요한 질문은 하나로 좁혀집니다:
“이 종양에서, 정상 조직이나 다른 암과 비교했을 때, 비정상적으로 발현이 튀는 유전자는 무엇인가?”
이 질문에 답하는 작업이 아웃라이어 유전자(outlier gene) 탐색입니다.
아이디어. 가로축은 한 유전자의 발현량(TPM), 세로축은 그 발현량을 가진 샘플 수. 대부분의 종양(All tumors)은 MDM2가 낮게 몰려 있는데, 이 환자(Index tumor) 는 저 멀리 오른쪽: 육종(Sarcoma) 샘플들과 비교해도 튄다. 오른쪽 작은 그림의 ACTB(housekeeping 유전자)는 대조로, 누구에게나 높아 분포 한가운데 위치한다: 이런 건 아웃라이어가 아니다.
튄다고 곧 원인은 아니다: 해석의 함정
섹션 제목: “튄다고 곧 원인은 아니다: 해석의 함정”발현량이 높게 튀는 걸 발견해도, 그 자체로 “이 유전자가 암을 일으켰다”고 말할 수는 없습니다. 높은 발현은 다음 중 무엇이든 될 수 있습니다:
| 성격 | 의미 |
|---|---|
| 원인(cause) | 종양 발생·진행에 직접 기여 |
| 결과(consequence) | 종양 관련 pathway가 켜지면서 부차적으로 함께 오른 발현 |
| 상태(state) | 종양 미세환경 변화(저산소·영양 부족·전이 능력 획득 등)에 적응하며 바뀐 발현 |
지금 데이터에서 확실히 아는 건 “발현이 튄다” 뿐입니다. 인과를 밝히려면 아래가 함께 필요합니다:
WGS 변이 분석 + pathway 분석 + 문헌조사: 예컨대 과발현이 이 유전자 자체의 변이 때문인지, 아니면 변이는 없는데 다른 조절 유전자의 영향으로 튀는 것인지를 가르고, 후자라면 그 조절 고리를 끊는 접근을 생각해 볼 수 있습니다. 여기에 다른 오믹스(예: 전사인자·효소가 붙을 수 있는 열린 자리를 보는 ATAC-seq)를 얹어 층층이 확인합니다.
그럼에도 아웃라이어 탐색이 유용한 이유는, 정상·타 환자 대비 변화가 분명한 유전자를 추려 표적 치료 후보·후보 약물 탐색의 출발점으로 삼을 수 있기 때문입니다: Sid 케이스가 MDM2 과발현을 표적으로 잡은 것처럼.
2. 비교 기준: 레퍼런스 데이터와 검증
섹션 제목: “2. 비교 기준: 레퍼런스 데이터와 검증”한 샘플을 “비교”하려면 기준이 될 공개 레퍼런스가 필요합니다. Sid가 쓴 것과 같은, 바로 받아 쓸 수 있는 TPM 수준의 두 데이터를 썼습니다:
| 레퍼런스 | 규모 | 성격 |
|---|---|---|
| GTEx | 74,628 genes · 19,616 샘플 (31 조직) | 정상 조직 발현 |
| PCAWG (EBI 경유) | 56,717 genes · 1,350 샘플 (27 암종) | 다른 암 환자 발현 |
세 데이터셋(내 샘플 + GTEx + PCAWG)에 모두 존재하는 공통 유전자 45,698개만 추려서 비교했습니다.
이 교집합은 테이블 조인의 키를 맞추는 단계입니다. 데이터셋마다 유전자 ID 버전과 포함 범위가 다르므로 이름만 보고 열을 이어 붙이면 다른 유전자를 비교하거나 결측치를 대량으로 만들 수 있습니다. 유전자 ID 체계, 주석 버전, 발현량 단위와 전처리를 먼저 맞춰야 합니다.
🧩 그림의 암 레퍼런스는 EBI(유럽 생물정보학 연구소) 의 Expression Atlas에서 받은 PCAWG 데이터이며, 그림·표에는
ebi라는 라벨로 표기했습니다. Sid 사례에 나온 PCAWG와 같은 계열입니다.
검증 ①: Boston vs Tempus 상관
섹션 제목: “검증 ①: Boston vs Tempus 상관”본론에 들어가기 전에, 만든 count가 믿을 만한지 먼저 봅니다. 같은 검체를 두 회사에 맡긴 데이터이므로, 둘의 발현량은 비슷한 분포여야 합니다.
같은 환자·같은 시점(T0)을 두 회사에서 시퀀싱한 결과의 log2(TPM+1) 산점도. 완전한 직선은 아니지만 대각선을 따라 강한 상관을 보여, 발현량 정량이 재현성 있게 만들어졌음을 확인. 상관과 로그 변환의 의미는 분포와 요약 통계에서 설명한다.
검증 ②: 이미 알려진 유전자가 실제로 튀는가
섹션 제목: “검증 ②: 이미 알려진 유전자가 실제로 튀는가”Sid 케이스에서 질병 관련 유전자로 명시됐던 것들이, 내 데이터에서도 같은 패턴으로 튀는지 봅니다.
한 플롯이 유전자 하나에 대응한다. 가로축은 log2(TPM+1), 세로축은 샘플 수. 파란색은 GTEx(정상), 초록색은 EBI(암) 분포이고, 주황 점선이 내 종양 값이다. FAP는 원래 싱글셀 데이터에서 찾아낸 표적이라 total RNA에서는 신호가 아주 도드라지진 않지만(그래도 과발현 쪽), 주요 표적으로 거론된 다른 유전자들은 분포를 한참 벗어나 찍힌다.
읽는 법: 대부분의 샘플은 히스토그램 중심 근처에 발현량이 몰려 있고, 내 샘플(주황 선)이 그 중심에서 얼마나 오른쪽으로 벗어났는지가 곧 “튀는 정도”입니다. 이건 Sid가 공개 뷰어에서 보여준 그림과 같은 형식입니다.
3. 아웃라이어를 고르는 지표
섹션 제목: “3. 아웃라이어를 고르는 지표”유전자 4.5만 개를 눈으로 다 볼 수는 없으니, 수치로 아웃라이어를 추리는 지표가 필요합니다. 종양 TPM을 레퍼런스 분포와 견주는 네 가지를 씁니다.
각 지표는 서로 다른 정렬 키입니다. 백분위는 분포 안의 위치, z-score는 평균으로부터의 표준화된 거리, fold change는 기준값과의 배수, max comparison은 관측 범위를 넘어섰는지를 나타냅니다. 한 열만 정답으로 두기보다 여러 열을 함께 보고 후보를 좁힙니다.
| 지표 | 무엇을 보나 | 기준 예시 |
|---|---|---|
| 백분위 순위(percentile) | 내 TPM이 레퍼런스 분포에서 상위 몇 % 인가 | max · 상위 1% · 5% 등 |
| Z-score | 레퍼런스 평균에서 표준편차 몇 배 떨어졌나 | ∣Z∣ ≥ 2 |
| Fold change | 레퍼런스 중앙값(median)의 몇 배 인가 | log2FC ≥ 2 → 4배 이상 |
| Max comparison | 정상에서 관찰된 최댓값보다도 종양이 높은가 | True / False (직접 비교) |
Z-score는 레퍼런스의 평균에서 표준편차 몇 칸 떨어졌는지 나타냅니다. Fold change는 기준보다 몇 배 큰지 나타내며, log2FC 2는 4배, log2FC 9는 512배입니다. 백분위·평균·중앙값·표준편차·로그 변환의 계산과 한계는 분포와 요약 통계를 참고하세요.
Max comparison은 논문에 널리 통용되는 표준 지표라기보다 Sid 사례를 재현하기 위한 직접 비교입니다. “정상에서 관찰된 최댓값까지 넘었는가”를 True 또는 False로 기록합니다.
실제 계산: 알려진 표적 유전자들
섹션 제목: “실제 계산: 알려진 표적 유전자들”암(EBI/PCAWG) 레퍼런스로 계산한 결과의 일부입니다. ebi_median·ebi_max는 암 환자 약 1,350명의 그 유전자 발현량의 중앙값·최댓값입니다.
| gene | tumor TPM | ebi median | ebi max | z-score | log2FC | ≈ 배수 |
|---|---|---|---|---|---|---|
| PANX3 | 619.43 | 0.138 | 19.84 | 9.14 | ≈ 560× | |
| MDM2 | 1235.72 | 4.954 | 2731 | 4.47 | 5.32 | ≈ 40× |
| FAP | 52.76 | 2.0 | 301 | 1.97 | 3.75 | ≈ 13× |
PANX3를 보면: 암 환자 중앙값이 0.138인데 내 종양은 619.43. 백분위로 최댓값마저 넘고, z-score는 19.84(평균에서 표준편차 20배 거리!), 중앙값 대비 약 560배(log2FC 9.14) 발현입니다. 교과서적인 아웃라이어죠.MDM2는 Sid가 실제 표적으로 삼은 유전자로, 종양 TPM 1235.72: 공개 뷰어에서는 1500 언저리로 표시되는 값입니다.FAP는 신호가 상대적으로 약합니다(z ≈ 2). total RNA에서는 도드라지지 않는, 싱글셀에서 잡힌 표적이라는 점과 일치합니다.
1차 필터: 정상/암 최댓값 대비 2배 선
섹션 제목: “1차 필터: 정상/암 최댓값 대비 2배 선”여러 지표를 조합해 후보를 좁힙니다. 예컨대 종양 TPM을 EBI 암 max TPM과 직접 산점도로 놓고, 대각선(y=x)과 그 2배 선(log2 기준 +1)을 그으면:
점 하나가 유전자 하나(27개 암종·1,350 샘플의 max 기준). 대각선(y=x)은 “암 최댓값과 같음”, 위쪽 점선은 그 2배 선. PANX3·MDM2처럼 2배 선을 넘어 위로 벗어난 유전자가 1차 아웃라이어 후보다. max 기준이라 종양에서만 튀는 완벽한 그림은 아니고(샘플 간 변이가 있어), 그래서 여러 지표를 함께 본다.
이렇게 레퍼런스 max·중앙값 대비 2배 이상 벗어난 유전자를 1차 아웃라이어로 추려 냅니다. 한 지표로 순위를 매길 수도, 여러 지표를 조합할 수도 있습니다.
한 샘플의 값은 혼자서는 높고 낮음을 말할 수 없습니다. 같은 유전자의 정상조직·다른 암 분포 위에 놓고 거리를 측정해야 합니다.
발현량 행렬 한 열 → 공통 유전자 선택 → 레퍼런스 분포와 비교 → 백분위·z-score·fold change로 후보 정렬
아웃라이어는 추가로 조사할 후보이지 암의 원인이나 치료 표적이라는 결론이 아닙니다. 다음 유전자 목록의 공통 기능 찾기에서는 후보 유전자 목록에 어떤 기능이 몰려 있는지 해석합니다.
- 정렬 도구: STAR · Samtools
- 정량: RSEM · featureCounts (Subread)
- 정상 조직 발현 레퍼런스: GTEx Portal
- 암 발현 레퍼런스: PCAWG (Nature 2020) · EBI Expression Atlas
- 기능 강화 분석: g:Profiler · NetworkAnalyst
- Sid 공개 발현 뷰어: osteosarc.com/rnaseq