遺伝的背景で、遺伝子発現量を補正する

こんにちは。ここでは研究に役立つ情報をシェアしていきます。

こちらのGTExの論文では、800人以上の被験者から組織を採取し、それらの組織における遺伝子発現量を調べています。特に、本論文では遺伝子ごとに発現を解析し、その遺伝子発現が低すぎる人や高すぎる人("outlier")を抽出しています。その過程で、遺伝子発現量(TPM)を補正しているのですが、本ページではその補正方法を紹介します。GTExのデータはアクセス制限があるので、今回は適当に作ったテストデータで行います。

上記論文の補正方法は、ずばり、こちらです。こちらはSNPの主成分(PC)、つまり、遺伝的背景や性別等を共変量としたリニアーモデルを遺伝子発現量に回帰させ、その残差を補正済みの発現量として使用しています。

詳しくみてみましょう。

Yは実際の発現量(transform済)です。μは発現量の平均値です。GがPCで、βがGの係数(パラメータ)です。Linear modelの回帰ではこのβを推定します。βの絶対値が大きいほど、Yとμとの違い(=平均発現量と特定サンプルでの発現量の差)がGで説明されることとなります。PやQも同様の考え方です。

εは残差で、Yとμとの違い(=平均発現量と特定サンプルでの発現量の差)の内、GやPやQで説明できなかった値です。これらは、遺伝的背景やその他のバッチエフェクト等で説明できない値です。この残差には遺伝子近傍に存在する変異による影響を含むと考えられます。上記の論文では、遺伝近傍のレアバリアントがこの残差に寄与すると仮説立てています。

*****

早速こちらをテストしてみましょう。

まずはゲノムデータをダウンロードします。今回は1000人ゲノムプロジェクトのデータを用います。1000人ゲノムプロジェクトではおおよそ2500人のゲノムデータを使用できます。学術目的で無料でアクセスできます。今回はテストなので、chr22のみを使います。実際の解析では全てのChromosomeを用います。方法はこちらなどが参考になります。

chr22をダウンロード。

# download chr22
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/release/20181203_biallelic_SNV/ALL.chr22.shapeit2_integrated_v1a.GRCh38.20181129.phased.vcf.gz

次に、Quality control (QC)をします。PLINKと呼ばれるソフトウェアを用います。ハーディー・ワインベルグ平衡の状態にないSNP等を解析から除きます。

# QC
plink2 --threads 4 \
--vcf ALL.chr22.shapeit2_integrated_v1a.GRCh38.20181129.phased.vcf.gz \
--const-fid \
--allow-extra-chr 0 \
--geno 0.05 \
--hwe 1e-6 \
--maf 0.01 \
--make-pgen \
--out chr22_QCed

次にPruningをします。Pruningとは「LD pruning is a procedure of filtering genetic markers by linkage disequilibrium leaving for the analysis only tagging SNPs that are representatives of the genetic haplotype blocks.」(こちらより引用)です。

# pruning
plink2 --threads 4 \
--pfile chr22_QCed \
--maf 0.05 --indep-pairwise 100kb 0.5 \
--make-pgen \
--out chr22_pruned

そして、PCAを行います。

# PCA
plink2 --threads 4 \
--pfile chr22_pruned \
--pca \
--out chr22_pca

ここまで来れば、「chr22_pca.eigenvec」というファイルがあるはずです。「eigen」は「アイゲン」と発音します(アメリカ発音)。こちらにTop10個のPCが記載されています。今回はこれらのPCを用い、遺伝子発現量を補正してみましょう。

次に、テストデータの遺伝子発現量を作ります。ここからPython3.7です。「id_to_pop_crossref_1kGP.txt」は筆者が作成したファイルです。1000人ゲノムプロジェクトのサイトからダウロードしたファイルを変形したものです。サンプル名とそのサンプルのSuper-populationが記載されています。このファイルはこちらにアップロードしています。

# Python3.7

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt
from scipy import stats as st
from sklearn.linear_model import LinearRegression

# make test data
f='id_to_pop_crossref_1kGP.txt'
out=[]
with open(f) as infile:
   for line in infile:
       ls=line.split()
       if ls[2] == 'EUR':
           if np.random.rand() < 0.1:
               cpm=np.random.normal(150, 7.5, 1)[0]
           else:
               cpm=np.random.normal(100, 5, 1)[0]
       elif ls[2] == 'AFR':
           cpm=np.random.normal(200, 10, 1)[0]
       else:
           cpm=np.random.normal(100, 5, 1)[0]
       out.append('%s\t%s\t%f\n' % (ls[0], ls[2], cpm))
with open('CPM.txt', 'w') as outfile:
   outfile.write(''.join(out))

今回作ったテストデータを可視化してみましょう。

# load test data
df_cpm=pd.read_csv('CPM.txt', index_col=0, header=None, sep='\t', names=('ancestry', 'cpm'))


# plot test data
fig=plt.figure(figsize=(3, 3))  # (x, y)
ax=fig.add_subplot(111)
sns.swarmplot(data=df_cpm, x='ancestry', y='cpm', size=1)
plt.suptitle('CPM')
plt.savefig('plot_orig_cpm.pdf')
plt.close()

1000人ゲノムプロジェクトには5つのSuper-populationが含まれています。
今回は、EAS、AMR、SASと比較しAFRが2倍のCPMを持ち、EURの10%が1.5倍のCPMを持つデータを作成しました。

今回のデータには2つの特徴があります。
   1) AFRの発現量が高い
   2) EURの一部の発現量が高い
上記の2つは何が原因でしょうか？遺伝的要因でしょうか？それともシーケンシングプラットフォームの違いでしょうか？
→ もし、AFRの発現量が高い原因が、AFRに共通した(主成分分析で捉えることのできる)遺伝的要因で説明がつくのであれば、今回のテストのようにSNPのPCを用いることで補正できるはずです。

それでは実際に、テストデータ(CPM)を遺伝的背景で補正してみましょう。
今回の補正では、以下のようにモデルを作ります。
　Y: 遺伝子発現量
　µ: 遺伝子発現量の平均
　G: PC(主成分)
　β: 共変量Gの係数 (回帰でこれを推定します)
　ε: 残差 (最終的にこれを知りたい)

実際のコードでは以下のように回帰を行っています。σは切片です。

回帰をした後、残差を計算します。

# load test data
df_cpm=pd.read_csv('CPM.txt', index_col=0, header=None, sep='\t', names=('ancestry', 'cpm'))

# load PC
df_pc=pd.read_csv('chr22_pca.eigenvec', index_col=0, sep='\t')

# join CPM and PC
df=df_cpm.join(df_pc, how='inner')

# convert data to zscore
X=st.zscore(df[df_pc.columns])
y= df['cpm'] - df['cpm'].mean()
y=st.zscore(y).reshape(-1, 1)

# linear regression
reg=LinearRegression().fit(X, y)

# calc. epsilon
eps= y.flatten() - (np.dot(X, reg.coef_.flatten()) + reg.intercept_)
eps=st.zscore(eps)
df_normalized=pd.DataFrame(eps, index=df.index, columns=['norm_cpm'])
df_normalized=df_normalized.join(df_cpm, how='inner')

さて、これで補正ができました。補正にはscikit-learnのLinearRegressionを用いました。補正前に全ての値をzscoreに変換しています。zscoreの変換は目的に応じてすると良いと思います。
補正した値をプロットしてみましょう。

# plot normalized value
fig=plt.figure(figsize=(3, 3))  # (x, y)
ax=fig.add_subplot(111)
sns.swarmplot(data=df_normalized, x='ancestry', y='norm_cpm', size=1)
plt.suptitle('Normalized CPM, zscore')
plt.savefig('plot_norm_cpm.pdf')
plt.close()

AFRの値が他の集団と近くなっています。EURの一部は補正されず、いまだ高い値のままです。

今回の補正では、以下の２点に着目します。
    1) AFRの発現量が下がった(他の集団と同じくらいになった)
    2) EURの一部の発現量が高い傾向はそのままである。
AFRの発現量が補正されたことから、AFRの発現量の高さはAFRに共通する遺伝的要因によることが考えられます。一方、EURの10%の発現量の高さは補正されませんでした。これはSNPのPCでは捉えられない遺伝的要因、もしくは遺伝的要因以外の何かにより遺伝子発現の違いが生まれていることを示唆します。

今回は、Top10個のPCを共変量とし、Linear regressionを行いました。それらのPCの係数とを見てみましょう。

df_coef=pd.DataFrame(reg.coef_.flatten(), index=df_pc.columns, columns=['coef'])
print(df_coef)

          coef
PC1  -0.945504
PC2   0.096534
PC3   0.025329
PC4  -0.005085
PC5   0.012873
PC6   0.001535
PC7   0.001018
PC8   0.018514
PC9  -0.023759
PC10  0.029497

係数を見ると、最も絶対値が大きい係数はPC1です。このことから、PC1が発現量の補正に大きく寄与していることがわかります。PC1がヒト集団の遺伝的背景をどのように捉えているのかをみてみましょう。

fig=plt.figure(figsize=(3, 3))  # (x, y)
ax=fig.add_subplot(111)
sns.scatterplot(data=df, x='PC1', y='PC2', hue='ancestry', s=3, alpha=0.3)
plt.suptitle('PC1 vs PC2')
plt.savefig('plot_PC1_PC2.pdf')
plt.close()

PC1は主にAFRとそれ以外の集団の遺伝的背景の違いを捉えています。テストデータはAFRのみ平均値を大きく変えたのでPC1が補正に大きく寄与したのでしょう。

*****

今回は、多様なヒト集団由来のビッグデータの補正を紹介しました。今後、ビッグデータが増えていく中、ヒト集団のgenetic ancestryを考慮し解析することは、False positiveを減らすとともによりクリアな結果を得るために重要だと思います。
本手法はシーケンスプラットフォームの違い等の補正にも用いられています。上記のGTEx論文では「PEER factor」という値を用い、"unknown batch effect"等を補正しています。こちらも近いうちに紹介できればと考えています。

*****

今回のコードはこちらからもダウンロード可能です。
実行環境：
Python 3.7.4
Ubuntu 18.04.2 LTS

*****

本記事の著者について

本記事はウイルス学若手ネットワークのメンバーが書いています。ウイルス学内外の若手の研究のお役に立てればと思っています。
ご意見やご感想、間違い等がありましたら以下のツイッターのDMにご連絡ください。もちろんフォローも忘れずにお願いします！！

Tweets by VirologyWakate

*****

以下、免責事項です。ご一読ください。
本ページの内容の正確性は一切保証しません。本情報は著者の所属機関を代表する見解ではありません。本情報を掲載された情報・資料を利用、使用、ダウンロードするなどの行為に関連して生じたあらゆる損害等に関しても、その理由に関わらず一切責任を負いません。本ページの内容は研究手法の解説および研究を行う上での知識共有が目的であり、それ以外の意図は一切ありません。本ページのコンテンツは、予告なく変更や公開の取り消しする場合があります。