analysis

앞으로 태어날 우리 아기, 데이터로 미리 만나볼까~

by

여기에서 노트북 형식으로도 볼 수 있습니다. 주피터 노트북에서 직접 써 보고 싶은 분은 여기에서 받으실 수 있습니다.

장혜식
http://highthroughput.org/

[알림 1]
이 글에서 소개하는 정보와 소스코드, 방법론 등은 모두 연구 또는 학술 목적으로만 사용 가능합니다. 의료용으로 사용하실 때는 의료기관을 통해 법률에서 정하는 절차로 사용하셔야 합니다.

[알림 2] 이 노트북은 유전학은 잘 모르지만 데이터 분석에 관심이 많은 분들과, 데이터 분석은 잘 모르지만 유전학은 잘 아는 분들도 큰 사전지식 없이 읽으실 수 있게 쓰려고 노력했습니다. 그래서 정확성을 희생하고 간단한 말로 대체한 부분이 많으니 제대로 된 정의에 대해서는 논문들을 참고하시길 부탁드립니다. 사실 저도 인간유전체학을 하진 않기 때문에 정확히는 잘 모릅니다. 하핫;;

[알림 3] 이 노트북에서 사용된 시퀀싱 데이터는 테라젠바이오연구소에서 만들어주셨습니다.

NIPT (non-invasive prenatal testing)는 임신 초기에 아기의 유전형을 검사하는 가장 최근에 나온 방법이다. 이름이 비침습적이라고 진짜 아무 곳도 안 찌르고 들어가는 건 아니고, 엄마 피를 뽑아서 분석한다. 그래도, 아기는 안 건드리기 때문에 양수 검사나 융모막 검사처럼 문제가 생겼을 때 아기에게 해가 가지는 않는다.

엄마 피로 도대체 아기의 유전자를 어떻게 검사하나? 놀랍게도(!) 임신 10-14주 정도 되는 엄마 피에 떠다니는 DNA의 대략 5-20% 정도는 아기 DNA라고 한다. 정확하게 말하면 엄마 피에 있는 DNA 중에 백혈구 같은 세포 안에 들어있는 DNA 말고, 둥둥 밖에 떠다니는 DNA 중의 5-20% 이다. 이 아기 DNA를 잘 분석하면 아기가 어떤 유전형을 가졌는지 아기를 건드리지 않고도 확인할 수 있다는 것이 NIPT의 기본 개념이다.

이론 상으로는 이런 방법으로 아기의 정말 자세한 유전형을 알아낼 수도 있지만, 엄마 DNA가 무려 90% 가까이 차지하고, 나머지 10% 가량만 아기 DNA라서 현실적으론 어렵다. 아기 DNA에 대한 커버리지만 충족시키는데도 기본적인 전체 유전체 시퀀싱(WGS)보다 10배 이상 드는데, 게다가 엄마DNA와 섞여있다보니 어느게 누구 것인지 통계적으로 확실히 구분할 수 있는 수준까지 가자면 일반적인 WGS의 수 백배가 필요하다. 다른 것 다 빼고 시퀀싱 시약 비용만 억대가 넘어간다.

그래서 지금 상용화되어 있는 NIPT에서는 다른 것은 모두 포기하고 염색체 개수가 이상한 이수성(aneuploidy)만 검사한다. 이수성으론 가장 유명한 것이 21번 염색체가 2개가 아니라 3개가 들어가는 3염색체성(trisomy 21)인데, 이 문제가 생기면 거의 대부분 다운 증후군이 나타난다. Trisomy 21은 생각보다 자주 발생하는 현상이라 완전 남 얘기가 아니다. 40대 초반 산모에서는 1%까지 빈도가 올라간다.

NIPT로 aneuploidy는 어떻게 검사할까? 현재 업계에서 쓰는 방법에선 0.1-0.5X 정도의 커버리지로 전체 유전체를 임의로 찔끔찔끔 읽어낸다. 즉, 엄마 피에 돌아다니는 DNA 조각 중에 아무거나 골라서 그게 뭔지 A/C/G/T가 연속되는 문자열로 읽어내는데(sequencing), 이걸 이미 알고 있는 사람 DNA 전체 문자열(서열)에 딱 맞는 부분을 찾아서 배열했을 때 대충 10%(0.1X)에서 50%(0.5X)정도 덮일 정도로 살짝 읽어내는 것이다. 많이 읽을 수록 커버리지가 올라가고 데이터도 더 좋아지지만, 돈이 많이 든다.

커버리지가 워낙 낮다보니 뭐 다른 건 분석할 수 있는 게 거의 없고, 사람 DNA 전체에서 일정 간격으로 자른 구간에서 DNA가 몇 번 읽혔는지 세서 그 숫자를 잘 비교하는 수 밖에 없다. 아기 DNA가 10%, 엄마 DNA가 90% 인 샘플이라면, 아기 DNA가 정상적인 염색체에 비해 trisomy가 있는 염색체는 5% 정도 더 읽히게 된다. 염색체가 1벌만 들어가는 monosomy가 발생했다면 5% 적게 읽힌다.

밤톨이의 등장

원래 넷을 낳겠다는 아내의 원대한 계획에 발 맞춰 얼마 전 둘째가 생겼다. 장모님이 밤이 후두둑 떨어지는 걸 받는 꿈을 꿨다고 해서 태명이 어쩌다가 밤톨이가 됐다. 원래 병원을 별로 좋아하지 않는 아내지만, 집 앞에 있는 병원의 시스템이 특별히 ESC를 누르지 않으면 저절로 OK가 계속 눌리는 프로그램처럼 딱히 얘기하지 않아도 자꾸 검사를 하게 돼 있어서 초기 검사를 꽤 했다. 12주엔 목둘레 투명대 검사(nuchal scan)를 했는데 한참을 재더니 재는 각도에 따라서 3.0mm이 넘기도 하고 아무리 좁게 잡아도 2.7mm은 항상 넘는다고 했다. 이 두께로 나오면 여러가지 염색체 이상이 발생했거나, 앞으로 신경관/심장 발달에 문제가 있을 가능성이 높다고 양수검사나 융모막검사를 하자고 했다.

그런데 이 검사들은 침습적 검사라서 실패하면 유산 가능성이 있고, 그 가능성이 1%를 넘기도 한다. 너무 위험한 검사라서 안 그래도 하루 종일 걱정을 하는 아내가 더욱 걱정을 하기 시작했다. 아내는 며칠동안 하루 종일 목둘레 검사 결과가 안 좋았지만 정상적으로 출산한 산모들과 양수검사를 하다가 사고가 난 엄마들의 인터넷 글들을 검색해서 읽었다. 정말 끝도 없이 나왔다. -O-; 여러모로 알아보다 기저율이 그렇게 높지 않은 나이이니까 문제가 없을 가능성이 훨씬 높다고 생각하고, 위험을 무릅쓸 필요는 없겠다 싶어서 그냥 추가 검사는 하지 않고 이전 검사 결과는 그냥 무시하기로 했다.

그러던 중에, 테라젠에서 NIPT검사에 관련된 논문을 냈다는 소식을 들었다. 예전에 같이 연구했던 친구들이 요즘 회사에서 NIPT사업화 관련해서 연구를 많이 하고 있다는 얘기도 듣고, 가장 활발하게 시퀀싱 사업을 하고 있는 Illumina에서도 요즘 NIPT 프로토콜과 제품을 많이 개발하고 있던게 생각났다. 그래서 NIPT 방법에 대해서 공부도 좀 해 보고, 인간 유전체 분석은 한 번도 안 해 봤기에, 그것도 어떤 방법으로 하는지 궁금해서, 테라젠의 김태형 이사님께 연락해서 아내의 NIPT 검사 데이터를 만들었다. 아직 태어나지도 않은 아이를 DNA 데이터로 먼저 만나다니! 정말 공부도 쏙쏙 되고, 뭔가 획기적인 결과를 하나 만들 수 있을 것만 같았다. ㅋㅋ

기술적으론, 아내의 피에서 뽑은 세포 밖 DNA를 Ion Proton으로 0.3X 정도 읽어서 각 염색체에 매핑되는 리드 개수를 몇 가지 보정한 후에 서로 다른 염색체들끼리의 리드 개수 비율을 구한다. 그 비율은 염색체별로 시퀀싱 효율 차이가 있기 때문에, 똑같은 DNA 준비 과정과 똑같은 시퀀싱 플랫폼을 거친 다른 임신부들의 검사 결과를 백그라운드 분포로 잡고 비교하는 방식으로 이수성의 가능성을 계산한다.

리드 매핑하고 bin에 들어가는 리드 세기

먼저 시퀀싱센터에서 받은 데이터를 매핑해서 지놈 전체를 100kbp 단위 빈으로 나눠서 리드를 센다. 빈 크기는 50k부터 300kbp까지 다양하게 많이 쓰이지만, 이번 데이터 경우에는 100kbp가 노이즈는 적으면서도 충분한 빈 개수를 확보할 수 있는 적당한 수준이었다. 사실 빈 크기가 결과에 큰 영향을 주지는 않는다.

시퀀싱 기계에서 읽은 DNA 서열을 이미 알려져 있는 사람 DNA 서열 어느 부분에 해당하는지 찾아보는 매핑 과정과 각 빈에 들어가는 리드 수를 세는 과정은 snakemake로 파이프라인을 간단하게 만드는 게 더 쉬우니까 이걸로 일단 돌린다. 매핑은 웬만한 프로그램과 비교해선 진짜 말 그대로 100배 이상 느린 것으로 유명한 GSNAP으로 한다. 보통은 BWA나 bowtie2 같은 것을 많이 사용하는데, 리드 개수도 얼마 안 되고 CPU와 메모리가 아주 펑펑 써도 충분할 정도로 남아 돌아서.. ㅎㅎ

우선 데이터가 있는 곳으로 뿅

In [2]:
%cd ~/tmp/niptraw

/atp/hyeshik/tmp/niptraw

아래는 Snakefile 얼개 그림과 내용. Snakefile 맨 윗쪽에서 쓰고 있는 powersnake 모듈은 여기에서 받을 수 있다.

In [31]:
!snakemake --rulegraph | dot -Tpng > rule.png

from IPython.display import Image
Image(filename='rule.png')
Out[31]:
In [32]:
!cat Snakefile

from powersnake import *

SAMPLES = ['bamtori']
BINWIDTH = 100000
GENOME_FASTA_URL = 'http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz'
GENOME_CHROMSIZES_URL = 'http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes'
SNP_TABLE_URL = 'http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp146.txt.gz'

rule all:
    input:
        expand('{sample}.sorted.bam{sufx}', sample=SAMPLES, sufx=['.bai', '']),
        expand('{sample}.bincov.txt', sample=SAMPLES),
        'snp146.bed.gz'

rule download_genome_reference_sequence:
    output: 'reference/hg38.fa.gz'
    shell: 'wget -O {output} {GENOME_FASTA_URL}'

rule download_genome_size_table:
    output: 'reference/hg38.chrom.sizes'
    shell: 'wget -O {output} {GENOME_CHROMSIZES_URL}'

rule build_gsnap_index:
    input: 'reference/hg38.fa.gz'
    output: 'reference/hg38/hg38.genomecomp'
    shell: 'gmap_build --gunzip --circular=chrM,chrMT  \
                -D reference -d hg38 {input}'

rule align_to_genome:
    input: ref='reference/hg38/hg38.genomecomp', fastq='{sample}.fastq.gz'
    output: temp('{sample}.bam')
    threads: 32
    shell: 'gsnap -D reference -d hg38 -m 0.05 -t {threads} -A sam -B4 \
                --gunzip {input.fastq} | \
            samtools view -b -q 3 -o {output} -'

rule sort_alignments:
    input: '{sample}.bam'
    output: '{sample}.sorted.bam'
    threads: 32
    shell: 'samtools sort -@ {threads} -o {output} {input}'

rule index_alignments:
    input: '{sample}.bam'
    output: '{sample}.bam.bai'
    shell: 'samtools index {input}'

rule make_depth_table:
    input: '{sample}.sorted.bam'
    output: '{sample}.bincov.txt'
    shell: 'bedtools makewindows -g reference/hg38.chrom.sizes -w {BINWIDTH} | \
            bedtools coverage -g reference/hg38.chrom.sizes -b {input} -a - \
                -nonamecheck > {output}'

rule download_snp_table:
    output: temp('snp146.txt.gz')
    shell: 'wget -O {output} {SNP_TABLE_URL}'

rule convert_snp_to_bed:
    input: 'snp146.txt.gz'
    output: 'snp146.bed.gz'
    shell: 'zcat {input} | cut -f2-7 | gzip -c - > {output}'

# vim: syntax=snakemake

레퍼런스를 준비해서 매핑하고 빈에 들어간 리드를 센다.

In [38]:
!snakemake -j -- 2>&1 > niptrun.log

리드 개수 데이터 훑어보고 전처리하기

작업 준비가 모두 끝났으니, 이상한 것이 없는지 한 번 살펴보고 쓰기 좋게 손질을 좀 해 둔다.

우선 테이블을 불러들인다. bedtools coverage에서 나오는 포맷이다.

In [33]:
coverage = pd.read_table('bamtori.bincov.txt',
                         names=['chrom', 'start', 'end', 'count', 'nzpos', 'length', 'nzratio'])

염색체 번호 순서대로 정렬할 수 있게 키를 만들어서 정렬한다. 크기가 균일하지 않은 맨 마지막 빈들은 제거한다. Y 염색체도 온전한 빈이 550개가 넘게 나오기 때문에 굳이 궁상맞게 따로 처리할 필요는 없다. 물리적으로 실재하는 염색체 말고 가상의 염색체들(alternative chromosome이나 unplaced scaffolds 등)은 제대로 하자면 복잡하니까 미리 제거한다.

In [35]:
chrom_human_order = lambda x: 'chr{:03d}'.format(int(x[3:])) if x[3:].isdigit() else x

binsize = coverage['length'].value_counts().index[0]

covs = coverage[(coverage['length'] == binsize) & (coverage['chrom'].apply(len) <= 6)].copy()
covs['chrom_sortkey'] = covs['chrom'].apply(chrom_human_order)
covs = covs.sort_values(by=['chrom_sortkey', 'start']).reset_index(drop=True)
covs.head()
Out[35]:
chrom start end count nzpos length nzratio chrom_sortkey
0 chr1 0 100000 175 18882 100000 0.18882 chr001
1 chr1 100000 200000 87 9354 100000 0.09354 chr001
2 chr1 200000 300000 59 6184 100000 0.06184 chr001
3 chr1 300000 400000 34 4553 100000 0.04553 chr001
4 chr1 400000 500000 41 5205 100000 0.05205 chr001

새 컬럼 chrom_sortkey는 chrom으로 문자열 정렬하면 보통 생각하는 숫자를 기준으로 한 염색체 순서를 맞추기가 어려워서 따로 0을 덧붙여서 만들었다.

컬럼 count에 각 빈에 들어가는 리드 개수가 들어 있다. nzpos나 nzratio는 이번엔 특별히 쓸 일은 없다.

이런 저런 개인적으로 자주 쓰는 라이브러리들을 쓰기 위해 패스를 추가한다. tailseeker 패키지는 GitHub 저장소에서 받을 수 있다.

In [37]:
import sys
sys.path.append('/atp/hyeshik/p/tailseeker')
import readline
from tailseeker import plotutils, stats
from matplotlib import style
style.use('ggplot-hs') # ggplot-hs 없으면 ggplot으로 바꾸세요

matplotlib 스타일 ggplot-hs는 내장 스타일인 ggplot에서 레이블 글꼴 크기, 색깔만 조금 바꾼 것이다. readline은 conda와 시스템 기본 libreadline이 버전이 다르다보니 rpy2 모듈을 부를 때 R과 Python이 서로 다른 libreadline을 부르려고 해서 어쩔 수 없이 미리 readline을 불러야만 rpy2를 쓸 수 있게 돼서 우선 이렇게 쓰고 있다.; 다음 업그레이드할 때는 고쳐야겠다.

일단 read count가 대충 어떤 분포로 나오는지 살펴보자. 요즘 생물학에선 바이올린 플랏이 대유행이다. 유행을 좇아 본다.

In [41]:
plt.figure(figsize=(9, 3))
vp = plt.violinplot([np.log2(bins['count'] + 1) for chrom, bins in covs.groupby('chrom')],
                    showextrema=False, showmedians=True, widths=0.8)
plt.setp(vp['bodies'], facecolor='k', alpha=.65)
plt.setp(vp['cmedians'], lw=2)

# chromosome 번호로 X 틱 레이블을 단다.
chromosomes = covs['chrom'].unique()
plt.xticks(np.arange(1, len(chromosomes) + 1),
           [c[3:] for c in chromosomes])

plt.xlabel('chromosome')
plt.ylabel('log2 read count/bin')
Out[41]:
<matplotlib.text.Text at 0x7f717c7155c0>

특별한 처리 없이도 대충 7.3 근처에서 모인다. 그런데! Y염색체에서 대충 중간값이 4.2 정도 되는 한 무리가 나온다. 이 수치를 다른 상염색체와 비교하면 대략 Y 염색체의 비율이 $2^{4.2-7.3}$ 쯤 된다고 할 수 있다. 즉, 11.7% 쯤 되는데, Y염색체는 엄마 DNA에는 없으므로 전적으로 아기 DNA라고 할 수 있다. 상염색체 기준으로 하면 23.4%에 육박하는 수치이므로 보통 아기 DNA가 10-20% 나오는 걸 감안하면 Y염색체가 있다는 걸 확실히 알 수 있다. IGV로 봐도 repeat 지역 여부에 관계없이 고루 퍼져있는 것을 볼 수 있었다. Y는 반복서열이 워낙 많아서 매핑이 어렵고 centromere의 비중이 상대적으로 높기 때문에 아래에 0 근처에 깔리는 빈이 많은 것은 예상할 수 있는 결과다. 아…. 귀염둥이 딸과 노는 꿈은 물거품이.. Y_Y 그래도 아들 둘이면 같이 잘 놀테니까! ^_^

분포만 그리면 재미가 없으니까 outlier가 대충 어떻게 분포하는지 멋지게 점으로 찍어보자. 염색체들을 이어 붙여서 찍을 것이므로 다른 염색체들끼리 구분할 수 있게 무지개색을 준비한다. Y염색체는 칙칙하니까 회색으로 칠한다.

In [44]:
chroms = sorted(covs['chrom_sortkey'].unique())
chromcolors = {chrom: color for chrom, color
               in zip(chroms, plotutils.colormap_lch(len(chroms), end=360))}
chromcolors['chrY'] = '#3d3d3d'

자 이제 찍는다. 각 염색체별로 median은 보기 좋게 바이올린 플랏과 마찬가지로 줄로 긋는다. 염색체 이름은 색깔 맞춰서 중간쯤에 그리는데, 뒤로 가면 글자끼리 겹치니까 위, 아래로 번갈아 표시한다.

In [48]:
fig, ax = plt.subplots(1, 1, figsize=(10, 4))

chromnameY = 8.5
wiggle_y = .2

for i, (chrom, rows) in enumerate(covs.groupby('chrom_sortkey')):
    ax.scatter(rows.index, np.log2(rows['count']),
               edgecolor='none', s=4, alpha=0.2,
               c=chromcolors[chrom])

    center_y = np.median(np.log2(rows['count']).dropna())
    ax.plot([rows.index[0], rows.index[-1]],
            [center_y, center_y],
            c='black')
    
    center_x = np.mean(rows.index)
    ax.annotate(chrom[3:].lstrip('0'),
                (center_x, chromnameY + (wiggle_y * (i % 2) * 2 - wiggle_y)),
                color=chromcolors[chrom], ha='center')

plotutils.apply_dropped_spine(ax)

ax.set_ylabel('Reads per 100k bin (log2)')
ax.set_xlim(-100, 15800*2)
ax.set_ylim(2, 10)
Out[48]:
(2, 10)

centromere 지역에서 집중적으로 빠지거나 튀는 부분이 눈에 띈다. 뭐 크게 문제는 아니다. 특별히 상염색체들끼리는 돋보이는 차이는 없고, 19번과 22번이 다른 것들보다 좀 떨어진다. 19번은 유전자 밀도가 워낙 높아서 GC ratio도 독보적으로 높고, 그 덕에 시퀀싱이 잘 안 되는 것으로 유명한다. 자세히 보지는 않았지만 아마도 그 원인이겠지? ; 22번도 유전자 밀도가 좀 높은 편이다.

사실 이수성(aneuploidy)이 있으면 log2 스케일로 0.15 정도 위아래로 왔다갔다 할 것이다. trisomy가 있다면 22번이 아래로 내려온 정도로 위로 튀면 된다. trisomy가 가장 잘 발생하는 21번 (다운 증후군), 18번 (에드워즈 증후군), 13번 (파타우 증후군) 염색체들은 다행히도 대충 대세를 따른다. 사실 이렇게 한 데이터만 놓고서 분석해서는 정확히 알기는 힘들고, 같은 플랫폼에서 검사한 다른 사람들의 데이터가 많이 모여야 정상치의 범위를 알 수 있는데.. 일단은 데이터가 1개 밖에 없으니 이걸로 대충 보고 만족한다.

GC 비율로 보정하기

사실 웬만한 것은 상용화될 수준까지 가면 간단하지 않은 문제가 돼 버린다. NIPT를 이용한 이수성 검사에서도 다른 사람들의 데이터와 비교하는 방식으로 통계적 유의성을 구하지만, 데이터를 좀 더 균일하게 만드는 여러가지 편향 보정법이 알려져있다.

가장 대표적인 것은 GC 비율 보정인데, 간단하다보니 대부분 업체에서 이 보정은 거의 반드시 채택하고 있다. DNA를 구성하는 AT/GC 중에 염색체에서 GC가 많이 들어있는 지역은 DNA 2가닥끼리 서로 붙어있으려는 힘이 훨씬 강하다. 유전체 전체를 시퀀싱할 때는 DNA를 그냥 쓰면 너무 길어서 다루기가 힘들어서, 초음파를 쏘는 등의 방법으로 DNA를 물리적으로 갈기갈기 찢은 다음에 읽는다. (여기서 쓰는 초음파는 산부인과에서 쓰는 초음파보다 훨씬 세다. 물 속에 담그고 쓰는데도 어쩌다 들으면 귀가 멍멍하다.) 그런데 GC가 많은 지역은 서로 붙어있으려는 힘이 강하다보니 그냥 흔드는 방법으로는 훨씬 덜 찢어져서, 시퀀싱까지 살려서 가기 어렵다. 그래서 GC가 많은 지역은 수가 적게 나올 수 밖에 없다.

그럼 우리도 한 번 해 보자! 먼저 사람 유전체에서 앞에서 계속 쓰던 100kbp 구간별로 GC 비율을 구한다. 이 비율은 굳이 귀찮게 프로그램 만들지 않아도 bedtools느님이 기본제공해 주신다. 하지만 랜덤액세스를 하기 때문에 gzip 압축된 상태로는 못하니까 일단 풀고 한다.

In [73]:
!zcat reference/hg38.fa.gz > reference/hg38.fa
In [57]:
genome_nuc_comp = !bedtools makewindows -g reference/hg38.chrom.sizes -w 100000 | bedtools nuc -fi reference/hg38.fa -bed -
In [58]:
import io
nuccomp = pd.read_table(io.StringIO('\n'.join(genome_nuc_comp))).iloc[:, :5]
nuccomp.columns = ['chrom', 'start', 'end', 'pct_at', 'pct_gc']
nuccomp.head()
Out[58]:
chrom start end pct_at pct_gc
0 chr1 0 100000 0.51793 0.38207
1 chr1 100000 200000 0.54185 0.45815
2 chr1 200000 300000 0.28508 0.19460
3 chr1 300000 400000 0.27479 0.24553
4 chr1 400000 500000 0.59808 0.40192

몇 개 빈에서 AT와 GC를 합쳐도 1이 안 되는 게 있다. 베이스를 모르는 N때문에 그런데, N은 AT인지 GC인지 모르니까 그냥 얘네들은 빼고 다시 계산한다.

In [59]:
nuccomp['pct_gc'] = nuccomp['pct_gc'] / (nuccomp['pct_gc'] + nuccomp['pct_at'])
del nuccomp['pct_at']

이제 앞에서 구해뒀던 빈별 커버리지 데이터와 GC비율을 합치고, 예상했던대로 GC비율과 리드 카운트가 관계가 있는지 보자.

In [66]:
covs_with_gc = pd.merge(covs, nuccomp, left_on=['chrom', 'start', 'end'],
                        right_on=['chrom', 'start', 'end']).dropna()
covs_with_gc.plot(kind='scatter', x='pct_gc', y='count',
                  edgecolor='none', s=5, c='black')
Out[66]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f7178ebc710>

위에 엄청나게 튄 outlier들때문에 뭘 볼 수가 없다. 이게 보이더라도 너무 겹쳐서 밀도가 잘 안 보일 것 같으니, 점의 밀도를 KDE로 색깔로 같이 찍어본다.

In [68]:
from scipy.stats import gaussian_kde
gc_count_mtx = covs_with_gc[['pct_gc', 'count']].as_matrix().T
density = gaussian_kde(gc_count_mtx, bw_method=0.1)(gc_count_mtx)
In [70]:
fig, axes = plt.subplots(1, 2, figsize=(10, 4))

covs_with_gc.plot(kind='scatter', x='pct_gc', y='count', edgecolor='none', s=5,
                  c=density, vmin=0, cmap='jet', ax=axes[0])
covs_with_gc.plot(kind='scatter', x='pct_gc', y='count', edgecolor='none', s=5,
                  c=density, vmin=0, cmap='jet', ax=axes[1])
axes[1].set_ylim(0, 500)
axes[1].set_xlim(0.3, 0.65)

plt.tight_layout()

오! GC 비율이 36-37%정도 되는 빈들이 가장 잘 나오고, 거기서 더 올라가면 떨어진다. 대략 GC 60%까지 가면 리드가 40% 정도 줄어드는 것 같다. AT가 적을 때도 약간 떨어지는 듯한 모양인데, 아마도 조각이 너무 잘아져서 크기로 정제하는 부분에서 없어지지 않을까 싶다. 이제 이걸로 LOESS regression을 하는 것이 생물학에서 보통 쓰는 레퍼토리인데.. GC 비율이 35%에서 42% 쯤 되는 지역의 아랫쪽에 아웃라이어가 꽤 있어서 약간 울퉁불퉁하게 나오기가 쉬울 것 같다. 아웃라이어들이 어떤 놈들인지 찾아봐서 제거하면 좋겠지만, 주업이 아니므로 대충 그냥 마구 잘라서 바로 LOESS에 집어넣는다. 아마도 centromere, heterochromatin이거나 reference와 다른 부분 등이 아닐까 싶다.

In [118]:
from sklearn import svm

gc_count_mtx_training = gc_count_mtx[:, gc_count_mtx[1] > 90]
model = svm.OneClassSVM(nu=0.06, kernel="rbf", gamma=0.5).fit(gc_count_mtx_training.T)
isinlier = model.predict(gc_count_mtx.T) == 1

gc_count_inliers = gc_count_mtx[:, isinlier]
fit_x, fit_y = stats.loess_fit(gc_count_inliers[0], gc_count_inliers[1],
                               np.linspace(0, 1, 200), alpha=.4)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

covs_with_gc.plot(kind='scatter', x='pct_gc', y='count', edgecolor='none', s=5,
                  c=isinlier, vmin=-.1, vmax=1.1, cmap='RdBu_r', ax=axes[0])
covs_with_gc.plot(kind='scatter', x='pct_gc', y='count', edgecolor='none', s=5,
                  c=isinlier, vmin=-.1, vmax=1.1, cmap='RdBu_r', ax=axes[1])
axes[1].set_ylim(0, 500)
axes[1].set_xlim(0.3, 0.65)
axes[1].plot(fit_x, fit_y, c='magenta', lw=3)
axes[0].plot(fit_x, fit_y, c='magenta', lw=3)

plt.tight_layout()

빨간색이 one class SVM으로 잡은 인라이어. 아래에 깔리는 아웃라이어들은 모델 파라미터에 넣어서 없애려니 금방 안 돼서 그냥 입력값을 미리 잘라서 넣어버렸다 -ㅇ-; 가운데 심지지역의 모서리 부분은 어차피 밀도가 중간에 비해서 많이 낮은 편이라 큰 영향은 안 주리라 보고 그냥 이대로 간다. 하하하하. =3 마젠타색 선이 LOESS 중심선인데, 대충 모양에 맞게 나왔다.

이제 원래 데이터를 LOESS 중심선에 맞춰서 끌어올린 다음에 원래 데이터와 보정한 데이터를 비교해 본다. log가 다루기 편하니, 이제 중심선에 비교해서 얼마나 enrich되었는지 log로 계산해서 쓴다.

In [125]:
gc_inliers_center = np.vstack(stats.loess_fit(gc_count_inliers[0], gc_count_inliers[1],
                                              gc_count_inliers[0], alpha=.4))

gc_logenr = gc_count_inliers.copy()
gc_logenr[1] = np.log2(gc_count_inliers[1] / gc_inliers_center[1])

orig_density = gaussian_kde(gc_count_inliers, bw_method=0.1)(gc_count_inliers)
adj_density = gaussian_kde(gc_logenr, bw_method=0.1)(gc_logenr)
In [129]:
fig, axes = plt.subplots(1, 2, figsize=(10, 4))

axes[0].scatter(gc_count_inliers[0], gc_count_inliers[1], c=orig_density, edgecolor='none', s=5)
axes[0].set_xlabel('pct_gc')
axes[0].set_ylabel('mapped reads')

axes[1].scatter(gc_logenr[0], gc_logenr[1], c=adj_density, edgecolor='none', s=5)
axes[1].set_xlabel('pct_gc')
axes[1].set_ylabel('adjusted mapped reads')
axes[1].set_ylim(-1, 1)

plt.tight_layout()

왼쪽은 GC비율로 보정하기 전, 오른쪽은 보정 후. GC에 맞춰서 파도치던 분포가 올곧게 딱 잡혔다. 후후훗.

이제 원래 작업하던 테이블에 보정된 값을 넣어둔다.

In [148]:
_, expected = stats.loess_fit(gc_count_inliers[0], gc_count_inliers[1],
                              covs_with_gc['pct_gc'], alpha=.4)
covs_w_logenr = covs_with_gc.copy()
covs_w_logenr['log2enrichment'] = np.log2(covs_w_logenr['count'] / expected)
covs_w_logenr = covs_w_logenr[covs_w_logenr['count'] > 0].dropna()

GC보정된 리드 밀도 살펴보기

혹시나 하는 마음에 Y염색체를 다시 살며시 들여다 본다.

In [174]:
chrY_enrichment = covs_w_logenr[covs_w_logenr['chrom'] == 'chrY']['log2enrichment']
xv = np.linspace(-8, 4, 100)
yv = gaussian_kde(chrY_enrichment, bw_method=0.2)(xv)
plt.plot(xv, yv)
plt.xlabel('log2 enrichment')
plt.ylabel('density')
plt.title('chromosome Y')
xv[np.argmax(yv)]
Out[174]:
-3.1515151515151514

대략 -3.15 정도에서 피크가 나온다. $2^{-3.15}=0.1127$ 로.. 여전히 아들이 확실해보인다.

앞에서 그렸던 빈 단위 커버리지를 이번에 GC보정한 log2 enrichment로 그리면 다른 것도 좀 더 깔끔히 나올까?

In [180]:
fig, ax = plt.subplots(1, 1, figsize=(10, 5))

chromnameY = -1.2
wiggle_y = .05

for i, (chrom, rows) in enumerate(covs_w_logenr.groupby('chrom_sortkey')):
    ax.scatter(rows.index, rows['log2enrichment'],
               edgecolor='none', s=4, alpha=0.4,
               c=chromcolors[chrom])

    center_y = np.median(rows['log2enrichment'].dropna())
    ax.plot([rows.index[0], rows.index[-1]],
            [center_y, center_y],
            c='black')
    
    center_x = np.mean(rows.index)
    ax.annotate(chrom[3:].lstrip('0'),
                (center_x, chromnameY + (wiggle_y * (i % 2) * 2 - wiggle_y)),
                color=chromcolors[chrom], ha='center')

plotutils.apply_dropped_spine(ax)

ax.set_ylabel('Reads per 100k bin (log2)')
ax.set_xlim(-100, 31000)
ax.set_ylim(-4, 1)
Out[180]:
(-4, 1)

오! 아까 GC보정 전엔 22번 염색체도 꽤 튀었는데, 이제 괜찮다. 19번 염색체도 훨씬 덜 튀어보인다. X 염색체는 GC보정 전에는 특별히 다른 상염색체와 달라보이지 않았는데, 이제 떨어진 게 확실히 보인다. 아기 DNA에 Y 때문에 X가 한 카피 적게 들어있으므로 대략 5-10% 정도 (log2로 0.07 정도) 빠지는 게 정상이다. 아들이 확실하다.;

그나저나 19번 염색체가 튀는 건 어떤 이유 때문인지 유전자가 없는 intergenic region만 골라서 해 본다던지 좀 더 자세히 들여다 봐야 할 것 같다.

엄마와 아빠 유전형으로 아기 DNA 더 세밀하게 구분하기

여기까지는 그냥 각 염색체별 리드 개수만 열심히 세고, 서열에서 나오는 상세한 정보는 따로 활용하지는 않았다. 시퀀싱의 멋은 그냥 세는 게 아니라 서열 내용을 자세히 분석하는 것 아닌가!
엄마 아빠의 유전형을 알면 읽어낸 서열로 더 자세한 분석이 가능한데, 아기 DNA는 엄마 DNA 반, 아빠 DNA 반을 물려 받으니까, 엄마와 아빠가 다른 유전형을 갖고 있는 부분을 봐서, 아빠와 유전형이 맞는 것을 고르면 아기 DNA의 반에 대해 거의 정확한 비율을 추산할 수 있게 된다.

상용 NIPT에서는 비용 문제 때문에 엄마와 아빠의 유전형 정보를 쓰지는 않지만, 우리 부부는 유전형 정보를 미리 갖고 있었기 때문에, 우리의 유전형 정보로 아기 DNA를 좀 더 세분화 해보기로 했다. 23andMe 데이터를 이런 곳에 처음으로 유용하게 써 볼 줄은 몰랐네..;

미리 준비를 좀 해야 한다. 23andMe는 GRCh37 기준 좌표계를 쓰는데, 우리는 여기서 GRCh38을 쓰고 있어서 좌표 변환이 필요하다. 굳이 번거롭게 변환하지 말고 UCSC Genome Browser가 정리해 놓은 dbSNP 데이터에서 GRCh38 좌표를 바로 얻어온다. dbSNP는 메모리에 띄워놓고 있기에 부담스러울 정도로 방대하니, 23andMe에 있는 것만 추린다.

In [182]:
snptbl = pd.read_table('snp146.bed.gz', compression='gzip',
                       names=['chrom', 'start', 'end', 'name', 'score', 'strand'])
snpcalls_m = pd.read_table('genotype-mother.txt.gz', compression='gzip',
                       comment='#', names=['name', 'chrom', 'position', 'genotype'],
                       low_memory=False)
snpcalls_f = pd.read_table('genotype-father.txt.gz', compression='gzip',
                       comment='#', names=['name', 'chrom', 'position', 'genotype'],
                       low_memory=False)
snpcalls = pd.merge(
    snpcalls_m[['name', 'genotype']], snpcalls_f[['name', 'genotype']],
    how='inner', left_on='name', right_on='name', suffixes=['_mo', '_fa'])
snpmerged = pd.merge(snptbl, snpcalls, how='inner', left_on='name', right_on='name')
snpmerged.to_csv('snpcalls.bed.gz', compression='gzip', header=False, index=False, sep='\t')
del snptbl, snpcalls_m, snpcalls_f, snpcalls

23andMe에서 제공하는 SNP으로 제한한 bed가 준비됐다. 지놈에 정렬된 서열들에서 이 SNP loci에 정렬된 녀석들을 골라온다.

In [184]:
!samtools mpileup -l snpcalls.txt.gz bamtori.sorted.bam | cut -f1,2,4,5,6 | gzip -c - > snpmatched.txt.gz

[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
In [186]:
serum_snps = pd.read_table('snpmatched.txt.gz', compression='gzip',
                          names=['chrom', 'position', 'reads', 'basecalls', 'basequalities'])
serum_snps['position'] -= 1

parental_snps =  pd.read_table('snpcalls.bed.gz', compression='gzip',
                          names=['chrom', 'start', 'end', 'name', 'score', 'strand',
                                 'mother', 'father'])

snps = pd.merge(parental_snps, serum_snps,
                how='inner', right_on=['chrom', 'position'],
                left_on=['chrom', 'start']).dropna()
snps.head()
Out[186]:
chrom start end name score strand mother father position reads basecalls basequalities
0 chr1 841165 841166 rs12124819 0 + AA AA 841165 1 A @
1 chr1 864626 864627 rs6681049 0 + CT CT 864626 1 T <
2 chr1 989517 989518 rs6665000 0 + AA AA 989517 1 a :
3 chr1 1095184 1095185 rs6687776 0 + CC CC 1095184 1 C-1N 7
4 chr1 1127257 1127258 rs9442373 0 + AC AC 1127257 1 A <

워낙 기본적인 리드 커버리지가 낮기 때문에 SNP별로 뭘 해 보기는 어렵고, 그냥 일치되는 게 있으면 모두 독립적인 SNP인 것처럼 따로 데이터를 취합한다. 엄마와 아빠가 유전형이 같은 곳에서는 아무 정보도 없으니까 엄마와 아빠가 유전형이 다른 곳만 일단 쓴다. 다른 곳은 이렇게 세 가지 유형이 있을 수 있다.

| 번호 | 엄마    | 아빠    |   예     |
|-----|--------|---------|---------|
|  1  | homo   | homo    | GG / AA |
|  2  | homo   | hetero  | GG / AG |
|  3  | hetero | homo    | AG / GG |

2번과 3번의 경우엔 확률 계산이 복잡해지니까 깨끗하게 떨어지는 1번을 우선 쓰자. 양쪽이 모두 homozygote이고 다른 타입이므로 직접적으로 비율을 계산할 수 있다. 아기는 무조건 heterozygote가 된다.

In [193]:
import re

called_homozygote = lambda gt: gt[0] == gt[1] if len(gt) == 2 and gt[0] != '-' else False
def classify_homozygotic_parental_snps(basereads, maternal, paternal,
                                       pat=re.compile('([+-][0-9]+)?([ACGTNacgtn])')):
    r = [0, 0, 0] # maternal matched, paternal matched, matched to neither of them
    for indelmark, base in pat.findall(basereads):
        if indelmark or base in 'nN':
            pass # skip it to simplify the analysis.

        ubase = base.upper()
        if ubase == maternal:
            r[0] += 1
        elif ubase == paternal:
            r[1] += 1
        else:
            r[2] += 1
    return r

snps_mofa_diffs = snps[snps['mother'] != snps['father']].copy()
snps_homo_homo = snps_mofa_diffs[
    snps_mofa_diffs['mother'].apply(called_homozygote) &
    snps_mofa_diffs['father'].apply(called_homozygote)]

snpcounts_homo_homo = pd.DataFrame(
    snps_homo_homo.apply(lambda row:
                         classify_homozygotic_parental_snps(row['basecalls'], row['mother'][0],
                                                            row['father'][0]), axis=1).tolist(),
    index=snps_homo_homo.index,
    columns='mother father mismatched'.split())
In [196]:
snpcounts_homo_homo.sum(axis=0)
Out[196]:
mother        10698
father          743
mismatched      113
dtype: int64

엄마와 아빠가 서로 다른 타입으로 homozygote인 SNP중에 엄마와 같은 것이 10698개, 아빠와 같은 것이 743개다. 아빠에 맞는 것이 아기DNA의 대략 반이므로, 아빠에 맞는 것의 2배를 아기 DNA 비율로 추산할 수 있다.

In [198]:
743/(743+10698)*100*2
Out[198]:
12.98837514203304

아기 DNA의 양은 좀 더 정확한 추산치는 이제 대략 13% 정도다. 오오.

그리고, 여기서 하나 더 알 수 있는 것은, 엄마와 아빠의 대립 SNP 중에 엄마 SNP에 압도적으로 많이 매치되는 것을 보아, 다른 사람의 피나 다른 사람의 샘플이 섞여서 엉뚱하게 나온 게 아니란 것이다.

그렇다면, 엄마에도 매치되지 않고 아빠에도 매치되지 않는 113개는 무엇일까? 거의 대부분 SNP은 A/G나 C/G 같이 두 염기가 대립하지 세 가지가 나오는 것은 드물다. 저것들은 아마도 시퀀싱에러가 아닐까 싶은데.. 확실히 해 보기 위해서 오독했을 가능성이 높은지 확인해 보자. 오독률을 알 수 있는 “리드 퀄리티”는 phred score로 보통 나타낸다. 10점이면 90% 정도 확률로 정확함, 20점이면 99% 확률로 정확함, 30점이면 99.9% 확률로 정확함을 나타낸다. 모두 기계에서 주는 추정값이다.

In [201]:
def iter_snpcall_reads(tbl):
    parse_basecalls = re.compile('([+-][0-9]+)?([ACGTNacgtn])')
    for _, row in tbl.iterrows():
        maternaltype = row['mother'][0]
        paternaltype = row['father'][0]
        basecalls = parse_basecalls.findall(row['basecalls'])
        for (indel, baseread), phredqual in zip(basecalls, row['basequalities']):
            if not indel and baseread not in 'nN':
                yield row['chrom'], row['position'], baseread.upper(), ord(phredqual) - 33, maternaltype, paternaltype

homosnpreads = pd.DataFrame.from_records(iter_snpcall_reads(snps_homo_homo),
                                         columns=['chrom', 'position', 'baseread', 'quality', 'maternal', 'paternal'])
assert (homosnpreads['baseread'] == homosnpreads['maternal']).sum() == 10665
In [208]:
maternal_reads = \
    homosnpreads[(homosnpreads['baseread'] == homosnpreads['maternal'])]
paternal_reads = \
    homosnpreads[(homosnpreads['baseread'] == homosnpreads['paternal'])]
mismatched_reads = \
    homosnpreads[(homosnpreads['baseread'] != homosnpreads['maternal']) &
                 (homosnpreads['baseread'] != homosnpreads['paternal'])]

mat_x, mat_y = plotutils.prepare_cumulative(maternal_reads['quality'])
pat_x, pat_y = plotutils.prepare_cumulative(paternal_reads['quality'])
mismatched_x, mismatched_y = plotutils.prepare_cumulative(mismatched_reads['quality'])

plt.plot(mat_x, mat_y, label='Maternal-matching loci')
plt.plot(pat_x, pat_y, label='Paternal-matching loci')
plt.plot(mismatched_x, mismatched_y, label='Mismatched loci', ls='--', c='magenta')
plt.legend(loc='best', fontsize=12)
Out[208]:
<matplotlib.legend.Legend at 0x7f7031098a20>

오.. 예상했던대로 엄마 것도 아니고 아빠 것도 아닌 SNP들은 리드 퀄리티가 현저히 떨어진다. 중간값이 대략 20쯤 된다.

자 그렇다면 염색체별로 차이가 있을 수 있으니 아빠:엄마 비율을 자세히 살펴보자.

In [211]:
r = []
for chrom, grp in homosnpreads.groupby('chrom'):
    if len(chrom) > 5:
        continue
    maternal_reads =  grp[(grp['baseread'] == grp['maternal'])]
    paternal_reads = grp[(grp['baseread'] == grp['paternal'])]
    mismatched_reads = grp[(grp['baseread'] != grp['maternal']) &
                           (grp['baseread'] != grp['paternal'])]
    r.append((chrom, len(maternal_reads), len(paternal_reads), len(mismatched_reads)))
    
homosnpmatches = pd.DataFrame(r, columns=['chrom', 'maternal', 'paternal', 'mismatched'])
homosnpmatches['p2m_ratio'] = homosnpmatches['paternal'] / (homosnpmatches['maternal'] + homosnpmatches['paternal'])
homosnpmatches.sort_values(by='p2m_ratio')
Out[211]:
chrom maternal paternal mismatched p2m_ratio
22 chrX 2 0 0 0.000000
14 chr22 186 7 0 0.036269
6 chr15 343 16 6 0.044568
19 chr7 602 29 2 0.045959
16 chr4 606 31 2 0.048666
4 chr13 390 23 2 0.055690
10 chr19 164 10 0 0.057471
15 chr3 758 48 3 0.059553
11 chr2 931 60 2 0.060545
20 chr8 620 40 2 0.060606
3 chr12 538 35 3 0.061082
5 chr14 290 19 1 0.061489
0 chr1 851 57 3 0.062775
1 chr10 615 44 11 0.066768
18 chr6 652 49 4 0.069900
21 chr9 537 43 1 0.074138
12 chr20 244 20 2 0.075758
2 chr11 571 47 4 0.076052
9 chr18 404 34 0 0.077626
8 chr17 229 20 2 0.080321
17 chr5 664 59 3 0.081604
7 chr16 278 25 2 0.082508
13 chr21 157 20 1 0.112994

X 염색체는 아빠가 1쪽 밖에 없어서 그런지 여기서 비교할 수 있는 것 자체가 적게 나왔다. (총 2개) 대부분 전체 유전체 평균 수치인 6.5% 근처로 나오는데, 22번 염색체가 유난히 적고, 21번 염색체가 유난히 많다. 21번 염색체는 수가 20개 밖에 되지 않지만, 아빠 SNP이 다른 염색체의 거의 2배가 나왔다. 이항분포로 보면,

In [222]:
%load_ext rpy2.ipython
%R pbinom(20, 157+20, 0.065, lower.tail=FALSE)

The rpy2.ipython extension is already loaded. To reload it, use:
  %reload_ext rpy2.ipython
Out[222]:
array([ 0.00574275])

흔히 나올 만한 숫자는 아니다. 21번 염색체가 더 많이 들어간 걸까? 좀 더 숫자를 늘리기 위해서 엄마는 homozygote, 아빠는 heterozygote인 경우도 해 보자.

In [226]:
snps_homo_hetero = snps_mofa_diffs[
    snps_mofa_diffs['mother'].apply(called_homozygote) &
    ~snps_mofa_diffs['father'].apply(called_homozygote)]

def iter_snpcall_reads_father_hetero(tbl):
    parse_basecalls = re.compile('([+-][0-9]+)?([ACGTNacgtn])')
    for _, row in tbl.iterrows():
        maternaltype = row['mother'][0]
        paternaltype = ''.join(set(row['father']) - set(maternaltype))
        if maternaltype == '-' or paternaltype == '-':
            continue
        basecalls = parse_basecalls.findall(row['basecalls'])
        for (indel, baseread), phredqual in zip(basecalls, row['basequalities']):
            if not indel and baseread not in 'nN':
                yield row['chrom'], row['position'], baseread.upper(), ord(phredqual) - 33, maternaltype, paternaltype

snpreads_fhetero = pd.DataFrame.from_records(iter_snpcall_reads_father_hetero(snps_homo_hetero),
                                         columns=['chrom', 'position', 'baseread', 'quality', 'maternal', 'paternal'])
In [227]:
r = []
for chrom, grp in snpreads_fhetero.groupby('chrom'):
    if len(chrom) > 5:
        continue
    maternal_reads =  grp[(grp['baseread'] == grp['maternal'])]
    paternal_reads = grp[(grp['baseread'] == grp['paternal'])]
    mismatched_reads = grp[(grp['baseread'] != grp['maternal']) &
                           (grp['baseread'] != grp['paternal'])]
    r.append((chrom, len(maternal_reads), len(paternal_reads), len(mismatched_reads)))
    
snpmatches_fhetero = pd.DataFrame(r, columns=['chrom', 'maternal', 'paternal', 'mismatched'])
snpmatches_fhetero['p2m_ratio'] = snpmatches_fhetero['paternal'] / (snpmatches_fhetero['maternal'] + snpmatches_fhetero['paternal'])
snpmatches_fhetero.sort_values(by='p2m_ratio')
Out[227]:
chrom maternal paternal mismatched p2m_ratio
22 chrX 4258 3 52 0.000704
9 chr18 1026 26 3 0.024715
17 chr5 2082 60 11 0.028011
21 chr9 1626 48 6 0.028674
4 chr13 1341 43 4 0.031069
0 chr1 2912 95 18 0.031593
1 chr10 2033 68 24 0.032366
15 chr3 2456 85 19 0.033451
5 chr14 1172 41 5 0.033800
11 chr2 2745 97 13 0.034131
18 chr6 2105 76 15 0.034846
12 chr20 871 32 4 0.035437
19 chr7 1912 71 10 0.035804
20 chr8 1853 75 11 0.038900
3 chr12 1660 68 7 0.039352
16 chr4 2143 88 9 0.039444
10 chr19 509 21 5 0.039623
7 chr16 966 40 9 0.039761
14 chr22 408 17 5 0.040000
2 chr11 1764 75 14 0.040783
8 chr17 877 40 8 0.043621
6 chr15 1137 52 12 0.043734
13 chr21 427 22 2 0.048998

오. 여기서 본 SNP 중에 X염색체에 있는 4000개가 넘는 SNP가 거의 몽땅 엄마에 매치된다. 역시 아들이다. ㅎㅎㅎ; (딸이면 아빠에게서도 X염색체를 물려받아서, 이 비율이 다른 상염색체와 거의 같게 나와야 한다.)

자 여기서도 21번이 다른 것들보다 유난히 높다. 아빠가 hetero이고 엄마가 homo이기 떄문에, 아기 DNA의 4분의 1은 아빠와 엄마가 다른 SNP loci에서 엄마와 맞지 않고 아빠와 맞게 된다. 따라서, 이 비율에 2배 하면 아빠한테서 물려받은 아기 DNA 비율을 추산할 수 있다.

21번 염색체의 수치인 4.9%는 아빠한테서 9.8% 정도를 물려받았다는 것이 된다. 앞에서 계산한 아기 DNA의 비율이 13% 정도 되기 때문에, 21번 염색체의 경우 거의 75%를 아빠한테서 물려받은 게 된다. 리드 수를 이용한 파이프라인에서는 trisomy 21 가능성은 거의 없다고 했다. 앞에서 homo-homo SNP들의 경우엔 다른 염색체들의 2배 가까이 아빠에 매치되었으니까, 종합하면 21번 염색체 두 카피를 모두 아빠한테서 받고 엄마한테서는 받지 않은 paternal uniparental disomy (patUPD)의 가능성도 무시할 수 없다.

하지만 논문을 열심히 뒤져보니, 21번 염색체의 patUPD가 워낙 드문 현상인데다, 21번 염색체는 다행히도 아직 UPD에 의한 표현형(즉 눈에 띄는 질병 또는 특징)이 보고된 것이 없다. 완전 정상이라고 보고된 사례만 몇 있다. 그래서 21번 염색체 patUPD도 가능성이 없는 것은 아니지만, 기저율 덕에 실질적인 가능성은 낮을 것 같고. 실제로 patUPD라고 하더라도 특별히 문제있을 것 같지 않다.

SNP 데이터로 뭘 좀 더 해보면 좋을텐데, 아내와 내가 다른 유전형이면서, 재미있는 정보가 알려진 것이 있고, 리드가 그래도 2개는 되는 걸 찾아보니 하나도 없었다. 커버리지가 이수성을 보는 것에 최적화되어 있다보니 어쩔 수 없는 것 같다. 몇 년 지나서 시퀀싱이 아주 싸지게 되면 전가족 시퀀싱을 한 번 해 보면 재미있겠다~

결론

여러가지 재미있게 분석을 해서 얻은 결론은

  • 아들임이 확실하다는 독립적인 증거를 셋이나 찾음. (Y염색체 비율, X염색체 비율 감소, X염색체 SNP)
  • 다른 사람 샘플과 뒤바뀌었을 가능성은 거의 없음.
  • 21번 염색체 patUPD가 있을 수도 있지만 가능성은 매우 낮고, 있더라도 아기에 큰 영향은 없을 듯하다.
  • 커버리지 0.3X는 생각보다 더 리드 구경하기가 힘들다. 대부분 SNP은 한 번도 읽히지 않고, 읽힌 것도 대부분 딱 1번 읽혀서, SNP단위로 뭘 해 보기는 어렵다.
  • 이수성만 검사하는 경우라면 커버리지를 훨씬 더 내려서 가격을 내릴 수 있을 것 같다.

말로만 듣던 NIPT 데이터를 대충이나마 직접 한 번 다뤄보니 감촉이 약간 느껴지는 것 같기도 하고 어떤 것인지 이해가 좀 됐다. 아직 눈으로 보지도 못한 둘째를 데이터로 미리 만나보니 그것도 기분이 참 좋았다. ^^; 셋째가 나오기 전에 이번 데이터로 뭔가 새로운 방법이나 개선점을 하나 만들 수 있으면 좋을텐데~

터미널엔 롯데리아, 동네에 맥도날드, 강남에 버거킹..?

by

2015년 1월 27일 개정 1판 안내

처음 공개된 분석과 글에 있었던 여러 오류를 바로 잡았습니다. 다음 문제점들을 수정하였습니다. 코멘트 주신 분들께 감사드립니다. ^^

  1. 사각형 cartogram 지역 재배치
    1. 서울시 노원구, 도봉구 위치 바로 잡음. (오유 자리****범 님)
    2. 강원도 춘천, 홍천 위치 바로 잡음. (오유 모* 님)
    3. 서울시 구로구, 금천구 위치 바로 잡음.
  2. 1월 26일 기준 데이터로 업데이트. (맥도날드와 롯데리아 매장이 많이 추가되었습니다.)
  3. KFC 석촌역점이 처리 도중에 누락됐던 문제 바로 잡음.
  4. 상관계수 분석을 잘못 뭉뚱그려 설명한 부분을 수정하여 구체적인 것으로 교체. (트위터 for*********er 님)
  5. 인구밀도-버거지수 스캐터 차트에서 분포가 좀 더 잘 보이게 LOESS 트렌드선 추가.
  6. 그림에서 일관성이 떨어지거나 부정확한 용어들이 섞여있던 것 수정.

그리고, 한 가지 덧붙이고 싶은 말은, 이 글은 진지한 어조와 형식으로 되어 있긴 하지만, 유머 또는 광고(IPython/pandas) 분류를 의도하고 쓴 글입니다. 결과에 대해 경영/사회 측면에서 심각하게 고민하실 정도의 정밀한 연구가 아니니, 구체적인 것은 해당 분야 전문가들이 제대로 분석하신 뒤에 보셔야 합니다. ^_^;;

읍내 버스터미널에 가면 롯데리아가 언제나 반갑게 맞이한다. 그렇게 보기 쉽던 롯데리아가 희한하게도 강남이나 신촌 같은 곳에서는 뜻밖에 찾아보기 힘들다. 작년 가을, 트위터 @RioterOfMiku 님은 이 현상을 한 눈에 이해되는 아주 강렬한 도시 발전 지수로 정리했다.

https://twitter.com/godtsune_miku/status/513648274406789120

누구든 대략 듣자마자 “그래? 그런 것 같네~” 정도론 동의하게 되는 간단하면서도 정보가 있는 지수인데. 어느 날 연구실에서 오랜만에 햄버거를 시켜 먹다가 문득 생각나서 “진짜 그럴까!?” 생각에 좀 더 구석구석 뒤져보기로 했다.

구체적인 것은 분석 과정을 담은 IPython 노트북을 올려두었으니 관심있는 분들은 여기를 좀 더 보도록 합시다. ^_^

전국 시군구 버거지수

우선, 전국 시군구에 대한 버거지수를 진짜로 계산해 본 적은 없으니 실제로 계산해 보자. 각각 롯데리아, 버거킹, 맥도날드, KFC 홈페이지에서 2015년 1월 26일 기준 데이터를 긁어서 계산했다.

정말 느끼는 대로 강남, 서초는 하늘을 찌르는 “버거지수”를 보여준다. 즉, 강남과 서초에는 롯데리아에 비해 버거킹, 맥도날드, KFC가 훨씬 많다. 다른 지역들은 이 비율이 대략 0.5를 넘는 곳도 드문 와중에, 4배를 넘어선다. 다른 지역에서도 주로 “시내”권인 중구 지역이 버거지수가 좀 더 높다. 하나 의외인 곳은 강원도 홍천군이 2.0점! 여기는 강남, 서초 다음으로 전국에서 3번째로 높은데… 매장 목록을 보면 대명비발디파크에 있는 버거킹 2개 때문이다.; 강원도 정선군은 색은 상당히 짙은 게 높아 보이지만, 롯데리아와 KFC가 1개 씩으로 롯데리아 수가 적어서 지수가 높게 나온 영향이다.

그렇다면 왜 이런 현상이 일어날까? 롯데리아가 시골을 좋아하는 걸까, 아니면 버거킹 등이 도심을 유난히 좋아하는 걸까? 아니면 또 다른 이유가 있을까? 좀 더 자세히 파 보자.

각 브랜드 별 매장 분포

우선 어느 버거 브랜드 매장이 어디에 있는지 수를 봐야한다.

역시 롯데리아 매장 수가 압도적으로 많다. 버거킹, KFC, 맥도날드 다 합쳐도 롯데리아 매장 수에 미치지 못한다. 단지 롯데리아 매장 수가 많기 때문에 좀 더 인구가 적은 지역까지 커버할 수 있다는 가설이 설득력이 생기는 차이다. 각 브랜드별 매장 분포를 좀 더 자세히 지역별로 살펴보자.

역시.. 롯데리아 색은 아주 예쁘게 고르게 분포한 반면에, 버거킹+맥도날드+KFC (줄여서 BMK) 모두 합쳐도 매장이 하나도 없는 곳이  많다. 롯데리아는 서울의 한 구보다 많은 곳이 지방에도 많은 반면에, BMK는 광역시 이상 도시보다 매장이 많은 지방이 거의 없다. 광역시 내에서도 다른 구들에 비해 “중구”에 쏠리는 현상이 BMK가 롯데리아보다 훨씬 뚜렷하다.

브랜드별로 각각 어느 브랜드가 어느 것과 비슷한지 좀 더 자세히 살펴보자.

피어슨 연관성 계수로 봤을 때, 의외의 결과! 맥도날드는 KFC나 버거킹보다 롯데리아에 더 가깝다. KFC와 버거킹은 서로 분포가 비슷하다. 맥도날드가 완전 도심보다 서울 안에서도 주로 주택가 상가에 많이 분포하는 걸 생각해보면 그럴 수도 있겠다 싶다. 버거지수 분자에서 맥도날드를 빼는게 더 원래 의도에 맞을 것 같다.

롯데리아 매장은 좀 더 인구가 적은 곳에서도 열 수 있는 것일까?

인구밀도와 버거지수의 관계

롯데리아가 인구밀도가 낮아도 먹고 살 수 있다면, 버거지수와 인구밀도가 어느 정도 상관관계가 있을 것이다. 비교해보면 진짜 그렇다! (다음 그림)

인구밀도가 높을 수록 (오른쪽으로 갈 수록), 버거지수가 높다. 켄달의 타오값도 0.618이나 된다. +_+ 원인은 분모 (롯데리아)에 있을 수도 있고 분자 (BMK)에 있을 수도 있겠다. 어느 것이 더 원인일까?

인구수에 정확하게 비례하는 모델로 매장을 흩어놓았을 때에 비해 얼마나 실제 매장이 많은지 (Y축)을 보면, 아주 강한 상관관계가 나오지는 않지만, 롯데리아는 어느 정도 이상 되면 인구밀도가 높을 수록 매장 수가 약간 줄어드는 경향이 있다. BMK의 경우에는 너무 흩어져 있어서 확실히 뭐라고 할 수는 없지만, 인구밀도가 낮은 지역에서 매장이 아예 없는 점이 상당히 많다. 전체 매장 수 자체가 적다보니 인구가 적은 지역을 커버하지 못하는 것으로도 볼 수 있겠다. 롯데리아가 인구밀도가 높을 수록 약간 줄어드는 것은 매장 간의 최소 거리를 고려하면 포화된 것으로도 해석할 수 있는데, 정확하게 보려면 좀 더 실제 데이터를 보고, 유동인구, 세대별 인구 같은 걸 봐야겠지만… 오바 같아서 그만 -ㅁ-;;

매장수 차이의 효과

마지막으로, 그냥 BMK가 각각 매장수가 적다보니 인구밀도가 적은 지역까지 진출하지 못한 효과가 버거지수를 만들 수 있는지 보기 위해, 각 시군구를 분할하면서 인구가 가장 많은 곳에 새로 매장을 만드는 식으로 시뮬레이션 했을 때 버거지수가 어떻게 나올지 확인해 보자. (자세한 것은 코드 확인.)

시뮬레이션에서 인구가 적어서 버거지수가 낮은 것으로 나온 곳은 실제 버거지수도 낮은 경향이 있었다. 하지만, 높게 예측된 곳 중에서는 인구가 큰 정보를 주지는 못했다. 특이한 점들을 몇 군데 보면, 다음 지역은 인구만 이용한 시뮬레이션에서 버거지수가 낮게 예측됐지만, 실제로는 매우 높은 점수가 나왔다. 과천시 빼고는 모두 인구에 비해 유동인구가 아주 많은 지역들이다.

시군구 롯데리아 버거킹 맥도날드 KFC
강원도 홍천군 1 2 0 0
인천광역시 중구 6 2 5 2
부산광역시 중구 4 1 1 3
경기도 과천시 3 1 1 1
강원도 정선군 1 0 0 1
대구광역시 중구 7 1 2 4

다음 지역들은 인구를 기반으로 한 시뮬레이션에서는 아주 낮지는 않은 버거지수가 나왔지만, 실제 버거지수는 매우 낮은 지역이다. BMK의 입지 조건에 맞는 곳이 없으려나..? 이유는 좀 더 살펴봐야겠다. -O-

시군구 롯데리아 버거킹 맥도날드 KFC
충청남도 보령시 3 0 0 0
경상북도 영천시 2 0 0 0
대구광역시 달성군 6 0 0 0

결론

롯데리아는 인구밀도가 많은 곳에서 오히려 인구 대비 매장수가 약간 적은 경향이 약간 있었고, BMK의 매장수가 전국을 덮을 수 있을 정도로 많지 않은 것이 또한 버거지수를 가르는 요인이었다. 다만, 이 글에서는 연관성에 집중했을 뿐, 이 현상의 원인에 대한 구체적인 인과관계를 알 수는 없으므로 해석은 좀 더 조심스러워야 한다. 또한, 시군구경계로 아주 단순화했고 유동인구나 인구 구조 등의 특징을 충분히 고려하지 못했기 때문에, 상권이나 구체적 인구 자료를 더 살펴봐야 할 것 같다. 이것은 소스와 데이터를 공개했으니 후행연구로 누군가… ㅎㅎ;