🙆‍♀️

eQTLの実装 (Omics)

2024/10/12に公開

はじめに

これまでプロテオーム・トランスクリプトーム・エピゲノムのデータ解析を行うことが多かったが、ゲノム解析は全く経験が無かった。ゲノムとトランスクリプトームのマルチオミックスな解析手法であるeQTL解析の実装を通してゲノム解析の一端に触れてみることにした。

eQTLの概要

eQTLは遺伝子の発現と有意に関連のあるSNPを同定する方法である。遺伝子発現量の変化は、細胞内・生体内の機能に影響を与える中間
形質（endophenotypge）であるため、eQTL解析を行うことで、遺伝子多型
と疾患発症のつながりを検討できると考えられています。

具体的な用途としては、複数の転写産物に影響を及ぼすDNA領域を決定したり、表現型のサブクラスを見つけるために利用するらしい。一口にeQTLといっても様々な手法が存在するが、線形回帰モデルやANOVAモデルや線形混合モデルなどを利用して解析を行うのが一般的である。2012年のMatrix eQTLの原著論文では、データの次元数が増えたとしてもPlinkなどと比較してMatrix eQTLの速度はそれほど低下しない（2024年現在の状況は調べていません）。

eQTLの実装

MatrixEQTLというRパッケージがeQTL解析ではよく使用される。そこで、このRパッケージを利用してeQTL解析を実演する。遺伝子とSNP（単一塩基多型）間の関連性を線形回帰モデルを使って解析する。

パッケージのインストールとロード

pkgs <- c("MatrixEQTL", "tidyverse")
library(BiocManager)
for (pkg in pkgs) if (!requireNamespace(pkg)) {
  BiocManager::install(pkg, update = FALSE)
}
for (pkg in pkgs) {
  library(pkg, character.only = TRUE)
}

基本的なファイル設定

base.dir変数でMatrixEQTLパッケージのインストール場所を取得し、次にSNPデータや遺伝子発現データ、共変量データのファイル名を設定する。useModelには解析モデルとして線形回帰（modelLINEAR）を選択した。共変量データを使わない場合は、covariates_file_nameを空文字列にすることで無視できる。

base.dir <- find.package("MatrixEQTL")
useModel <- modelLINEAR
SNP_file_name <- paste0(base.dir, "/data/SNP.txt")
expression_file_name <- paste0(base.dir, "/data/GE.txt")
covariates_file_name <- paste0(base.dir, "/data/Covariates.txt")
output_file_name <- tempfile()

p値のしきい値設定と誤差共分散行列

p値の閾値を設定する。解析結果の出力ファイルには、p値がこの閾値以下であったSNP-遺伝子ペアが保存される。errorCovarianceには誤差の共分散行列を指定できるが、通常は利用することが少なくnumeric()としておけば問題ない。

pvOutputThreshold <- 1e-2
errorCovariance <- numeric()

ファイルのロード

snps <- SlicedData$new()
snps$fileDelimiter = "\t"
snps$fileOmitCharacters = "NA"
snps$fileSkipRows = 1
snps$fileSkipColumns = 1
snps$fileSliceSize = 2000
snps$LoadFile( SNP_file_name )

gene <- SlicedData$new()
gene$fileDelimiter = "\t"
gene$fileOmitCharacters = "NA"
gene$fileSkipRows = 1
gene$fileSkipColumns = 1
gene$fileSliceSize = 2000
gene$LoadFile( expression_file_name )

cvrt = SlicedData$new()
cvrt$fileDelimiter = "\t"
cvrt$fileOmitCharacters = "NA"
cvrt$fileSkipRows = 1
cvrt$fileSkipColumns = 1
cvrt$fileSliceSize = 2000
cvrt$LoadFile( covariates_file_name )

15サンプル、15SNPの遺伝子型情報（0, 1, 2にエンコードされている）を解析に用いる。

> snps
SlicedData object. For more information type: ?SlicedData
Number of columns: 16 
Number of rows: 15 
Data is stored in 1 slices
Top left corner of the first slice (up to 10x10):
       Sam_01 Sam_02 Sam_03 Sam_04 Sam_05 Sam_06 Sam_07 Sam_08 Sam_09
Snp_01      2      0      2      0      2      1      2      1      1
Snp_02      0      1      1      2      2      1      0      0      0
Snp_03      1      0      1      0      1      1      1      1      0
Snp_04      0      1      2      2      2      1      1      0      0
Snp_05      1      1      2      1      1      2      1      1      0
Snp_06      2      2      2      1      1      0      1      0      2
Snp_07      1      1      2      2      0      1      1      1      1
Snp_08      1      0      1      0      1      0      0      1      1
Snp_09      2      1      2      2      0      1      1      0      2
Snp_10      1      1      0      0      0      2      2      1      1
       Sam_10
Snp_01      1
Snp_02      1
Snp_03      1
Snp_04      0
Snp_05      1
Snp_06      1
Snp_07      0
Snp_08      1
Snp_09      1
Snp_10      2

15サンプル、10遺伝子の遺伝子発現データを解析に用いる。

> gene
SlicedData object. For more information type: ?SlicedData
Number of columns: 16 
Number of rows: 10 
Data is stored in 1 slices
Top left corner of the first slice (up to 10x10):
        Sam_01 Sam_02 Sam_03 Sam_04 Sam_05 Sam_06 Sam_07 Sam_08 Sam_09
Gene_01   4.91   4.63   5.18   5.07   5.74   5.09   5.31   5.29   4.73
Gene_02  13.78  13.14  13.18  13.04  12.85  13.07  13.09  12.83  14.94
Gene_03  12.06  12.29  13.07  13.65  13.86  12.84  12.29  13.03  13.13
Gene_04  11.63  11.88  12.74  12.66  13.16  11.99  11.97  12.81  11.74
Gene_05  14.72  14.66  14.63  15.91  15.46  14.74  15.17  15.01  14.05
Gene_06  12.29  12.23  12.47  13.21  12.63  12.18  12.44  12.45  13.22
Gene_07  12.56  12.71  12.49  13.41  13.60  12.35  12.27  13.42  14.73
Gene_08  12.27  12.58  12.61  13.02  12.86  12.32  12.30  13.19  15.21
Gene_09   9.82   9.29   8.95   8.18   8.11   9.43   9.17   9.25  10.95
Gene_10  14.24  14.52  14.62  13.65  13.46  14.04  13.97  13.73  12.29
        Sam_10
Gene_01   5.72
Gene_02  13.86
Gene_03  14.93
Gene_04  13.29
Gene_05  14.90
Gene_06  12.70
Gene_07  15.47
Gene_08  14.60
Gene_09   9.82
Gene_10  12.01

今回は15サンプルのgenderとageのデータを共変量に用いる。

> cvrt 
SlicedData object. For more information type: ?SlicedData
Number of columns: 16 
Number of rows: 2 
Data is stored in 1 slices
Top left corner of the first slice (up to 10x10):
       Sam_01 Sam_02 Sam_03 Sam_04 Sam_05 Sam_06 Sam_07 Sam_08 Sam_09
gender      0      0      0      0      0      0      0      0      1
age        36     40     46     65     69     43     40     54     44
       Sam_10
gender      1
age        70

メインのeQTL解析関数

Matrix_eQTL_engine関数を用いてeQTL解析を実行する。SNPデータ・遺伝子発現データ・共変量データを渡し、指定されたモデルやp値の閾値に基づいて解析が行われる。ここでは、遺伝子発現量を従属変数として線形回帰分析を実施する。

me = Matrix_eQTL_engine(
    snps = snps,
    gene = gene,
    cvrt = cvrt,
    output_file_name = output_file_name,
    pvOutputThreshold = pvOutputThreshold,
    useModel = useModel,
    errorCovariance = errorCovariance,
    verbose = TRUE,
    pvalue.hist = TRUE,
    min.pv.by.genesnp = FALSE,
    noFDRsaveMemory = FALSE)

# Processing covariates
# Task finished in 0.02 seconds
# Processing gene expression data (imputation, residualization)
# Task finished in 0 seconds
# Creating output file(s)
# Task finished in 0.01 seconds
# Performing eQTL analysis
# 100.00% done, 5 eQTLs
# Task finished in 0 seconds

解析結果

解析結果はmeオブジェクト（リスト）の中に以下の通り出力された

> show(me$all$eqtls)
    snps    gene statistic       pvalue          FDR       beta
1 Snp_05 Gene_03 38.812160 5.515519e-14 8.273279e-12  0.4101317
2 Snp_13 Gene_09 -3.914403 2.055817e-03 1.541863e-01 -0.2978847
3 Snp_11 Gene_06 -3.221962 7.327756e-03 2.853368e-01 -0.2332470
4 Snp_04 Gene_10  3.201666 7.608981e-03 2.853368e-01  0.2321123
5 Snp_14 Gene_01  3.070005 9.716705e-03 2.915011e-01  0.2147077

全pvalueのヒストグラムを作成する

plot(me)

cis/trans eQTL解析

その他にも、Matrix eQTLは遺伝子近傍のSNP（cis-eQTL）と、離れた領域のSNP（trans-eQTL）を区別して解析できる。これらについて解析を行いたい場合は、前述のMatrix_eQTL_engine()関数ではなくてMatrix_eQTL_main()関数を用いるとともにいくつかのパラメータを設定する必要がある。具体的には、以下の４つのパラメータを定義した際にgene-SNPのペア間の距離がcisDist以下であった場合、そのペアはcisであると判定される。

pvOutputThreshold.cis > 0
cisDist
snpspos
genepos

# source("Matrix_eQTL_R/Matrix_eQTL_engine.r");
library(MatrixEQTL)

## Location of the package with the data files.
base.dir = find.package('MatrixEQTL');

## Settings
# Linear model to use, modelANOVA, modelLINEAR, or modelLINEAR_CROSS
useModel = modelLINEAR; # modelANOVA, modelLINEAR, or modelLINEAR_CROSS

# Genotype file name
SNP_file_name = paste0(base.dir, "/data/SNP.txt");
snps_location_file_name = paste0(base.dir, "/data/snpsloc.txt");

# Gene expression file name
expression_file_name = paste0(base.dir, "/data/GE.txt");
gene_location_file_name = paste0(base.dir, "/data/geneloc.txt");

# Covariates file name
# Set to character() for no covariates
covariates_file_name = paste0(base.dir, "/data/Covariates.txt");

# Output file name
output_file_name_cis = tempfile();
output_file_name_tra = tempfile();

# Only associations significant at this level will be saved
pvOutputThreshold_cis = 2e-2;
pvOutputThreshold_tra = 1e-2;

# Error covariance matrix
# Set to numeric() for identity.
errorCovariance = numeric();
# errorCovariance = read.table("Sample_Data/errorCovariance.txt");

# Distance for local gene-SNP pairs
cisDist = 1e6;

## Load genotype data
snps = SlicedData$new();
snps$fileDelimiter = "\t";      # the TAB character
snps$fileOmitCharacters = "NA"; # denote missing values;
snps$fileSkipRows = 1;          # one row of column labels
snps$fileSkipColumns = 1;       # one column of row labels
snps$fileSliceSize = 2000;      # read file in slices of 2,000 rows
snps$LoadFile(SNP_file_name);

## Load gene expression data
gene = SlicedData$new();
gene$fileDelimiter = "\t";      # the TAB character
gene$fileOmitCharacters = "NA"; # denote missing values;
gene$fileSkipRows = 1;          # one row of column labels
gene$fileSkipColumns = 1;       # one column of row labels
gene$fileSliceSize = 2000;      # read file in slices of 2,000 rows
gene$LoadFile(expression_file_name);

## Load covariates
cvrt = SlicedData$new();
cvrt$fileDelimiter = "\t";      # the TAB character
cvrt$fileOmitCharacters = "NA"; # denote missing values;
cvrt$fileSkipRows = 1;          # one row of column labels
cvrt$fileSkipColumns = 1;       # one column of row labels
if(length(covariates_file_name)>0) {
  cvrt$LoadFile(covariates_file_name);
}

## Run the analysis
snpspos = read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE);
genepos = read.table(gene_location_file_name, header = TRUE, stringsAsFactors = FALSE);

me = Matrix_eQTL_main(
    snps = snps, 
    gene = gene, 
    cvrt = cvrt,
    output_file_name      = output_file_name_tra,
    pvOutputThreshold     = pvOutputThreshold_tra,
    useModel = useModel, 
    errorCovariance = errorCovariance, 
    verbose = TRUE, 
    output_file_name.cis  = output_file_name_cis,
    pvOutputThreshold.cis = pvOutputThreshold_cis,
    snpspos = snpspos, 
    genepos = genepos,
    cisDist = cisDist,
    pvalue.hist = TRUE,
    min.pv.by.genesnp = FALSE,
    noFDRsaveMemory = FALSE);

unlink(output_file_name_tra);
unlink(output_file_name_cis);

## Results:

cat('Analysis done in: ', me$time.in.sec, ' seconds', '\n');
cat('Detected local eQTLs:', '\n');
show(me$cis$eqtls)
cat('Detected distant eQTLs:', '\n');
show(me$trans$eqtls)

## Make the histogram of local and distant p-values

plot(me)

> show(me$cis$eqtls)
    snps    gene statistic       pvalue          FDR      beta
1 Snp_05 Gene_03 38.812160 5.515519e-14 5.515519e-12 0.4101317
2 Snp_04 Gene_10  3.201666 7.608981e-03 3.804491e-01 0.2321123
> show(me$trans$eqtls)
    snps    gene statistic      pvalue       FDR       beta
1 Snp_13 Gene_09 -3.914403 0.002055817 0.1027908 -0.2978847
2 Snp_11 Gene_06 -3.221962 0.007327756 0.1619451 -0.2332470
3 Snp_14 Gene_01  3.070005 0.009716705 0.1619451  0.2147077

解析結果の可視化

せっかくなので可視化してみた。

expression_file_name <- paste0(base.dir, "/data/GE.txt")
df_ge <- read.table(expression_file_name, header = TRUE, row.names = 1)
sampleids <- colnames(df_ge)
df_ge_long <- df_ge %>% mutate(geneid = rownames(.)) %>% 
  pivot_longer(cols = all_of(sampleids), names_to = "sample", values_to = "expression")

snp_file_name <- paste0(base.dir, "/data/SNP.txt")
df_snp <- read.table(snp_file_name, header = TRUE, row.names = 1)
df_snp_long <- df_snp %>% mutate(snp = rownames(.)) %>% 
  pivot_longer(cols = all_of(sampleids), names_to = "sample", values_to = "genotype")

df_gene_03 <- merge(df_ge_long, df_snp_long, by = "sample", all = F) %>% 
  filter(geneid == "Gene_03" & snp == "Snp_05")
df_gene_03$genotype <- as.factor(df_gene_03$genotype)
p <- ggplot(df_gene_03, aes(x = genotype, y = expression, group = genotype)) +
  geom_boxplot(aes(fill = genotype)) +
  geom_dotplot(binaxis = "y") +
  scale_color_brewer()
print(p)

参考

Matrix eQTLの各種ページへのリンク
Matrix eQTLの原著論文
Github - MatrixEQTL
遺伝統計学・夏の学校講義実習資料　GenomeDataAnalysis２
Introduction to eQTL analysis ←公式ドキュメントを個人で解説した資料
A pipeline for RNA-seq based eQTL analysis with automated quality control procedures ← 本記事では触れないが、Matrix eQTLを実行する際の前処理の参考になりそうな論文