遺伝子発現解析

問題: パン酵母 (出芽酵母) の遺伝子発現の解析

目標は、一般にパン酵母または醸造用酵母として知られる出芽酵母の遺伝子発現について理解することです。出芽酵母は、パンを焼いたり、ブドウからワインを発酵させたりするのに使用する菌類です。

出芽酵母をブドウ糖が豊富な培地に加えると、ブドウ糖はエタノールに変換されます。酵母は、まず "発酵" と呼ばれる代謝プロセスによってブドウ糖をエタノールに変換します。ただし、ブドウ糖の供給が尽きると、酵母は、ブドウ糖の嫌気性発酵からエタノールの好気呼吸に移行します。このプロセスは、ジオキシック シフトと呼ばれます。これは、遺伝子発現の大きな変化が伴うため、非常に興味深いプロセスです。

この例では、DNA マイクロアレイ データを使用し、ジオキシック シフト中の出芽酵母におけるほぼすべての遺伝子の一時的な遺伝子発現を調べます。

この例を実行するには、Bioinformatics Toolbox™ が必要です。

if ~nnDependency.bioInfoAvailable
    errordlg('This example requires Bioinformatics Toolbox.');
    return;
end

データ

この例では、次のデータを使用します。DeRisi, JL, Iyer, VR, Brown, PO. "Exploring the metabolic and genetic control of gene expression on a genomic scale." Science. 1997 Oct 24;278(5338):680-6. PMID: 9381177

データセット全体は、Gene Expression Omnibus の Web サイト https://www.yeastgenome.org からダウンロードできます。

最初に、データを MATLAB® に読み込みます。

load yeastdata.mat

遺伝子発現レベルはジオキシック シフト中に 7 つの時間点で測定されました。変数 times には、実験で発現レベルが測定された時間が含まれます。変数 genes には、発現レベルが測定された遺伝子の名前が含まれます。変数 yeastvalues には、"VALUE" データか、LOG_RAT2N_MEAN (実験における 7 つのタイム ステップから取得した CH2DN_MEAN と CH1DN_MEAN の比率の log2) が含まれます。

データのサイズを把握するには、numel(genes) を使用してデータ セットに含まれる遺伝子の数を表示します。

numel(genes)

ans = 6400

genes は遺伝子の名前の cell 配列です。MATLAB の cell 配列のインデックス付けを使用してエントリにアクセスできます。

genes{15}

ans = 
'YAL054C'

これは、変数 yeastvalues の 15 行目に ORF YAL054C の発現レベルが含まれていることを示しています。

遺伝子のフィルター処理

データセットはかなり大きく、情報の多くは、実験中に興味深い変化を示さない遺伝子に対応しています。興味深い遺伝子を簡単に見つけるためには、まず最初に興味深い変化を示さない発現プロファイルの遺伝子を削除して、データセットのサイズを小さくします。6400 個の発現プロファイルがあります。複数の手法を使用し、これを最も重要な遺伝子を含むいくつかのサブセットに減らすことができます。

遺伝子リストに目を通すと、'EMPTY' としてマークされたいくつかのスポットを確認できます。これらは配列の空のスポットであり、データが関連付けられている可能性もありますが、この例ではこれらのポイントをノイズと見なすことができます。これらのポイントを関数 strcmp を使用して検出し、インデックス付けのコマンドでデータ セットから削除できます。

emptySpots = strcmp('EMPTY',genes);
yeastvalues(emptySpots,:) = [];
genes(emptySpots) = [];
numel(genes)

ans = 6314

yeastvalues のデータで、発現レベルが NaN としてマークされている場所もいくつか確認できます。これは、特定のタイム ステップでこのスポットのデータが収集されなかったことを示します。これらの欠損値を処理する 1 つの方法は、特定の遺伝子のデータの平均値または中央値を使用して、欠損値を経時的に補完することです。この例では、1 つ以上の発現レベルが測定されなかった遺伝子のデータを破棄するだけという、あまり厳密ではないアプローチを使用します。

関数 isnan を使用して欠損データをもつ遺伝子を特定し、インデックス付けのコマンドを使用して欠損データをもつ遺伝子を削除します。

nanIndices = any(isnan(yeastvalues),2);
yeastvalues(nanIndices,:) = [];
genes(nanIndices) = [];
numel(genes)

ans = 6276

残りのプロファイルすべての発現プロファイルをプロットすると、ほとんどのプロファイルはフラットで、他のプロファイルと大きく異ならないことを確認できます。このフラット データは、これらのプロファイルと関連付けられている遺伝子が、ジオキシック シフトの影響をあまり受けないことを示す点で確かに役立ちます。しかし、この例では、ジオキシック シフトに伴う発現の変化が大きい遺伝子に注目します。Bioinformatics Toolbox™ のフィルター処理関数を使用し、代謝の変化の影響を受ける遺伝子に関する有用な情報が得られない、さまざまなタイプのプロファイルをもつ遺伝子を削除できます。

関数 genevarfilter を使用し、時間の経過に沿って分散が小さい遺伝子を除外できます。この関数は、可変遺伝子と同じサイズの logical 配列を返します。1 は分散が 10 番目の百分位数よりも大きい yeastvalues の行に、0 はこのしきい値を下回る yeastvalues の行に対応します。

mask = genevarfilter(yeastvalues);
% Use the mask as an index into the values to remove the filtered genes.
yeastvalues = yeastvalues(mask,:);
genes = genes(mask);
numel(genes)

ans = 5648

関数 genelowvalfilter は、非常に低い絶対発現値をもつ遺伝子を削除します。遺伝子フィルター関数は、フィルター処理されたデータおよび名前の自動計算もできることに注意してください。

[mask, yeastvalues, genes] = ...
    genelowvalfilter(yeastvalues,genes,'absval',log2(3));
numel(genes)

ans = 822

geneentropyfilter を使用してプロファイルのエントロピーが低い遺伝子を削除します。

[mask, yeastvalues, genes] = ...
    geneentropyfilter(yeastvalues,genes,'prctile',15);
numel(genes)

ans = 614

主成分分析

扱いやすい遺伝子リストができたので、プロファイル間の関係を調べることができます。

データの標準偏差と平均を正規化することで、ネットワークは、各入力をその値の範囲全体で等しく重要なものとして扱うことができます。

主成分分析 (PCA) は、マイクロアレイ解析などからの大規模なデータセットの次元を削減するために使用できる便利な手法です。この手法は、データセットの主成分を分離し、データセットの変動への影響が最も少ない成分を削除します。

2 つの設定変数を使用して mapstd と processpca を新しいデータに適用し、整合性を保つことができます。

[x,std_settings] = mapstd(yeastvalues');  % Normalize data
[x,pca_settings] = processpca(x,0.15);    % PCA

まず、入力ベクトルがゼロ平均と単位分散をもつように mapstd を使用して正規化されます。processpca は、PCA アルゴリズムを実装する関数です。processpca に渡される 2 番目の引数は 0.15 です。これは、processpca によって、データセットの全変動への影響が 15% 未満の主成分が排除されることを意味します。ここで、変数 pc は yeastvalues データの主成分を含みます。

主成分は関数 scatter を使用して可視化できます。

figure
scatter(x(1,:),x(2,:));
xlabel('First Principal Component');
ylabel('Second Principal Component');
title('Principal Component Scatter Plot');

クラスター分析: 自己組織化マップ

自己組織化マップ (SOM) クラスタリング アルゴリズムを使用して主成分をクラスター化できるようになりました。

関数 selforgmap は、自己組織化マップ ネットワークを作成します。その後、関数 train を使用してこのネットワークの学習を行うことができます。

ネットワークはまだ入力データに一致するように構成されていないため、入力のサイズは 0 です。ネットワークの学習時にはこのようになります。

net = selforgmap([5 3]);
view(net)

これでネットワークの学習の準備が整いました。

Neural Network Training Tool を使用すると、学習するネットワークと、学習に使用されているアルゴリズムが表示されます。さらに、学習中には学習の状態が表示され、学習を停止した条件が緑で強調表示されます。

下部にあるボタンを使用すると、便利なプロットを開くことができます。これらのプロットは、学習中および学習後に開くことができます。アルゴリズム名およびプロット ボタンの隣のリンクを使用すると、これらに関するドキュメンテーションを開くことができます。

net = train(net,x);

plotsompos を使用し、データの最初の 2 つの次元に関する散布図にネットワークを重ねて表示します。

figure
plotsompos(net,x);

データセットの各ポイントに最も近いノードを見つけることにより、SOM を使用してクラスターを割り当てることができます。

y = net(x);
cluster_indices = vec2ind(y);

plotsomhits を使用し、マップの各ニューロンに割り当てられているベクトルの数を確認します。

figure
plotsomhits(net,x);

クラスター分析には、階層クラスタリングや k-means など、Statistics and Machine Learning Toolbox™ で利用可能な他のクラスタリング アルゴリズムも使用できます。

用語

ORF – オープン リーディング フレーム (ORF) は、終止配列によって中断されない、塩基配列を含む遺伝子配列の一部であり、タンパク質を符号化する可能性があります。