edgeRドキュメント第二章（一部抜粋） 日本語訳
バイオサプリ2.1 用語
edgeRは、事前に定義されたゲノム特徴の差次的豊富さの分析を行います。厳密には必要ありませんが、通常、これらのゲノム特徴が重複しないことが望ましいです。単純化のために、今後はゲノム特徴を「遺伝子」と呼ぶことにしますが、原則としては転写産物、エクソン、一般的なゲノム区間、または他のタイプの特徴である可能性があります。ChIP-seq実験の場合、豊富さは転写因子の結合やヒストンマークの占有に関連しているかもしれませんが、今後は豊富さを遺伝子発現の観点から言及します。言い換えれば、このガイドの残りの部分はRNA-seq実験の遺伝子レベル分析の用語を使用しますが、方法論はそれよりも広く適用可能です。
2.2 ゲノムへのリードのアラインメント
RNA-seq分析の最初のステップは、通常、生のシーケンスリードを参照ゲノムにアラインすることですが、このプロセスには多くのバリエーションがあります。アラインメントは、リードが複数のエクソンにまたがっており、それらがゲノム上のよく分離された場所にアラインする可能性があることを考慮する必要があります。Subread-featureCountsパイプライン [24, 25] は非常に高速かつ効果的であると私たちは考えていますが、STAR-featureCounts、STAR-htseq、Bowtie-TopHat-htseqパイプラインも人気です。SubreadとfeatureCountsは、Bioconductor RパッケージRsubread [26] に実装されているため、特に便利です。
2.3 リードカウントのテーブルの作成
edgeRは、遺伝子に対応する行と独立したライブラリに対応する列があるリードカウントのテーブルで作業します。カウントは、各遺伝子（または他のゲノム位置）にアラインするリードの総数を示します。
このようなカウントは、アラインされたリードからさまざまなショートリードソフトウェアツールを使用して生成できます。私たちは、Rsubreadパッケージ [25, 26] のfeatureCounts関数が特に効果的で便利であると考えていますが、GenomicRangesパッケージのfindOverlapsやPythonソフトウェアのhtseq-countsなど、他のツールも利用可能です。
リードは、いくつかの方法でカウントできます。遺伝子レベルの分析を行う場合、カウントは遺伝子のゲノムスパン内のどこでもマッピングされるリードのものであるか、またはエクソンのみにカウントされる場合があります。通常、指定された遺伝子の任意のエクソンに重なるリードをカウントし、UTRを最初のエクソンの一部として含めます [25]。
プールされたshRNA-seqやCRISPR-Cas9遺伝子スクリーンからのデータの場合、processAmplicons関数 [11] を使用して、FASTQファイルから直接カウントを取得できます。edgeRは実際のリードカウントと共に動作するように設計されています。予測された転写産物の豊富さを実際のカウントの代わりにedgeRに入力することはお勧めしません。
2.4 ファイルからカウントを読む
カウントのテーブルがファイルに書き込まれている場合、通常、分析の最初のステップはこれらのカウントをRセッションに読み込むことになります。
カウントデータが複数の列（サンプルごとに1つの列）を持つ単一のタブ区切りまたはカンマ区切りのテキストファイルに含まれている場合、最も簡単な方法は、read.delimなどの標準的なR読み込み関数を使用してファイルをRに読み込むことです。詳細は上記のクイックスタートや、このガイドの後半にあるLNCaP細胞または口腔癌のケーススタディを参照してください。
異なるサンプルのカウントが別々のファイルに保存されている場合、ファイルを個別に読み込み、まとめる必要があります。このために、edgeRのreadDGE関数が提供されています。ファイルには、カウント用の1列と遺伝子識別子用の1列が含まれている必要があります。
2.5 疑似アラインメントと準マッピング
kallistoとSalmonのソフトウェアツールは、参照ゲノムの代わりにトランスクリプトームにシーケンスリードをアラインし、各トランスクリプトごとの推定カウントを生成します。いずれのツールの出力も、tximportパッケージを介してedgeRに入力でき、これにより遺伝子レベルの推定カウントと関連するedgeRオフセットマトリックスが生成されます。あるいは、トランスクリプトレベルの分析を行う意図がある場合は、catchSalmonおよびcatchKallisto関数を使用してkallistoまたはSalmonの出力を直接edgeRに読み込むことができます。
2.6 DGEListデータクラス
edgeRは、DGEListと呼ばれる単純なリストベースのデータオブジェクトにデータを格納します。このタイプのオブジェクトは、Rの任意のリストのように操作できるため、使いやすいです。readDGE関数は直接DGEListオブジェクトを作成します。カウントのテーブルがすでにマトリックスまたはデータフレーム（xと呼ぶ）として利用可能である場合、次のようにしてDGEListオブジェクトを作成できます：
y <- DGEList(counts=x)
同時にグループファクターを追加することもできます：
> group <- c(1,1,2,2)
> y <- DGEList(counts=x, group=group)
DGEListオブジェクトの主なコンポーネントは、整数カウントを含むマトリックスcounts、サンプルやライブラリに関する情報を含むデータフレームsamples、および遺伝子やゲノム特徴の注釈を含むオプショナルなデータフレームgenesです。データフレームsamplesには、各サンプルのライブラリサイズまたはシーケンシング深度を示すlib.size列が含まれています。ユーザーによって指定されていない場合、ライブラリサイズはカウントの列の合計から計算されます。classic edgeRでは、データフレームsamplesは、各サンプルのグループメンバーシップを識別するためにgroup列も含んでいる必要があります。
2.7 フィルタリング
全ライブラリで非常に低いカウント数の遺伝子は、差次的発現の証拠をほとんど提供しません。生物学的な観点からは、遺伝子がタンパク質に翻訳される可能性があるか、生物学的に重要であると見なされるためには、ある最低限の発現レベルが必要です。さらに、これらのカウントの顕著な離散性は、後のパイプラインで使用されるいくつかの統計的近似に干渉します。したがって、これらの遺伝子は、さらなる分析に進む前にフィルタリングする必要があります。
一般的なルールとして、すべてのサンプルで発現があり得ない遺伝子は除外します。ユーザーは自分自身で遺伝子の発現定義を設定できます。通常、遺伝子はそのライブラリで発現と見なされるために、カウントが5〜10である必要があります。ユーザーは、カウントそのものではなく、カウント・パー・ミリオン（CPM）でフィルタリングを行うべきです。後者はサンプル間のライブラリサイズの違いを考慮していないからです。
以下に簡単な例を示します。DGEListオブジェクトyのサンプル情報が次のように示されているとします：
y$samples
group lib.size norm.factors
Sample1 1 10880519 1
Sample2 1 9314747 1
Sample3 1 11959792 1
Sample4 2 7460595 1
Sample5 2 6714958 1
次のコマンドを使用して低発現遺伝子をフィルタリングします：
> keep <- filterByExpr(y)
> y <- y[keep, , keep.lib.sizes=FALSE]
filterByExpr 関数は、最小限のサンプル数で有用なカウントを持つ行を保持します（この場合は最小グループサイズが2なので2つのサンプル）。この関数は、yに含まれるグループ因子にアクセスして最小グループサイズを計算しますが、フィルタリングはどのサンプルがどのグループに属するかに依存せずに実行されるため、バイアスが導入されることはありません。フィルタリング後にDGEListオブジェクトのライブラリサイズを再計算することが推奨されますが、下流の分析はこれが行われるかどうかに対して頑健です。グループ因子や実験設計行列は、DGEListオブジェクトに既に設定されていない場合、直接 filterByExpr 関数に渡すこともできます：
keep <- filterByExpr(y, group=group)
より一般的には、filterByExpr を任意の設計行列とともに使用できます：
keep <- filterByExpr(y, design)
この形式では、設計行列は完全に一般的であり、連続共変量を含むこともできます。
フィルタリングは、差次的発現検定で関与するグルーピング因子や処理因子に基づくべきであり、科学的に自体が興味のないブロッキング変数に基づくべきではありません。たとえば、同じ治療が複数の被験者または患者に適用されるペア比較実験を考えてみてください：
design <- model.matrix(~ Patient + Treatment)
この設計では、患者間のベースラインの違いを補正するために患者が設計行列に含まれていますが、患者間での差次的発現をテストするわけではありません。したがって、フィルタリングは治療に基づくべきであり、患者に基づくべきではありません。すなわち：
keep <- filterByExpr(y, group=Treatment)
とし、次のようにするべきではありません：
keep <- filterByExpr(y, design)
2.8 正規化
2.8.1 正規化はサンプル固有の影響に対してのみ必要
edgeRは発現レベルの定量化よりも差次的発現分析に関心があります。条件間の発現レベルの相対的変化に関心がありますが、絶対的な発現レベルの推定には直接関心がありません。これにより、技術的要因を考慮する際に大幅に簡略化されます。なぜなら、実験条件に無関係な技術的要因は、差次的発現分析の結果に影響を与えるべきではないからです。たとえば、リードカウントは通常、各トランスクリプトの長さおよび発現に比例すると期待されますが、edgeRは一般に遺伝子長に対して調整を行う必要はありません。なぜなら、遺伝子長は各RNAサンプルのリードカウントに対して同じ相対的影響を持つからです。このため、正規化の問題は技術的要因がサンプル固有の影響を持つ場合にのみ発生します。
2.8.2 シーケンシング深度
遺伝子発現レベル以外でリードカウントに影響を与える最も明らかな技術的要因は、各RNAサンプルのシーケンシング深度です。edgeRは、異なるライブラリサイズによって表される異なるシーケンシング深度に対して、いかなる差次的発現分析も調整します。これは基本的なモデリング手続きの一部であり、自動的にfold-changeやp値計算に組み込まれます。常に存在し、ユーザーの介入は不要です。
2.8.3 有効ライブラリサイズ
差次的発現に対する二番目に重要な技術的影響は、より明白でないものです。RNA-seqは各RNAサンプル内の各遺伝子の相対的な豊富さを測定しますが、細胞あたりの総RNA出力の測定は提供しません。言い換えれば、RNA-seqは絶対的な発現ではなく、相対的な発現を測定します。これは、一部のサンプルで少数の遺伝子が非常に高い発現を示すが他のサンプルではそうでない場合に、差次的発現分析で重要になります。特定のサンプルにおいて、高発現の遺伝子がライブラリサイズのかなりの部分を消費すると、残りの遺伝子がそのサンプルで過少サンプリングされることになります。この影響が調整されない場合、残りの遺伝子がそのサンプルで誤ってダウンレギュレートされたように見える可能性があります。
normLibSizes 関数は、ほとんどの遺伝子に対してサンプル間のログフォールド変化を最小化するようにライブラリサイズを正規化します。このスケーリングファクターを計算するためのデフォルト方法は、サンプル間のTMM（トリム平均M値）です。元のライブラリサイズとスケーリングファクターの積を有効ライブラリサイズと呼びます。すなわち、正規化されたライブラリサイズです。有効ライブラリサイズは、すべての下流の分析で元のライブラリサイズに置き換わります。
TMMは、サンプルのペア間で過半数（半分以上）の遺伝子が差次的発現していないと考えられるRNA-Seqデータのほとんどで推奨されます。もし y がDGEListオブジェクトであるなら、以下のコマンドでTMM正規化を行い、正規化ファクターを表示します：
y <- normLibSizes(y)
y$samples
group lib.size norm.factors
Sample1 1 10880519 1.17
Sample2 1 9314747 0.86
Sample3 1 11959792 1.32
Sample4 2 7460595 0.91
Sample5 2 6714958 0.83
DGEListのすべての正規化ファクターのセットは、1に乗算して、効果的なライブラリサイズの幾何平均が元のライブラリサイズの幾何平均と同じになるようにします。正規化ファクターが1未満の場合、高カウントの少数の遺伝子がシーケンシングを独占し、他の遺伝子のカウントがライブラリサイズに対して通常よりも低くなります。その結果、ライブラリサイズは縮小され、そのライブラリ内のカウントが上方にスケーリングされるのと同じです。逆に、1より大きなファクターはライブラリサイズを拡大し、カウントを下方にスケーリングするのと同じです。
2.8.4 GC含量
各遺伝子のGC含量はサンプル間で変わらないため、差次的発現分析に対してほとんど影響がないと期待されます。しかし、最近の研究では、GC含量に対するサンプル固有の影響が検出されたことが示されています。EDASeqおよびcqnパッケージは、edgeRと互換性のある方法でサンプル固有のGC含量の影響を調整する補正ファクターを推定します。各ケースで、観察固有の補正ファクターをedgeRのglm関数にオフセット行列として入力することができます。
2.8.5 遺伝子長
GC含量と同様に、遺伝子長もサンプル間で変わらないため、差次的発現分析に対してほとんど影響がないと期待されます。それでも、遺伝子長に対するサンプル固有の影響が検出されていますが、その証拠はGC含量ほど強くありません。
2.8.6 モデルベースの正規化、変換ではない
edgeRでは、正規化は統計モデルに入る補正ファクターの形を取ります。このような補正ファクターは通常、edgeR関数によって内部的に計算されますが、ユーザーが提供することも可能です。補正ファクターは、normLibSizes で計算されたライブラリサイズのスケーリングファクターの形式を取ることがあり、それが有効ライブラリサイズを計算するために使用されます。あるいは、遺伝子固有の補正ファクターをedgeRのglm関数にオフセットとして入力することもできます。この場合、オフセット行列はシーケンシング深度やRNA組成を含むすべての正規化の問題を考慮していると見なされます。
edgeRでの正規化はモデルベースであり、元のリードカウント自体が変換されるわけではありません。つまり、ユーザーはedgeRに入力する前にリードカウントを変換するべきではありません。例えば、RPKMやFPKMの値をリードカウントの代わりにedgeRに入力するべきではありません。このような数量は、edgeRがデータの平均-分散関係を正しく推定するのを妨げます。これは、edgeRの統計的戦略にとって重要です。同様に、ユーザーはリードカウントに人工的な値を追加してからedgeRに入力するべきではありません。
edgeRは、たとえばCufflinksによって出力される推定発現レベルでの作業には設計されていません。edgeRは、RSEMによって出力される期待カウントで作業することができますが、生のカウントが好まれます。
2.8.7 擬似カウント
古典的なedgeR関数 estimateCommonDisp と exactTest は、出力オブジェクトの一部として擬似カウントの行列を生成します。擬似カウントは、散布推定と古典的edgeRパイプラインでの正確なテストの条件付き尤度の計算を迅速化するために内部的に使用されます。擬似カウントは、フィットしたモデルを前提とした場合に、ライブラリサイズがすべて等しかったと仮定した場合の相当するカウントを表します。擬似カウントは特定の目的のために計算されており、その計算は実験設計およびライブラリサイズに依存するため、ユーザーは擬似カウントを一般的な正規化カウントとして解釈しないことが推奨されます。擬似カウントは主にedgeRパイプラインでの内部使用を目的としています。
曖昧さの解消
他のソフトウェアパッケージでは、擬似カウントという用語を、edgeRの事前カウントと類似の意味で使用することがあります。つまり、カウントがゼロのときに対数計算を行う際に欠損値を回避するために追加される開始値です。edgeRでは、擬似カウントは正規化されたカウントの一種であり、事前カウントは小さなカウントを補うために使用される開始値です。
2.9 負の二項分布モデル
2.9.1 はじめに
RNA-Seq実験の出発点は、通常、さまざまな処理条件に関連付けられたn個のRNAサンプルです。各サンプルはシーケンスされ、短いリードが適切なゲノムにマッピングされ、各ゲノムの特徴にマッピングされたリードの数が記録されます。サンプルiから遺伝子gにマッピングされたリードの数はygiで表されます。サンプルiの遺伝子ごとのカウントのセットは、そのサンプルの発現プロファイルまたはライブラリを構成します。各カウントの期待されるサイズは、ライブラリサイズとそのサンプルでのその遺伝子の相対的な豊富さの積です。
2.9.2 生物学的変動係数 (BCV)
RNA-Seqプロファイルはn個のRNAサンプルから形成されます。πgiをサンプルi内の全cDNA断片の中で遺伝子gに由来する割合とします。Gを遺伝子の総数とすると、各サンプルで∑(g=1 to G) πgi = 1 です。√ϕgをレプリケートi間のπgiの変動係数 (CV)（標準偏差を平均で割ったもの）とします。ライブラリiのマッピングされたリードの総数をNi、g番目の遺伝子にマッピングされたリードの数をygiとします。
この場合、
E(ygi) = µgi = Niπgi です。
同じRNAサンプルの反復シーケンス実行でカウントygiがポアソン分布に従うと仮定すると、混合分布の分散に関するよく知られた公式が示唆するのは、
var(ygi) = Eπ [var(y|π)] + varπ [E(y|π)] = µgi + ϕgµgi^2 です。
両辺をµgi^2 で割ると、
CV^2(ygi) = 1/µgi + ϕg です。
最初の項1/µgiはポアソン分布の二乗CVであり、第二項は未観測の発現値の二乗CVです。したがって、総合的なCV^2は、πgiが測定される技術的なCV^2と、真のπgiの生物学的CV^2の合計です。この記事では、ϕgを分散と呼び、p√ϕgを生物学的CVと呼びますが、厳密には、これはレプリケート間のライブラリ間の変動のすべての要因を捉え、ライブラリ準備などの技術的な原因からの寄与やサンプル間の真の生物学的変動も含む可能性があります。
RNA-Seq実験では、2つのレベルの変動が区別できます。まず、各遺伝子の相対的な豊富さはRNAサンプル間で変動し、主に生物学的な原因によるものです。次に、測定誤差があり、各サンプルでの各遺伝子の豊富さの推定に対する不確実性です。同じRNAサンプルのアリコートがシーケンスされる場合、特定の遺伝子のリードカウントはポアソン法則に従うべきです。シーケンスの変動がポアソン分布に従うなら、各生物学的レプリケートライブラリ間のカウントの二乗変動係数 (CV) は、技術的および生物学的変動のそれぞれの二乗CVの合計であることが示されます。
Total CV^2 = Technical CV^2 + Biological CV^2
生物学的CV (BCV) は、レプリケートRNAサンプル間で（未知の）真の遺伝子豊富さがどの程度変動するかを示す変動係数です。これは、シーケンシングの深さを無限に増加させた場合に生物学的レプリケート間に残るであろうCVを表します。技術的CVはカウントのサイズが増加するにつれて減少しますが、BCVはそうではありません。したがって、BCVは高カウント遺伝子にとって主要な不確実性の源である可能性が高いため、RNA-Seq実験での差次的発現の現実的な評価にはBCVの信頼できる推定が重要です。もし、各遺伝子の豊富さがレプリケートRNAサンプル間で、遺伝子ごとの標準偏差が遺伝子ごとの平均に比例するように変動するなら、これは物理的な量の測定によく見られる特性であり、BCVが遺伝子間でおおむね一定であると考えるのが合理的です。ただし、BCVが遺伝子間で変動し、遺伝子発現や予想カウントに対して体系的な傾向を示す可能性があることも考慮します。
BCVの大きさは、真の遺伝子豊富さが従う確率論的な法則の正確さよりも重要です。数学的な便宜のために、真の遺伝子豊富さがレプリケートRNAサンプル間でガンマ分布の法則に従うと仮定します。これにより、リードカウントが負の二項分布確率法則に従うことが示唆されます。
2.9.3 BCVの推定
負の二項分布モデルが適合されるときには、分析を実行する前にBCVを推定する必要があります。BCVは、負の二項分布モデル下での分散パラメーターの平方根であるため、負の二項分布モデルの分散を推定するのと同等です。
シーケンシングデータの並列性により、各遺伝子個別の推論を助けるために遺伝子の集合から情報を借りるいくつかの可能性があります。遺伝子間で情報を共有する最も簡単な方法は、すべての遺伝子が同じ平均-分散関係を持っていると仮定すること、つまり分散がすべての遺伝子で同じであるとすることです。
この「共通分散」アプローチの拡張は、分散関数内のパラメーターに平均依存の傾向を置くことであり、同じ期待カウントを持つすべての遺伝子が同じ分散を持つようにします。
しかし、遺伝子発現レベルは遺伝子間で非同一で依存する分布を持つため、上記の仮定は単純すぎる可能性があります。遺伝子ごとの分散関数を許可し、経験ベイズの調整を用いたより一般的なアプローチが数年前に導入され、一般化線形モデル（GLM）やより複雑な実験デザインに拡張されました。タグごとの分散を使用する場合、レプリケート間で一貫性のある遺伝子が、一貫性のない遺伝子よりも高くランク付けされます。多くのRNA-Seqデータセットでは、遺伝子固有の分散を許可することが、差次的発現が外れ値によって駆動されることを防ぐために必要です。したがって、タグごとの分散はモデル適合と差次的発現の検定において強く推奨されます。
edgeRでは、タグごとの分散をグローバル分散傾向や共通分散値に向けて絞り込むために、経験ベイズ戦略を適用します。絞り込みの量は、グローバル値に対する重みとタグごとの推定の精度によって決まります。これらの戦略は、通常、比較的少数のレプリケートがある場合に有用です。
2.10 edgeRの古典的なパイプライン：二群以上のグループ間のペアワイズ比較
2.10.1 分散の推定
edgeRは、単一因子の実験に対して、分位点調整条件付き最大尤度（qCML）法を使用します。広範なシミュレーション研究において、いくつかの他の推定量（例えば、最大尤度推定量、準尤度推定量など）と比較した結果、qCMLは広範な条件下でバイアスに対して最も信頼性が高く、特に多くの小さなサンプルで共通の分散を持つ状況において最も優れた性能を発揮します。このモデルは、次世代シーケンシングデータに適用可能です。DNAシーケンシングのコストが高いため、SAGEおよびRNA-seq実験で多数のレプリケートを行うことができないため、非常に小さなサンプルに焦点を当てています。
qCML法は、各タグの総カウントに条件付けして尤度を計算し、ライブラリサイズの調整後に擬似カウントを使用します。カウントのテーブルまたはDGEListオブジェクトがある場合、estimateDisp()関数を使用してqCMLの共通分散およびタグごとの分散を推定できます。別の方法として、estimateCommonDisp()関数を使用してqCMLの共通分散を推定し、次にestimateTagwiseDisp()関数を使用してqCMLのタグごとの分散を推定することもできます。
ただし、qCML法は単一因子設計のデータセットにのみ適用可能で、複雑な実験における複数因子の効果を考慮できません。実験に複数の因子が関与する場合、新しい分散推定方法を探す必要があります。
以下は、qCML法を使用して分散を推定する簡単な例です。DGEListオブジェクトyがある場合、次のコマンドを使用して分散を推定します。
共通分散とタグごとの分散を一度に推定する（推奨）:
y <- estimateDisp(y)
別々に推定する場合、共通分散を推定する:
y <- estimateCommonDisp(y)
次に、タグごとの分散を推定する:
y <- estimateTagwiseDisp(y)
共通分散を推定した後でタグごとの分散を推定する必要があることに注意してください。
2.10.2 DEG遺伝子の検定
次世代シーケンシングデータ解析において、最も関心があるのは、二群（またはそれ以上の群）間での差次的に発現する遺伝子/タグを見つけることです。負の二項分布モデルが適合し、分散の推定が得られたら、正確検定を使用して差次的発現を判断するための検定手順に進むことができます。
正確検定は、qCML法に基づいています。グループ内のカウントの合計に対する条件付き分布を知ることで、帰無仮説下での観測されたカウントの合計の確率よりも小さい確率を持つすべてのカウントの合計を合計することで、正確なp値を計算できます。ネガティブ・バイノミアル分布に対する正確な検定は、フィッシャーの正確検定と強い類似性があります。
前のセクションで説明したように、正確検定は単一因子の実験にのみ適用されます。検定はexactTest()関数を使用して行い、関数は共通分散アプローチとタグごとの分散アプローチの両方をサポートしています。例えば:
et <- exactTest(y)
topTags(et)
2.11 より複雑な実験（GLM機能）
2.11.1 一般化線形モデル
一般化線形モデル（GLM）は、非正規分布に従う応答データに対する古典的線形モデルの拡張です。GLMは、平均-分散の関係に応じて確率分布を指定します。例えば、上記で指定したリードカウントの二次平均-分散関係です。推定値がϕgの値に対して利用可能であれば、分散は任意のµgiの値について評価できます。GLM理論を用いて、次のようなログ線形モデルを適合させます:
log µgi = x^T_i βg + log Ni
ここで、xiはRNAサンプルiに適用された処置条件を指定する共変量のベクトルであり、βgは遺伝子gの共変量効果を媒介する回帰係数のベクトルです。二次分散関数はネガティブ・バイノミアルGLMの分布ファミリーを指定します。ネガティブ・バイノミアル分布の使用は、πgiをガンマ分布として扱うことに相当します。
2.11.2 分散の推定
一般的な実験（複数因子が関与する場合）では、edgeRはCox-Reidプロファイル調整尤度（CR）法を使用して分散を推定します。CR法は、前述のqCML法の限界を克服するために導出されました。CR法は、設計行列を用いて一般化線形モデル（GLM）を適合させることで複数因子を処理します。
CR法は、近似条件付き尤度のアイデアに基づいています。カウントのテーブルまたはDGEListオブジェクトと実験の設計行列が与えられた場合、一般化線形モデルを適合させます。これにより、すべての系統的変動源が考慮され、分散を有効に推定することができます。
CR法は、すべてのタグに対する共通分散、タグの豊富さに応じたトレンド分散、または個別タグの分散を計算するのに使用できます。これらは、指定された設計でestimateDisp()関数を呼び出すことで実行できます。別の方法として、estimateGLMCommonDisp()、estimateGLMTrendedDisp()、およびestimateGLMTagwiseDisp()を使用して、共通分散、トレンド分散、およびタグごとの分散を別々に推定できます。タグごとの分散アプローチは、多因子実験の場合に強く推奨されます。
以下は、GLM法を使用して分散を推定する簡単な例です。DGEListオブジェクトyと設計行列がある場合、次のコマンドを使用して分散を推定します。
共通分散、トレンド分散、タグごとの分散を一度に推定する（推奨）:
y <- estimateDisp(y, design)
別々に推定する場合、共通分散を推定する:
y <- estimateGLMCommonDisp(y, design)
トレンド分散を推定する:
y <- estimateGLMTrendedDisp(y, design)
タグごとの分散を推定する:
y <- estimateGLMTagwiseDisp(y, design)
タグごとの分散を推定する前に、共通分散またはトレンド分散のいずれかを推定する必要があることに注意してください。タグごとの分散を推定する際には、経験的ベイズ法が適用され、タグごとの分散が共通分散またはトレンド分散に向けて調整されます。両方が存在する場合、デフォルトではトレンド分散が使用されます。
より詳細な例については、セクション4.1（Tuchのデータ）、セクション4
2.12 複製がない場合の対処法
edgeRは、主に生物学的複製を含むデータでの使用を想定しています。しかし、RNA-SeqやChIP-Seqは高コストの技術であるため、各処理条件に対して一つのライブラリしか作成できないことがあります。この場合、生物学的変動を推定するための複製ライブラリが存在しません。このような状況では、データアナリストは以下の選択肢に直面しますが、どれも理想的ではありません。これらの選択肢を生物学的複製の満足できる代替手段として推奨するものではありません。むしろ、分析段階でできる最善の方法であり、選択肢2～4は生物学的変動がないと仮定するよりは良いかもしれません。
-
記述的な分析に満足する
MDSプロットやフォールド変化の分析を含むことがありますが、有意性の分析は試みない方が良いでしょう。これが最善のアドバイスかもしれません。 -
経験に基づいて妥当な分散値を選ぶ
類似のデータに基づいて妥当な分散値を選び、それをexactTestやglmFitに使用します。よく制御された実験から得られるデータセットの一般的な共通BCV（平方根分散）値は、人間のデータで0.4、遺伝的に同一のモデル生物のデータで0.1、技術的複製で0.01です。以下はシミュレーションデータの例です：
bcv <- 0.2
counts <- matrix( rnbinom(40,size=1/bcv^2,mu=10), 20,2)
y <- DGEList(counts=counts, group=1:2)
et <- exactTest(y, dispersion=bcv^2)
選択した分散値に対して得られるp値と有意な遺伝子の数は非常に敏感であり、バッチ効果などの未考慮の要因がある場合、実際にはここで示されたよりもはるかに大きな分散が存在する可能性があります。それにもかかわらず、生物学的変動を完全に無視するよりも、名目上の分散値を選ぶ方が現実的かもしれません。
-線形モデルから一つまたは複数の説明変数を削除する
これにより、いくつかの残差自由度が得られます。理想的には、最も重要でない因子を削除することですが、もし因子が一つでグループが二つだけの場合、これは設計行列全体を削除するか、インターセプト用の一つの列に減らすことを意味します。実験に複数の説明因子がある場合は、最も小さなフォールド変化を持つ因子を削除することができます。複数の処理条件がある場合は、最も似ている二つの条件を複製として扱うことを試みることができます。この縮小モデルから分散を推定し、これらの分散を含む完全な設計行列を持つデータオブジェクトに挿入し、glmFitやglmLRTを使ってモデルの適合とテストを進めます。このアプローチは、DE遺伝子の数が比較的小さい場合にのみ成功します。
この縮小設計行列とともに、estimateGLMCommonDisp（method="deviance", robust=TRUE, subset=NULL）を試すことができます。これは複製なしで分散を推定するための現在の最良の試みですが、カウントがあまりにも小さくなく、DE遺伝子が全体の少数である場合にのみ良好な結果を得られます。これは使えるものを返すための最良の試みであることをご理解ください。信頼できる分散の推定には通常、複製が必要です。
- DEでないべきコントロール転写産物が存在する場合、その分散を推定する
例えば、コントロール遺伝子として実験に使用された処理に応答しないハウスキーピング遺伝子を示すインデックス変数があるとします。まず、一つの処理グループだけのデータオブジェクトのコピーを作成します：
y1 <- y
y1$samples$group <- 1
次に、ハウスキーピング遺伝子とすべてのライブラリを一つのグループとして共通の分散を推定します：
y0 <- estimateDisp(y1[housekeeping,], trend="none", tagwise=FALSE)
これを完全なデータオブジェクトに挿入して続けます：
y$common.dispersion <- y0$common.dispersion
fit <- glmFit(y, design)
lrt <- glmLRT(fit)
十分に大きな数のコントロール転写産物が必要です。少なくとも数十個、理想的には数百個が必要です。
2.13 フォールド変化の閾値を超える差次的発現
上記のすべての検定方法は、差次的発現を統計的有意性に基づいて特定しますが、差がどれほど小さいかには関係ありません。一方で、特定の量だけ発現レベルが変化する遺伝子に興味がある場合があります。一般的に使用されるアプローチは、DEテストを実施し、フォールド変化のカットオフを適用し、その後、そのフォールド変化の閾値を超えるすべての遺伝子をp値でランク付けする方法です。場合によっては、まずp値カットオフに基づいて遺伝子を選び、その後、フォールド変化でソートします。これらのp値とフォールド変化の閾値の組み合わせは、どちらか一方を単独で使用するよりも、生物学的に意味のある遺伝子セットを提供するようです。しかし、これらはどちらもアドホックであり、フォールド変化の閾値に対する差次的発現のテストに対して意味のあるp値を提供しません。これらは、低発現で変動が大きい遺伝子を好み、一般的にFDRの制御を破壊します。
edgeRは、GLMフレームワークの下で閾値付き仮説に対する厳密な統計テストを提供します。これはTREAT [31] に類似していますが、元のTREATメソッドよりもはるかに強力です。フォールド変化（または対数フォールド変化）閾値が与えられた場合、glmTreat() 関数を使用してDGEGLMオブジェクトで閾値付きテストを実施できます。この関数は、glmFit() または glmQLFit() によって生成されたDGEGLMオブジェクトに対して呼び出します。
以下は例として、1対2のlog2-フォールド変化が1よりも有意に大きい遺伝子を検出する場合のコマンドです：
fit <- glmQLFit(y, design)
tr <- glmTreat(fit, coef=2, lfc=1)
topTags(tr)
glmTreat() におけるフォールド変化の閾値は、テスト結果から予期される最小値ではありません。遺伝子が統計的に有意とされるためには、この閾値を超える必要があります。この閾値は「私たちが遺伝子に対して興味がない下限のフォールド変化」と解釈する方が良いです。DE遺伝子が非常に多く存在する場合は、関心のある遺伝子に検索を絞るために、比較的大きなフォールド変化の閾値が適切かもしれません。一方で、DE遺伝子が存在しない場合は、小さなフォールド変化の閾値、あるいは閾値なしの使用を検討します。
詳細な例については、セクション4.4（Fuのデータ）のケーススタディを参照してください。
原文