学習内容

1. 尤度比検証（LRT）を仮説検定に応用する方法
2. LRTで生成された結果をWald検証の結果と比較する
3. LRTの有意な遺伝子リストから共通の発現パターンを特定する

尤度比検証

2つ以上のレベルにおける発現変化を評価する際、DESeq2はWald検証に代わる方法として尤度比検証（LRT）を提供します。重要と判定される遺伝子は、因子レベル間でいずれの方向にも発現が変化する遺伝子です。

通常、この検証は個別の対比比較よりも多くの遺伝子を生成します。LRTは因子の任意のレベル間の差の有意性検証ですが、Wald検証を使用した遺伝子セットの和集合と完全に一致することは期待できません（ただし、高度に重なることは期待できます）。

結果抽出

dds_lrtオブジェクトから結果を抽出するには、Wald検証と同じresults()関数を使用できます。対比比較を行っていないため、引数は不要です。

# LRTの結果を抽出
lrt_results <- extract_findings(lrt_model)

結果テーブルを確認してみましょう：

# LRTの結果を表示
lrt_results  

出力はWald検証の結果に似ており、以前観察されたものと同じ列を持っています。

• なぜLRT検証では倍数変化を報告するのか？

尤度比検証を使用した分析では、p値は完全なモデル式と簡略化されたモデル式間の偏差差によってのみ決定されます。個々のlog2倍数変化は結果テーブルに印刷されますが、他の結果テーブル出力と一貫性を持たせるためであり、実際の検証とは関係ありません。

LRT検証に関連する列：

• baseMean：すべてのサンプルの正規化カウントの平均値
• statistic：簡略化モデルと完全モデル間の偏差差
• pvalue：統計値をカイ二乗分布と比較してp値を生成
• padj：BH調整後のp値

追加列：

• log2FoldChange：log2倍数変化
• lfcSE：標準誤差

重要な遺伝子の特定

LRTから重要な遺伝子をフィルタリングする際には、padj列にしきい値のみを設定します。padj < 0.05の場合、どれくらいの遺伝子が有意ですか？

# LRT結果用のtibbleを作成
lrt_table <- lrt_results %>%
  convert_to_dataframe() %>%
  add_rownames(column="gene") %>% 
  create_tibble()

# padj < 0.05の遺伝子をサブセット化
significant_lrt_genes <- lrt_table %>% 
  filter(padj < significance_threshold)

# 有意な遺伝子の数を取得
count_genes(significant_lrt_genes)

# Wald検証で得た数と比較
count_genes(overexpressed_genes)
count_genes(knockdown_genes)

LRTから観察された重要な遺伝子の数はかなり多いです。このリストには、3つの因子レベル（対照、ノックアウト、過発現）のいずれかの方向で変化する遺伝子が含まれます。重要な遺伝子の数を減らすには、FDRしきい値（significance_threshold）を厳しくすることができます。

• 共通の発現パターンを持つ遺伝子クラスターの特定

約7Kの重要な遺伝子のリストがあり、これらの遺伝子は3つの異なるサンプルグループで何らかの方法で変化することがわかっています。次に何をすべきでしょうか？

次のステップは、サンプルグループ（レベル）間で共通の発現変化パターンを持つ遺伝子グループを特定することです。これには、"DEGreport"パッケージから"findExpressionClusters"という名前のクラスタリングツールを使用します。"findExpressionClusters"ツールは、遺伝子間のペア相関に基づく階層的クラスタリング方法を使用し、次に階層木をカットして類似した発現パターンを持つ遺伝子グループを生成します。このツールは、クラスター間の変動性>クラスター内の変動性になるように、木をカットすることでクラスター多様性を最適化します。

クラスタリングを開始する前に、まず"rlog"変換正規化カウントをサブセット化して、差発現遺伝子（padj < 0.05）のみを保持します。この例では、7Kの遺伝子でクラスタリングを実行するには時間がかかるため、デモ目的で、調整p値でソートされた上位1000の遺伝子のみをサブセット化します。

# よ速いクラスタ検索のための結果サブセット化（デモ目的）
cluster_candidates <- significant_lrt_genes %>%
  sort_by_padj() %>%
  select_top_genes(n=1000)


# それらの有意な遺伝子のrlog値を取得
clustered_data <- rlog_transformed_matrix[cluster_candidates$gene, ]

有意な遺伝子の"rlog"変換カウントといくつかの追加パラメータが"findExpressionClusters"に入力されます：

• sample_info：サンプルに対応するメタデータdataframe
• condition：メタデータ内の文字列列名、変化の変数として使用
• group：サンプルを分離するためのメタデータ内の文字列列名

# 'DEGreport'パッケージの'findExpressionClusters'関数を使用してサンプルグループ間の遺伝子クラスターを表示
expression_clusters <- findExpressionClusters(clustered_data, metadata = sample_metadata, condition = "sampletype", group=NULL)

クラスタリングの実行が完了すると、コンソールにコマンドプロンプトが返され、プロットウィンドウにグラフが表示されるはずです。これらの遺伝子は4つの異なるグループに分類されます。各遺伝子グループについて、異なるサンプルグループ間の発現変化を示すボックスプロットがあります。発現変化の傾向を示す折れ線グラフが重ねて表示されています。

サンプルで発現の減少と過発現の増加を示す遺伝子に興味があると仮定しましょう。この図によると、この発現パターンを共有する遺伝子は275個あります。これらの遺伝子を特定するために、出力を探索してみましょう。クラスタリング出力のデータ構造はどのタイプですか？

# 'expression_clusters'の出力のデータ構造はどのタイプですか？
data_type(expression_clusters)

名前（クラスター）を使用してリストにどのオブジェクトが格納されているかを確認できます。データフレームが格納されています。これは主要な結果であり、確認してみましょう。最初の列には遺伝子が含まれ、2番目の列にはそれらが属するクラスター番号が含まれています。

# 'df'コンポーネントに何が格納されているか見てみましょう
head(expression_clusters$df)

1番目のグループに興味があるので、データフレームをフィルタリングしてそれらの遺伝子のみを保持します：

# グループ1の遺伝子を抽出
group_1 <- expression_clusters$df %>%
          filter(cluster == 1)

遺伝子グループを抽出した後、アノテーションパッケージを使用して追加情報を取得できます。また、これらの遺伝子リストを下流の機能分析ツールの入力として使用して、さらなる生物学的洞察を得、ゲノムが特定の機能を共有しているかどうかを確認することもできます。