Summary
RNAシーケンシングの微分発現解析方法の詳細なプロトコルが提供されました: リンマ, エッジ, DESeq2.
Abstract
RNAシーケンシング(RNA-seq)は、遺伝子改変と複雑な生物学的プロセスとの関係を明らかにし、腫瘍の診断、予後、および治療に大きな価値を持つため、トランスクリプトミクスで最も広く使用されている技術の1つです。RNA-seqデータの微分分析は異常転写を同定するために非常に重要であり、リンマ、EdgeRおよびDESeq2は微分分析のための効率的なツールである。しかし、RNA-seq微分解析には、R言語を持つ特定のスキルと、医学教育のカリキュラムに欠けている適切な方法を選択する能力が必要です。
本明細書において、リンマ、DESeq2およびEdgeRを介して、脈管癌(CHOL)と正常組織の間で、それぞれ微分発現遺伝子(DEG)を同定するための詳細なプロトコルを提供し、その結果は火山プロットおよびベン図に示される。リンマ、DESeq2およびEdgeRの3つのプロトコルは類似しているが、分析のプロセス間で異なったステップを有する。たとえば、線形モデルは limma の統計に使用され、負の二項分布は edgeR および DESeq2 で使用されます。さらに、正規化された RNA-seq カウント データは、EdgeR と limma に必要ですが、DESeq2 には必要ありません。
ここでは、リンマ、EdgeR、DESeq2の3つの差分解析方法の詳細なプロトコルを提供します。3つの方法の結果は部分的に重なっています。3 つのメソッドには独自の利点があり、メソッドの選択はデータのみに依存します。
Introduction
RNA-シーケンシング(RNA-seq)は、多くの利点(例えば、高いデータ再現性)を有するトランスクリプトミクスで最も広く使用されている技術の1つであり、複雑な生物学的プロセス1,2の機能とダイナミクスについての理解を劇的に高めています。異なる生物学的文脈下での異常転写物の同定は、RNA-seq分析の重要なステップである。RNA-seqは、病因と関連する分子機構と生物学的機能を深く理解することを可能にします。したがって、微分分析は、腫瘍3、4、5の診断、予後および治療に対して価値があると考えられてきた。現在、RNA-seqの微分発現解析、特にリンマ、DESeq2およびEdgeR1、6、7のために、よりオープンソースのR/バイオ伝導体パッケージが開発されている。しかし、差分分析には、R言語を持つ特定のスキルと、医学教育のカリキュラムに欠けている適切な方法を選択する能力が必要です。
このプロトコルでは、がんゲノムアトラス(TCGA)から抽出された胆管癌(CHOL)RNA-seqカウントデータに基づいて、最も知られた方法の3つ(リンマ8、EdgeR9およびDESeq210)を、それぞれRプログラム11によって実施し、CHOLと正常組織との間のDEGsを同定した。リンマ、EdgeR、DESeq2の3つのプロトコルは似ていますが、分析のプロセス間で異なるステップがあります。例えば、正規化された RNA-seq カウントデータは、EdgeR および limma8,9に必要ですが、DESeq2 は独自のライブラリーの不一致を使用して、正規化10ではなくデータを修正します。さらに、edgeRはRNA-seqデータに特に適しており、リムマはマイクロアレイおよびRNA-seqに使用されます。線形モデルは、DEG12を評価するためにlimmaによって採用され、edgeRの統計は経験的ベイズ推定、正確なテスト、一般化線形モデルおよび準尤度検定を含む負の二項分布に基づいている。
要約すると、それぞれリンマ、DESeq2、EdgeRを用いたRNA-seq微分発現解析の詳細なプロトコルを提供します。この記事を参照することにより、ユーザーは簡単にRNA-seqの微分分析を行い、データに適切な微分分析方法を選択することができます。
Subscription Required. Please recommend JoVE to your librarian.
Protocol
注:R-studioプログラムを開き、Rファイル"DEGs.R"をロードし、ファイルは補助ファイル/スクリプトから取得することができます。
1. データのダウンロードと前処理
- がんゲノムアトラス(TCGA)からコリンジオカルシノーマ(CHOL)のハイスループットシーケンシング(HTSeq)カウントデータをダウンロードしてください。このステップは、次の R コードで簡単に実現できます。
- R パッケージをインストールするには、[ 実行 ] をクリックします。
- [ 実行 ] をクリックして R パッケージを読み込みます。
if(!必要な名前空間(「バイオックマネージャー」、静かに=真))
+ インストール.パッケージ(「バイオックマネージャー」)
バイオクマネージャー::インストール(c(「TCGAbiolinks」、"要約実験") - 作業ディレクトリを設定します。
図書館(TCGAbiolinks)
ライブラリ(要約実験)
setwd(「C:/ユーザー/リウシイ/デスクトップ」) - がんの種類を選択します。
がん< - "TCGA-CHOL" - "GDCquery.R" ファイルから R コードを実行して、データをダウンロードします。ファイル "GDCquery.R" は補助ファイル/スクリプトから取得できます。
ソース(「補足ファイル/スクリプト/GDCquery.R」)
ヘッド(cnt)
##TCGA-3X-AAVA-01A-11R-A41I-07
##ENSG00000000003 4262
##ENSG00000000005 1
##ENSG00000000419 1254
##ENSG00000000457 699
##ENSG00000000460 239
##ENSG00000000938 334
注: 実行後、CHOLHTSeq カウント データはダウンロードされ、"cnt" という名前が付けられ、行はアンサンブル遺伝子 ID を表し、列はサンプル ID を表します。サンプル ID の 14 から 15 桁目の数字に注意してください。01から09までの数字は腫瘍を示し、10から19の範囲は正常組織を示す。
- アンサンブル遺伝子 ID を遺伝子シンボルに変換します。
- 注釈ファイルをストレージ パスに従って R にインポートします。注釈ファイル (gencode.v22.annotation.gtf) は補助ファイルから取得できます。
gtf_v22 <- rtracklayer::インポート('補足ファイル/gencode.v22.annotation.gtf') - 「gtf_v22から R コードを実行します。補足ファイル/スクリプトから取得できる R" ファイル:
ソース("補足ファイル/スクリプト/gtf_v22。R") - アンサンブル遺伝子IDを遺伝子シンボルに変換する関数「ann」を適用します。
cnt=ann(cnt,gtf_v22)
- 注釈ファイルをストレージ パスに従って R にインポートします。注釈ファイル (gencode.v22.annotation.gtf) は補助ファイルから取得できます。
- 低発現遺伝子のフィルタリング
- R パッケージ "edgeR" をインストールするには、[ 実行 ] をクリックします。
バイオクマネージャー::インストール("エッジR") - R パッケージ 「edgeR」をロードするには、[ 実行 ] をクリックします。
ライブラリ(エッジR) - 次の R コードを実行して、少なくとも 2 つのサンプルで 1 を超える 100 万個あたりの数 (CPM) 値を持つ遺伝子を保持します。
<行を保持するsums(cpm(cnt)>1)>=2
<- as.matrix(cnt[keep,])
注: 100 万分の数 (CPM) 値は、異なるシーケンス深度によって生じる偏差を排除するために、読み取り数の代わりに使用されます。
- R パッケージ "edgeR" をインストールするには、[ 実行 ] をクリックします。
2. 「リンマ」による微分発現解析
- R パッケージ "limma" をインストールするには、[ 実行 ] をクリックします。
バイオクマネージャー::インストール(「リンマ」) - R パッケージ 「リンマ」、「エッジR」をロードするには、[ 実行 ] をクリックします。
ライブラリ(リンマ)
ライブラリ(エッジR) - 次の R コードを実行して、デザイン マトリックスを作成します。
群 <- substring(colnames(cnt),14,15) # Extract group information
グループ [グループ %in% "01"] <- "Cancer" # set '01' as tumor tissue
グループ [グループ %in% "11"] <- "Normal" # set '11' as normal tissue
群 <- factor (group, levels = c("Normal","Cancer"))- 設計マトリックスを作成します。
デザイン < - モデル.マトリックス (~グループ)
行名(デザイン) <- コルネーム(cnt) - DGEList オブジェクトを作成します。
dge <- DGEList(カウント = cnt、グループ = グループ) - データを正規化します。
dge <- calcNormFactors(dge,メソッド = "TMM") - 次の R コードを実行して、リンマトレンドメソッドベースの微分式分析を実行します。
dge
##An クラス "DGEList" のオブジェクト
##$counts
##TCGA-3X-AAVA-01A-11R-A41I-07
##TSPAN6 4262
##DPM1 1254
##SCYL3 699
##C1orf112 239
##FGR 334 - CPM 値を計算します。
<- cpm(dge、ログ =TRUE、前の.count=3) - [ 実行 ]をクリックして線形モデルに適合し、データを予測するか、変数間の関係を推測します。
フィット< - lmFit (ログ、デザイン) - ベイジアンに基づいて T 値、F 値、および対数オッズを計算します。
フィット< - eBayes(フィット、トレンド=真) - 結果表を抽出します。
res_limma< - データとして.フレーム(トップテーブル(フィット、n=Inf))
頭(res_limma)
## ログFC AveExpr t P.値の adj.P.ヴァルB
##RP11-252E2.2 -4.899493 -2.488589 -20.88052 2.386656e-25 4.931786e-21 47.28823
##BX842568.1 -4.347930 -2.595205 -20.14532 1.082759e-24 1.118706e-20 45.83656
##CTC-537E7.3 -5.154894 -2.143292 -19.59571 3.452354e-24 2.216114e-20 44.72001
##RP11-468N14.3 -6.532259 -2.029714 -19.49409 4.289807e-24 2.216114e-20 44.51056
##AP006216.5 -4.507051 -2.670915 -19.25649 7.153356e-24 2.956339e-20 44.01704
##RP11-669E14.4 -4.107204 -2.828311 -18.93246 1.448209e-23 4.987633-20 43.33543
#Theの結果、遺伝子ID、log2折りたたみ値(logFC)、実験中の遺伝子の平均log2発現量(AveExpr)、改変t統計(t)、再発明p値(P.Value)を含む「res_limma」に保存され、偽発見率(FDR)補正されたp値(adj)。P.Val)と微分発現遺伝子の対数オッズ(B)
注: "edgeR" 関数は、サンプルの準備やライブラリの構築とシーケンスによって引き起こされる影響を排除するためにデータを正規化するために使用されました。設計マトリックスの構築では、マトリックスのサンプルIDに実験計画(例えば、組織タイプ:正常または腫瘍組織)を一致させる必要があります。リンマトレンドは、シーケンシング深度が同じデータに適していますが、リンマ・ブームは適しています: (i) サンプルライブラリのサイズが異なる場合。(ii) TMMによって正規化されていないデータ。(iii) データに多くのノイズがあります。正のlogFCは、遺伝子が実験でアップレギュレートされることを意味し、負の数は遺伝子がダウンレギュレートされていることを意味します。 - DEG を識別します。
res_limma$sig <- as.factor()
ifelse(res_limma$adj。P.Val < 0.05と腹筋(res_limma$logFC)> 2,
ifelse(res_limma$logFC > 2 ,''up'')))#adj.p値は0.05と|log2FC<|>= 2 は、DEG を識別するためのしきい値です。
サマリー(res_limma $sig)
##down上がっていません
##1880 17341 1443 - 結果表をファイルに出力します。
書き込み.csv(res_limma、ファイル = 'result_limma.csv') - R パッケージ "ggplot2" をインストールするには、[ 実行 ] をクリックします。
パッケージをインストールする(「ggplot2」) - R パッケージ "ggplot2" をロードするには、[ 実行 ] をクリックします。
ライブラリ(ggplot2) - 「火山」から R コードを実行します。火山プロットを作成する場合は R" を選択します。ファイル "火山。R"は補助ファイルから取得できます。
ソース("補足ファイル/スクリプト/火山。R")
火山(res_limma"logFC","adj.P.Val",2,0.05)
注: 遺伝子は log2FC と adj-p 値に従って異なる位置にマッピングすることができ、アップレギュレーション DEG は赤で色分けされ、ダウンレギュレーション DEG は緑色で色分けされます。 - [ エクスポート ]をクリックして、火山プロットを保存します。
注:火山プロットは、生成され、異なる形式(例えば、pdf、TIFF、PNG、JPEG形式)でダウンロードすることができます。遺伝子は、log2FCとadj p 値に応じて異なる位置にマッピングすることができ、アップレギュレーションDEG(log2FC > 2、adj p < 0.05)は赤で着色され、ダウンレギュレートされたDEGs(log2FC < -2、adj p < 0.05)は緑色で着色され、非DEGは灰色で着色されます。
- 設計マトリックスを作成します。
3. 「edgeR」による微分表現解析
- R パッケージ 「edgeR」をロードするには、[ 実行 ] をクリックします。
ライブラリ(エッジR) - 次の R コードを実行して、デザイン マトリックスを作成します。
グループ<サブストリング(コルネーム(cnt),14,15)
グループ [グループ %in% "01"] < - "がん"
グループ [グループ %in% "11"] <- "通常"
グループ=因子(グループ、レベル= c(「正常」「癌」)
デザイン <-モデル.マトリックス(~グループ)
行名(デザイン) = コルネーム(cnt) - DGEList オブジェクトを作成するには、[ 実行 ] をクリックします。
< - DGEリスト(カウント=cnt) - データを正規化します。
dge <- calcNormFactors(dge,メソッド = "TMM") - [ 実行 ] をクリックして、遺伝子発現値の分散を推定します。
dge <- 見積Disp(dge、設計、堅牢性= T) - [ 実行 ] をクリックして、データをカウントするモデルに適合します。
フィット< - glmQLFit(dge、デザイン) - 統計的検定を実施します。
フィット< - glmQLFTest(フィット) - 結果表を抽出します。結果は、ログフォールド変更値、ログCPM、F、p値、およびFDR補正p値を含む「res_edgeR」に保存されます。
res_edgeR=as.data.frame(フィット、n=インフ))
頭(res_edgeR)
## ログFC ログ CPM F P値 FDR
##GCDH -3.299633 5.802700 458.5991 1.441773e-25 2.979280e-21
##MSMO1 -3.761400 7.521111 407.0416 1.730539e-24 1.787993e-20R
##CL1 -3.829504 5.319641 376.5043 8.652474e-24 5.516791e-20
##ADI1 -3.533664 8.211281 372.6671 1.067904e-23 5.516791e-20
##KCNN2 -5.583794 3.504017 358.6525 2.342106e-23 9.679455e-20
##GLUD1 -3.287447 8.738080 350.0344 3.848408e-23 1.194406e-19
#The結果は、ログ折りたたみ変更値(logFC)、ログCPM、F、p値、およびFDR補正p値を含む「res_edgeR」に保存されます。 - DEG を識別します。
res_edgeR$sig = as.factor(
ifelse(res_edgeR$FDR < 0.05と腹筋(res_edgeR $logFC) > 2,
ifelse(res_edgeR$logFC > 2 ,'up','down')、そうではありません)
サマリー(res_edgeR$sig)
##down上がっていません
##1578 15965 3121 - 結果表をファイルに出力します。
書き込み.csv(res_edgeR、ファイル = 'res_edgeR.csv') - 火山プロットを作成します。
火山(res_edgeR"logFC","FDR",2,0.05) - [ エクスポート ]をクリックして、火山プロットを保存します。
4. 「DESeq2」による微分発現解析
- R パッケージ "DESeq2" をインストールするには、[ 実行 ] をクリックします。
バイオクマネージャー::インストール(「DESeq2」) - R パッケージ "DESeq2" を読み込むには、[ 実行 ] をクリックします。
ライブラリ(DESeq2) - 次の R コードを実行して、グループ化係数を決定します。
グループ<サブストリング(コルネーム(cnt),14,15)
グループ [グループ %in% "01"] < - "がん"
グループ [グループ %in% "11"] <- "通常"
グループ=因子(グループ、レベル= c(「正常」「癌」) - オブジェクトを作成します。
dds <-DESeqDataSet マトリックス (cnt, データフレーム(グループ), デザイン = ~グループ)
dds
##class: デセクデータセット
##dim: 20664 45
##metadata(1): バージョン
##assays(1): カウント
##rownames(20664): TSPAN6 DPM1 ..RP11-274B21.13 LINC01144
##rowData名(0):
##colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ..
##colData名(1): グループ - 解析を実行します。
dds < - DESeq(dds) - 結果表を生成します。
res_DESeq2 < - データフレーム(結果(dds))
頭(res_DESeq2)
## ベースMean ログ2フォールド変更 lfcSE スタット p値パジ
##TSPAN6 4704.9243 -0.8204515 0.3371667 -2.433370 1.495899e-02 2.760180e-02
##DPM1 1205.9087 -0.3692497 0.1202418 -3.070894 2.134191e-03 4.838281e-03
##SCYL3 954.9772 0.2652530 0.2476441 1.071106 2.841218e-01 3.629059e-01
##C1orf112 277.7756 0.7536911 0.2518929 2.992109 2.770575e-03 6.101584e-03
##FGR 345.8789 -0.6423198 0.3712729 -1.730047 8.362180e-02 1.266833e-01
##CFH 27982.3546 -3.8761382 0.5473363 -7.081823 1.422708e-12 1.673241e-11
注: 結果は、正規化された読み取りカウント(baseMean)、ログフォールド変更値(log2FoldChange)、ログフォールド変更標準エラー(lfcSE)、Wald統計(stat)、元のp値(pvalue)、修正されたp値(padj)の平均を含む「res_DESeq2」に保存されます。 - DEG を識別します。
res_DESeq2$sig = as.factor()
ifelse(res_DESeq2$padj < 0.05 & abs(res_DESeq2$log2FoldChange) > 2,
ifelse(res_DESeq2$log2Foldchange > 2 ,'up','down')、そうではありません)
サマリー(res_DESeq2$sig)
##down上がっていません
##1616 16110 2938 - 結果表をファイルに出力します。
書き込み.csv(res_DESeq2、ファイル = 'res_DESeq2.csv') - 火山プロットを作成します。
火山(res_DESeq2"log2FoldChange"""padj",2,0.05) - [ エクスポート ]をクリックして、火山プロットを保存します。
5. ベン図
- R パッケージ "VennDiagram" をインストールするには、[ 実行 ] をクリックします。
パッケージをインストールする(「ベンダイアグラム」) - R パッケージ "VennDiagram" を読み込むには、[ 実行 ] をクリックします。
図書館 (ベンダイアグラム) - 規制されたDEGのベン図を作成します。
グリッド.新ページ()
グリッド.ドロー(ベンダイアグラム(リスト(リンマ=ローネーム(res_
リムマ[res_limma$シグ="アップ"、])
エッジR=行名(res_edgeR[res_edgeR$sig="up",))
DESeq2=行名(res_DESeq2[res_DESeq2$sig==
"アップ",[])
NULL、高さ = 3、幅 = 3、単位 = "in"、
col="黒"、lwd=0.3,fill=c(「#FF6666""、#FFFF00"、
"#993366"),
アルファ=c(0.5、0.5、0.5)、メイン=「アップ規制されたDEG」) - [ エクスポート ] をクリックしてベン図を保存します。
- ダウン規制されたDEGsのベン図を作成します。
グリッド.新ページ()
グリッド.ドロー(ベンダイアグラム(リスト(リンマ=ローネーム(res_
リンマ[res_limma$sig="down",))
エッジR=行名(res_edgeR[res_edgeR$sig==
「ダウン」、])
DESeq2=行名(res_DESeq2[res_DESeq2$sig="down",))
NULL、高さ = 3、幅 = 3、単位 = "in"、
col="黒"、lwd=0.3,fill=c(「#FF6666""、#FFFF00"、
"#993366"),
アルファ=c(0.5、0.5、0.5)、メイン=「ダウン規制されたDEG」) - [ エクスポート ] をクリックしてベン図を保存します。
Subscription Required. Please recommend JoVE to your librarian.
Representative Results
火山プロットとベン図が特に使用されている微分表現解析の結果を可視化するための様々なアプローチがあります。リンマは、cholと正常組織の間の3323 AGを|logFC|≥2およびadjで同定した。P.Val <0.05 しきい値として、その中で1880はCHOL組織でダウンレギュレートされ、1443 はアップレギュレートされた (図 1a)。一方、edgeR は 1578 のダウンレギュレーション DEG と 3121 のアップレギュレーション DEG を特定しました (図 1b)。DESeq2 は、1616 のダウンレギュレーション DEG と 2938 のアップレギュレーション DEG を特定しました (図 1c)。これら3つの方法の結果を比較すると、1431個のアップレギュレーションDEGと1531個のダウンレギュレートDEGが重なっていた(図2)。
図 1.CHOLと正常組織との間での遺伝子(DEG)の遺伝子の同定 (a-c)リンマ、エッジR、DESeq2によって獲得されたすべての遺伝子の火山プロットは、それぞれ 、adj p 値(-log10)がフォールド変化(log2)に対してプロットされ、赤色の点は、アップレギュレーションされたDEG(調整された p 値<0.05およびログ|FC|> 2)と緑の点は、ダウンレギュレーションされた DEG (調整された p 値< 0.05 およびログ|FC|< 2)。 この図の大きなバージョンを表示するには、ここをクリックしてください。
図 2.ベン図は、リンマ、エッジRおよびDESeq2から導出された結果の間で重複を示す。この図の大きなバージョンを表示するには、ここをクリックしてください。
Subscription Required. Please recommend JoVE to your librarian.
Discussion
癌における豊富な異常性転写物は、RNA-seq差動分析5によって容易に同定することができる。しかし、RNA-seq微分発現解析の適用は、R言語と適切な方法を選択する能力を持つ特定のスキルを必要とするため、しばしば制限されています。この問題に対処するために、我々は、3つの最も知られている方法(リンマ、EdgeRおよびDESeq2)とRNA-seqの微分発現解析を適用するためのチュートリアルの詳細な紹介を提供する。これにより、3つの方法すべてにおける類似点と相違点の理解が容易になり、個々のデータに適した方法を選択することができ、複雑な動的生物学的プロセスを理解できるようになります。
ここでは、リンマ、edgeR、DESeq2を通じたRNA-seq微分発現解析の詳細なプロトコルを、(i)データのダウンロードと前処理、(ii-iv)リンマ、edgeR、DESeq2を介した微分発現解析、(v)それぞれ、これらの3つの方法の結果をベン図を通じた比較に示す。
3つの方法は、微分式解析のプロセス間で類似したステップと異なるステップを有する。リニアモデルは、マイクロアレイ、RNA-セク、定量PCR8、13を含むすべての遺伝子発現技術に適用可能なリンマの統計に使用され、edgeRとDESeq2は負の二項分布9、10、およびedgeRおよびDESeq2に基づく統計方法論の範囲を実装するRNA-seqデータに適しています。さらに、正規化された RNA-seq カウント データは EdgeR と limma に必要ですが、DESeq2 は独自のライブラリの不一致を使用して正規化の代わりにデータを修正し、DESeq2 のデータは整数行列でなければなりません。正規化法には、TMM(M値のトリム平均)、TMMwsp、RLE(相対ログ式)およびアッパー四分位数が含まれ、その中でTMMはRNA-seqデータに最も一般的に使用される正規化方法である。3つの方法の結果は、DESeq2およびEdgeRがリンマよりも多くのDEGを得ることを示した。この違いの理由は、edgeR と DESeq2 が負の二項モデルに基づいているため、多数の誤検出に寄与します。逆に、limma-voomは分散関数のみを使用し、過剰な偽陽性を示さないため、リニアモデル解析に続くリニアモデル解析に続いて、リニアモデル解析がリムマ14、15、16になります。
3 つのメソッドには独自の利点があり、選択はデータの種類に依存します。例えば、マイクロアレイデータがある場合、リムマは優先して与えられるべきですが、次世代シーケンシングデータの場合、DESeq2とEdgeRは9、10、17が好ましい。要約すると、Rパッケージリンマ、エッジRおよびDESeq2をそれぞれ用いたRNA-seq微分発現解析のための詳細なプロトコルをここで提供する。3つの方法の出力結果は部分的に重なっており、これらの微分法にはそれぞれの利点があります。残念ながら、このプロトコルは、他のデータ型(例えば、マイクロアレイデータ)およびメソッド(例えば、EBSeq)18の技術的な詳細をカバーしていない。
Subscription Required. Please recommend JoVE to your librarian.
Disclosures
原稿は以前に出版されておらず、他の場所で出版される必要はありません。すべての著者は、重要な知的コンテンツのためにこの原稿の作成に貢献し、最終的な原稿を読んで承認しました。私たちは、利益相反がないことを宣言します。
Acknowledgments
この研究は、中国国立自然科学財団(グラント第81860276)と国家主要研究開発プログラム(グラント2018YFC1003200)の主要特別基金プロジェクトによって支援されました。
Materials
Name | Company | Catalog Number | Comments |
R | version 3.6.2 | free software | |
Rstudio | free software |
References
- Tambonis, T., Boareto, M., Leite, V. B. P. Differential Expression Analysis in RNA-seq Data Using a Geometric Approach. Journal of Computational Biology. 25, 1257-1265 (2018).
- Wang, Z., Gerstein, M., Snyder, M. RNA-Seq: a revolutionary tool for transcriptomics. Nature Reviews. Genetics. 10, 57-63 (2009).
- Anders, S., et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nature Protocols. 8, 1765-1786 (2013).
- McDermaid, A., Monier, B., Zhao, J., Liu, B., Ma, Q. Interpretation of differential gene expression results of RNA-seq data: review and integration. Briefings in Bioinformatics. 20, 2044-2054 (2019).
- Costa-Silva, J., Domingues, D., Lopes, F. M. RNA-Seq differential expression analysis: An extended review and a software tool. PloS One. 12, 0190152 (2017).
- Law, C. W., et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research. 5, (2016).
- Varet, H., Brillet-Guéguen, L., Coppée, J. Y., Dillies, M. A. SARTools: A DESeq2- and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data. PloS One. 11, 0157022 (2016).
- Ritchie, M. E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 43, 47 (2015).
- Robinson, M. D., McCarthy, D. J., Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 26, Oxford, England. 139-140 (2010).
- Love, M. I., Huber, W., Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 15, 550 (2014).
- Gentleman, R. C., et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biology. 5, 80 (2004).
- Law, C. W., Chen, Y., Shi, W., Smyth, G. K. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology. 15, 29 (2014).
- Smyth, G. K. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology. 3, (2004).
- Lund, S. P., Nettleton, D., McCarthy, D. J., Smyth, G. K. Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates. Statistical Applications in Genetics and Molecular Biology. 11, (2012).
- Reeb, P. D., Steibel, J. P. Evaluating statistical analysis models for RNA sequencing experiments. Frontiers in Genetics. 4, 178 (2013).
- Rocke, D. M., et al. Excess False Positive Rates in Methods for Differential Gene Expression Analysis using RNA-Seq Data. bioRxiv. , (2015).
- Agarwal, A., et al. Comparison and calibration of transcriptome data from RNA-Seq and tiling arrays. BMC genomics. 11, 383 (2010).
- Leng, N., et al. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 29, Oxford, England. 1035-1043 (2013).