Name: Three Differential Expression Analysis Methods for RNA Sequencing: limma, EdgeR, DESeq2 (Video) | JoVE Journal
Uploaded: 2021-09-18T14:00:00.000+00:00
Duration: 10 min 10 s
Description: Learn differential expression analysis for RNA sequencing using limma, EdgeR, and DESeq2 in RStudio. This protocol compares methods for identifying differentially expressed genes in Cholangiocarcinoma data.

ai_agent_IP

ai_agent

maintenance

instantTrialDurationMinutes

isBiopharma

enableNewMyJoveSectionDesign

educationFeatures

seoMetaApi

ai_copilot_institution_blacklist

ai_chat_institution_blacklist

freeAccessWindow

ssoRoleCapture

scheduledMaintenanceBanner

joveLabsFeatures

disable_autoplay_for_bots

extended_trial_duration

5 minutes

This work was supported by the&#8194;National&#8194;Natural&#8194;Science&#8194;Foundation&#8194;of&#8194;China&#8194;(Grant&#8194;No. 81860276) and Key Special Fund Projects of National Key R&#38;D Program (Grant&#8194;No. 2018YFC1003200).

NOTE: Open the R-studio program and load R file &#34;DEGs.R&#34;, the file can be acquired from Supplementary files/Scripts.
1. Downloading and pre-processing of data
<ol>
	<li>Download the high-throughput sequencing (HTSeq) count data of cholangiocarcinoma (CHOL) from The Cancer Genome Atlas (TCGA). This step can be easily achieved by the following R code.
	<ol>
		<li>Click Run to install R packages.</li>
		<li>Click Run to load R packages. 
		if(!requireNamespace(&#34;BiocManager&#34;, quietly=TRUE)) 
		&#160;+ install.packages(&#34;BiocManager&#34;) 
		BiocManager::install(c(&#34;TCGAbiolinks&#34;, &#34;SummarizedExperiment&#34;))</li>
		<li>Set the working directory. 
		library (TCGAbiolinks) 
		library(SummarizedExperiment) 
		setwd(&#34;C:/Users/LIUSHIYI/Desktop&#34;)</li>
		<li>Choose the cancer type. 
		cancer &#60;- &#34;TCGA-CHOL&#34;</li>
		<li>Run the R code from the &#34;GDCquery.R&#34; file to download the data. The file &#34;GDCquery.R&#34; can be acquired from Supplementary files/Scripts: 
		source(&#34;Supplementary files/Scripts/GDCquery.R&#34;) 
		head(cnt) 
		##TCGA-3X-AAVA-01A-11R-A41I-07 
		##ENSG00000000003 4262 
		##ENSG00000000005 1 
		##ENSG00000000419 1254 
		##ENSG00000000457 699 
		##ENSG00000000460 239 
		##ENSG00000000938 334 
		NOTE: After execution, the CHOLHTSeq count data will be downloaded and named &#34;cnt&#34;, where rows represent ensemble gene IDs and columns represent sample IDs. Please notice the numbers at positions 14-15 in the sample IDs; numbers ranging from 01 to 09 indicate tumors and ranging from 10 to 19 indicate normal tissues.</li>
	</ol>
	</li>
	<li>Convert ensemble gene IDs to gene symbols.
	<ol>
		<li>Import the annotation file into R according to its storage path. The annotation file (gencode.v22.annotation.gtf) can be acquired from Supplementary files. 
		&#8203;gtf_v22 &#60;- rtracklayer::import(&#39;Supplementary files/gencode.v22.annotation.gtf&#39;)</li>
		<li>Run the R code from the &#34;gtf_v22.R&#34; file, which can be acquired from Supplementary files/Scripts: 
		source(&#34;Supplementary files/Scripts/gtf_v22.R&#34;)</li>
		<li>Apply the function &#34;ann&#34; to convert the ensemble gene IDs to gene symbols. 
		cnt=ann(cnt,gtf_v22)</li>
	</ol>
	</li>
	<li>Filtering low-expressed genes
	<ol>
		<li>Click Run to install the R package &#34;edgeR&#34;. 
		&#8203;BiocManager::install(&#34;edgeR&#34;)</li>
		<li>Click Run to load the R package &#34;edgeR&#34;. 
		&#8203;library(edgeR)</li>
		<li>Run the following R code to keep genes with counts per million (CPM) values greater than one in at least two samples. 
		keep &#60;- rowSums(cpm(cnt)&#62;1)&#62;=2 
		cnt &#60;- as.matrix(cnt[keep,]) 
		&#8203;NOTE: The counts per million (CPM) value is used instead of the read count to eliminate the deviation caused by different sequencing depths.</li>
	</ol>
	</li>
</ol>
2. Differential expression analysis through &#34;limma&#34;
<ol>
	<li>Click Run to install the R package &#34;limma&#34;. 
	BiocManager::install(&#34;limma&#34;)</li>
	<li>Click Run to load the R packages &#34;limma&#34;, &#34;edgeR&#34;. 
	library(limma) 
	library(edgeR)</li>
	<li>Run the following R code to create the design matrix. 
	group &#60;- substring(colnames(cnt),14,15) # Extract group information 
	group [group %in% &#34;01&#34;] &#60;- &#34;Cancer&#34; # set &#39;01&#39; as tumor tissue 
	group [group %in% &#34;11&#34;] &#60;- &#34;Normal&#34; # set &#39;11&#39; as normal tissue 
	&#8203;group &#60;- factor (group, levels = c(&#34;Normal&#34;,&#34;Cancer&#34;))
	<ol>
		<li>Create the design matrix. 
		design &#60;- model.matrix (~group) 
		rownames(design) &#60;- colnames(cnt)</li>
		<li>Create the DGEList object. 
		dge &#60;- DGEList(counts = cnt, group = group)</li>
		<li>Normalize the data. 
		dge &#60;- calcNormFactors(dge, method = &#34;TMM&#34;)</li>
		<li>Run the following R code to perform the limma-trend method based differential expression analysis. 
		dge 
		##An object of class &#34;DGEList&#34; 
		##$counts 
		##TCGA-3X-AAVA-01A-11R-A41I-07 
		##TSPAN6 4262 
		##DPM1 1254 
		##SCYL3 699 
		##C1orf112 239 
		##&#8203;FGR 334</li>
		<li>Calculate the CPM value. 
		&#8203;logdge &#60;- cpm(dge, log=TRUE, prior.count=3)</li>
		<li>Click Run to fit a linear model to predict the data or infer the relationship between variables. 
		fit &#60;- lmFit (logdge, design)</li>
		<li>Calculate the T value, F value and log-odds based on Bayesian. 
		fit &#60;- eBayes(fit, trend=TRUE)</li>
		<li>Extract the result table. 
		res_limma&#60;- as.data.frame(topTable(fit,n=Inf)) 
		 
		head(res_limma) 
		## logFC AveExpr t P.Value adj.P.Val 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.987633e-20 43.33543 
		#The result of differential expression analysis is saved in &#34;res_limma&#34;, which includes the gene id, log2 fold change value (logFC), the average log2 expression level of the gene in the experiment (AveExpr), the modified t statistic (t), relavent p value (P.Value), the false discovery rate (FDR) corrected p value (adj.P.Val) and the log-odds of differentially expressed genes (B) 
		&#8203;NOTE: The function &#34;calcNormFactors()&#34; of the &#34;edgeR&#34; was used to normalize the data to eliminate the influence caused by sample preparation or library construction and sequencing. In the construction of design matrix, it is necessary to match experimental design (e.g., tissue type: normal or tumor tissues) to sample IDs of the matrix. limma-trend is suitable to data whose sequencing depth is the same, while limma-voom is suitable: (i) when the sample library size is different; (ii) data not normalized by TMM; (iii) there is a lot of &#34;noise&#34; in the data. A positive logFC means that gene is up-regulated in the experiment, while negative number means that gene is down-regulated.</li>
		<li>Identify the DEGs. 
		res_limma$sig &#60;- as.factor( 
		&#160;ifelse(res_limma$adj.P.Val &#60; 0.05 &#38; abs(res_limma$logFC) &#62; 2, 
		&#160; &#160;ifelse(res_limma$logFC &#62; 2 ,&#39;up&#39;,&#39;down&#39;),&#39;not&#39;)) # The adj.p Value &#60; 0.05 and the |log2FC| &#62;= 2 are thresholds to identify the DEGs 
		summary(res_limma$sig) 
		##down not up 
		##1880 &#8203;17341 1443</li>
		<li>Output the result table to a file. 
		write.csv(res_limma, file = &#39;result_limma.csv&#39;)</li>
		<li>Click Run to install the R package &#34;ggplot2&#34;. 
		install.packages(&#34;ggplot2&#34;)</li>
		<li>Click Run to load the R package &#34;ggplot2&#34;. 
		library(ggplot2)</li>
		<li>Run the R code from the &#34;volcano.R&#34; to create the volcano plot. The file &#34;volcano.R&#34; can be acquired from Supplementary files. 
		source(&#34;Supplementary files/Scripts/volcano.R&#34;) 
		volcano(res_limma,&#34;logFC&#34;,&#34;adj.P.Val&#34;,2,0.05) 
		NOTE: Genes can be mapped to different positions according to their log2FC and adj-p values, the up regulated DEGs are colored in red, and the down-regulated DEGs are colored in green.</li>
		<li>Click Export to save the volcano plot. 
		&#8203;NOTE: The volcano plots can be generated and downloaded in different formats (e.g., pdf, TIFF, PNG, JPEG format). Genes can be mapped to different positions according to their log2FC and adj p values, the up-regulated DEGs (log2FC &#62; 2, adj p &#60; 0.05) are colored in red, and the down-regulated DEGs (log2FC &#60; -2, adj p &#60; 0.05) are colored in green, non-DEGs are colored in grey.</li>
	</ol>
	</li>
</ol>
3. Differential expression analysis through &#34;edgeR&#34;
<ol>
	<li>Click Run to load the R package &#34;edgeR&#34;. 
	library(edgeR)</li>
	<li>Run the following R code to create design matrix. 
	group &#60;-substring(colnames(cnt),14,15) 
	group [group %in% &#34;01&#34;] &#60;- &#34;Cancer&#34; 
	group [group %in% &#34;11&#34;] &#60;- &#34;Normal&#34; 
	group=factor(group, levels = c(&#34;Normal&#34;,&#34;Cancer&#34;)) 
	design &#60;-model.matrix(~group) 
	rownames(design) = colnames(cnt)</li>
	<li>Click Run to create the DGEList object. 
	dge &#60;- DGEList(counts=cnt)</li>
	<li>Normalize the data. 
	dge &#60;- calcNormFactors(dge, method = &#34;TMM&#34;)</li>
	<li>Click Run to estimate the dispersion of gene expression values. 
	dge &#60;- estimateDisp(dge, design, robust = T)</li>
	<li>Click Run to fit model to count data. 
	fit &#60;- glmQLFit(dge, design)</li>
	<li>Conduct a statistical test. 
	fit &#60;- glmQLFTest(fit)</li>
	<li>Extract the result table. The result is saved in &#34;res_edgeR&#34;, which includes the log fold change value, log CPM, F, p value and FDR corrected p value. 
	res_edgeR=as.data.frame(topTags(fit, n=Inf)) 
	head(res_edgeR) 
	## logFC logCPM F PValue 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 result is saved in &#34;res_edgeR&#34;, which includes the log fold change value(logFC), log CPM, F, p value and FDR corrected p value</li>
	<li>Identify the DEGs. 
	res_edgeR$sig = as.factor( 
	ifelse(res_edgeR$FDR &#60; 0.05 &#38; abs(res_edgeR$logFC) &#62; 2, 
	ifelse(res_edgeR$logFC &#62; 2 ,&#39;up&#39;,&#39;down&#39;),&#39;not&#39;)) 
	summary(res_edgeR$sig) 
	##down not up 
	##1578 15965 3121</li>
	<li>Output the result table to a file. 
	write.csv(res_edgeR, file = &#39;res_edgeR.csv&#39;)</li>
	<li>Create the volcano plot. 
	volcano(res_edgeR,&#34;logFC&#34;,&#34;FDR&#34;,2,0.05)</li>
	<li>Click Export to save the volcano plot.</li>
</ol>
4. Differential expression analysis through &#34;DESeq2&#34;
<ol>
	<li>Click Run to install R packages &#34;DESeq2&#34;. 
	BiocManager::install(&#34;DESeq2&#34;)</li>
	<li>Click Run to load R packages &#34;DESeq2&#34;. 
	library(DESeq2)</li>
	<li>Run the following R code to determine the grouping factor. 
	group &#60;-substring(colnames(cnt),14,15) 
	group [group %in% &#34;01&#34;] &#60;- &#34;Cancer&#34; 
	group [group %in% &#34;11&#34;] &#60;- &#34;Normal&#34; 
	group=factor(group, levels = c(&#34;Normal&#34;,&#34;Cancer&#34;))</li>
	<li>Create the DESeqDataSet object. 
	dds &#60;-DESeqDataSetFromMatrix (cnt, DataFrame(group), design = ~group) 
	dds 
	##class: DESeqDataSet 
	##dim: 20664 45 
	##metadata(1): version 
	##assays(1): counts 
	##rownames(20664): TSPAN6 DPM1 ... RP11-274B21.13 LINC01144 
	##rowData names(0): 
	##colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ... 
	##colData names(1): group</li>
	<li>Perform the analysis. 
	dds &#60;- DESeq(dds)</li>
	<li>Generate the result table. 
	res_DESeq2 &#60;- data.frame(results(dds)) 
	 
	head(res_DESeq2) 
	## baseMean log2FoldChange lfcSE stat pvalue padj 
	##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 
	NOTE: The result is saved in &#34;res_DESeq2&#34;, which includes the mean of the normalized read count (baseMean), log fold Change value(log2FoldChange), log fold change standard error (lfcSE), the Wald statistic (stat), original p value (pvalue) and corrected p value (padj)</li>
	<li>Identify DEGs. 
	res_DESeq2$sig = as.factor( 
	&#160;ifelse(res_DESeq2$padj &#60; 0.05 &#38; abs(res_DESeq2$log2FoldChange) &#62; 2, 
	&#160; &#160;ifelse(res_DESeq2$log2FoldChange &#62; 2 ,&#39;up&#39;,&#39;down&#39;),&#39;not&#39;)) 
	summary(res_DESeq2$sig) 
	##down not up 
	##1616 16110 2938</li>
	<li>Output the result table to a file. 
	write.csv(res_DESeq2, file = &#39;res_DESeq2.csv&#39;)</li>
	<li>Create the volcano plot. 
	volcano(res_DESeq2,&#34;log2FoldChange&#34;,&#34;padj&#34;,2,0.05)</li>
	<li>Click Export to save the volcano plot.</li>
</ol>
5. Venn diagram
<ol>
	<li>Click Run to install the R package &#34;VennDiagram&#34;. 
	install.packages(&#34;VennDiagram&#34;)</li>
	<li>Click Run to load the R package &#34;VennDiagram&#34;. 
	library (VennDiagram)</li>
	<li>Make a Venn diagram of up regulated DEGs. 
	grid.newpage() 
	grid.draw(venn.diagram(list(Limma=rownames(res_ 
	limma[res_limma$sig==&#34;up&#34;,]), 
	&#160; edgeR=rownames(res_edgeR[res_edgeR$sig==&#34;up&#34;,]), 
	&#160; DESeq2=rownames(res_DESeq2[res_DESeq2$sig== 
	&#160; &#34;up&#34;,])), 
	&#160;NULL,height = 3,width = 3,units = &#34;in&#34;, 
	&#160;col=&#34;black&#34;,lwd=0.3,fill=c(&#34;#FF6666&#34;,&#34;#FFFF00&#34;, 
	&#160;&#34;#993366&#34;), 
	&#160;alpha=c(0.5, 0.5, 0.5),main = &#34;Up-regulated DEGs&#34;))</li>
	<li>Click Export to save the Venn diagram.</li>
	<li>Make a Venn diagram of down regulated DEGs. 
	grid.newpage() 
	grid.draw(venn.diagram(list(Limma=rownames(res_ 
	limma[res_limma$sig==&#34;down&#34;,]), 
	&#160; &#160;edgeR=rownames(res_edgeR[res_edgeR$sig== 
	&#34;down&#34;,]), 
	&#160; &#160;DESeq2=rownames(res_DESeq2[res_DESeq2$sig==&#34;down&#34;,])), 
	&#160;NULL,height = 3,width = 3,units = &#34;in&#34;, 
	&#160;col=&#34;black&#34;,lwd=0.3,fill=c(&#34;#FF6666&#34;,&#34;#FFFF00&#34;, 
	&#34;#993366&#34;), 
	&#160;alpha=c(0.5, 0.5, 0.5),main = &#34;Down-regulated DEGs&#34;))</li>
	<li>Click Export to save the Venn diagram.</li>
</ol>

<ol><li>Tambonis, T., Boareto, M., Leite, V. B. P. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Differential+Expression+Analysis+in+RNA-seq+Data+Using+a+Geometric+Approach%5BTitle%5D)+AND+%22Journal+of+Computational+Biology%22%5BJournal%5D)">Differential Expression Analysis in RNA-seq Data Using a Geometric Approach.</a> Journal of Computational Biology. 25, 1257-1265 (2018).</li>
<li>Wang, Z., Gerstein, M., Snyder, M. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((RNA-Seq%3A+a+revolutionary+tool+for+transcriptomics%5BTitle%5D)+AND+%22Nature+Reviews.+Genetics%22%5BJournal%5D)">RNA-Seq: a revolutionary tool for transcriptomics.</a> Nature Reviews. Genetics. 10, 57-63 (2009).</li>
<li>Anders, S., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Count-based+differential+expression+analysis+of+RNA+sequencing+data+using+R+and+Bioconductor%5BTitle%5D)+AND+%22Nature+Protocols%22%5BJournal%5D)">Count-based differential expression analysis of RNA sequencing data using R and Bioconductor.</a> Nature Protocols. 8, 1765-1786 (2013).</li>
<li>McDermaid, A., Monier, B., Zhao, J., Liu, B., Ma, Q. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Interpretation+of+differential+gene+expression+results+of+RNA-seq+data%3A+review+and+integration%5BTitle%5D)+AND+%22Briefings+in+Bioinformatics%22%5BJournal%5D)">Interpretation of differential gene expression results of RNA-seq data: review and integration.</a> Briefings in Bioinformatics. 20, 2044-2054 (2019).</li>
<li>Costa-Silva, J., Domingues, D., Lopes, F. M. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((RNA-Seq+differential+expression+analysis%3A+An+extended+review+and+a+software+tool%5BTitle%5D)+AND+%22PloS+One%22%5BJournal%5D)">RNA-Seq differential expression analysis: An extended review and a software tool.</a> PloS One. 12, 0190152(2017).</li>
<li>Law, C. W., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((RNA-seq+analysis+is+easy+as+1-2-3+with+limma%2C+Glimma+and+edgeR%5BTitle%5D)+AND+%22F1000Research%22%5BJournal%5D)">RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR.</a> F1000Research. 5, (2016).</li>
<li>Varet, H., Brillet-Guéguen, L., Coppée, J. Y., Dillies, M. A. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((SARTools%3A+A+DESeq2-+and+EdgeR-Based+R+Pipeline+for+Comprehensive+Differential+Analysis+of+RNA-Seq+Data%5BTitle%5D)+AND+%22PloS+One%22%5BJournal%5D)">SARTools: A DESeq2- and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data.</a> PloS One. 11, 0157022(2016).</li>
<li>Ritchie, M. E., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((limma+powers+differential+expression+analyses+for+RNA-sequencing+and+microarray+studies%5BTitle%5D)+AND+%22Nucleic+Acids+Research%22%5BJournal%5D)">limma powers differential expression analyses for RNA-sequencing and microarray studies.</a> Nucleic Acids Research. 43, 47(2015).</li>
<li>Robinson, M. D., McCarthy, D. J., Smyth, G. K. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((edgeR%3A+a+Bioconductor+package+for+differential+expression+analysis+of+digital+gene+expression+data%5BTitle%5D)+AND+%22Bioinformatics%22%5BJournal%5D)">edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.</a> Bioinformatics. 26, Oxford, England. 139-140 (2010).</li>
<li>Love, M. I., Huber, W., Anders, S. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Moderated+estimation+of+fold+change+and+dispersion+for+RNA-seq+data+with+DESeq2%5BTitle%5D)+AND+%22Genome+Biology%22%5BJournal%5D)">Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.</a> Genome Biology. 15, 550(2014).</li>
<li>Gentleman, R. C., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Bioconductor%3A+open+software+development+for+computational+biology+and+bioinformatics%5BTitle%5D)+AND+%22Genome+Biology%22%5BJournal%5D)">Bioconductor: open software development for computational biology and bioinformatics.</a> Genome Biology. 5, 80(2004).</li>
<li>Law, C. W., Chen, Y., Shi, W., Smyth, G. K. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((voom%3A+Precision+weights+unlock+linear+model+analysis+tools+for+RNA-seq+read+counts%5BTitle%5D)+AND+%22Genome+Biology%22%5BJournal%5D)">voom: Precision weights unlock linear model analysis tools for RNA-seq read counts.</a> Genome Biology. 15, 29(2014).</li>
<li>Smyth, G. K. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Linear+models+and+empirical+bayes+methods+for+assessing+differential+expression+in+microarray+experiments%5BTitle%5D)+AND+%22Statistical+Applications+in+Genetics+and+Molecular+Biology%22%5BJournal%5D)">Linear models and empirical bayes methods for assessing differential expression in microarray experiments.</a> Statistical Applications in Genetics and Molecular Biology. 3, (2004).</li>
<li>Lund, S. P., Nettleton, D., McCarthy, D. J., Smyth, G. K. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Detecting+differential+expression+in+RNA-sequence+data+using+quasi-likelihood+with+shrunken+dispersion+estimates%5BTitle%5D)+AND+%22Statistical+Applications+in+Genetics+and+Molecular+Biology%22%5BJournal%5D)">Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates.</a> Statistical Applications in Genetics and Molecular Biology. 11, (2012).</li>
<li>Reeb, P. D., Steibel, J. P. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Evaluating+statistical+analysis+models+for+RNA+sequencing+experiments%5BTitle%5D)+AND+%22Frontiers+in+Genetics%22%5BJournal%5D)">Evaluating statistical analysis models for RNA sequencing experiments.</a> Frontiers in Genetics. 4, 178(2013).</li>
<li>Rocke, D. M., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Excess+False+Positive+Rates+in+Methods+for+Differential+Gene+Expression+Analysis+using+RNA-Seq+Data%5BTitle%5D)+AND+%22bioRxiv%22%5BJournal%5D)">Excess False Positive Rates in Methods for Differential Gene Expression Analysis using RNA-Seq Data.</a> bioRxiv. , (2015).</li>
<li>Agarwal, A., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Comparison+and+calibration+of+transcriptome+data+from+RNA-Seq+and+tiling+arrays%5BTitle%5D)+AND+%22BMC+genomics%22%5BJournal%5D)">Comparison and calibration of transcriptome data from RNA-Seq and tiling arrays.</a> BMC genomics. 11, 383(2010).</li>
<li>Leng, N., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((EBSeq%3A+an+empirical+Bayes+hierarchical+model+for+inference+in+RNA-seq+experiments%5BTitle%5D)+AND+%22Bioinformatics%22%5BJournal%5D)">EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments.</a> Bioinformatics. 29, Oxford, England. 1035-1043 (2013).</li>
</ol>

The manuscript has not been published before and is not being considered for publication elsewhere. All authors have contributed to the creation of this manuscript for important intellectual content and read and approved the final manuscript. We declare there is no conflict of interest.

Abundant aberrate transcripts in cancers can be easily identified by RNA-seq differential analysis5. However, the application of RNA-seq differential expression analysis is often restricted as it requires certain skills with R language and the capacity to choose appropriate methods. To address this problem, we provide a detailed introduction to the three most known methods (limma, EdgeR and DESeq2) and tutorials for applying the RNA-seq differential expression analysis. This will facilitate the un...

RNA-sequencing (RNA-seq) is one of the most widely used technologies in transcriptomics with many advantages (e.g., high data reproducibility), and has dramatically increased our understanding of the functions and dynamics of complex biological processes1,2. Identification of aberrate transcripts under different biological context, which are also known as differentially expressed genes (DEGs), is a key step in RNA-seq analysis. RNA-seq makes it possible to get a deep understanding of pathogenesis related molecular mechanisms and biological functions. Therefore, differential analysis has been regarded as valuable for diagnostics, prognostics and therapeutics of tumors3,4,5. Currently, more open-source R/Bioconductor packages have been developed for RNA-seq differential expression analysis, particularly limma, DESeq2 and EdgeR1,6,7. However, differential analysis requires certain skills with R language and the ability to choose the appropriate method, which is lacking in the curriculum of medical education.
In this protocol, based on the cholangiocarcinoma (CHOL) RNA-seq count data extracted from The Cancer Genome Atlas (TCGA), three of the most known methods (limma8, EdgeR9 and DESeq210) were carried out, respectively, by the R program11 to identify the DEGs between CHOL and normal tissues. The three protocols of limma, EdgeR and DESeq2 are similar but have different steps among the processes of the analysis. For example, the normalized RNA-seq count data is necessary for EdgeR and limma8, 9, whereas DESeq2 uses its own library discrepancies to correct data instead of normalization10. Furthermore, edgeR is specifically suitable for RNA-seq data, while the limma is used for microarrays and RNA-seq. A linear model is adopted by limma to assess the DEGs12, while the statistics in edgeR are based on the negative binomial distributions, including empirical Bayes estimation, exact tests, generalized linear models and quasi-likelihood tests9.
In summary, we provide the detailed protocols of RNA-seq differential expression analysis by using limma, DESeq2 and EdgeR, respectively. By referring to this article, users can easily perform the RNA-seq differential analysis and choose the appropriate differential analysis methods for their data.

<table><tbody><tr><td>R</td><td></td><td>version 3.6.2</td><td>free software</td></tr><tr><td>Rstudio</td><td></td><td></td><td>free software</td></tr></tbody></table>

three differential expression analysis methods for rna sequencing: limma, edger, deseq2

Three Differential Expression Analysis Methods for RNA Sequencing: limma, EdgeR, DESeq2

RNA sequencing (RNA-seq) is one of the most widely used technologies in transcriptomics as it can reveal the relationship between the genetic alteration and complex biological processes and has great value in diagnostics, prognostics, and therapeutics of tumors. Differential analysis of RNA-seq data is crucial to identify aberrant transcriptions, and limma, EdgeR and DESeq2 are efficient tools for differential analysis. However, RNA-seq differential analysis requires certain skills with R language and the ability to choose an appropriate method, which is lacking in the curriculum of medical education.
Herein, we provide the detailed protocol to identify differentially expressed genes (DEGs) between cholangiocarcinoma (CHOL) and normal tissues through limma, DESeq2 and EdgeR, respectively, and the results are shown in volcano plots and Venn diagrams. The three protocols of limma, DESeq2 and EdgeR are similar but have different steps among the processes of the analysis. For example, a linear model is used for statistics in limma, while the negative binomial distribution is used in edgeR and DESeq2. Additionally, the normalized RNA-seq count data is necessary for EdgeR and limma but is not necessary for DESeq2.
Here, we provide a detailed protocol for three differential analysis methods: limma, EdgeR and DESeq2. The results of the three methods are partly overlapping. All three methods have their own advantages, and the choice of method only depends on the data.

There are various approaches to visualize the result of differential expression analysis, among which the volcano plot and Venn diagram are particularly used. limma identified 3323 DEGs between the CHOL and normal tissues with the |logFC|≥2 and adj.P.Val &lt;0.05 as thresholds, among which 1880 were down-regulated in CHOL tissues and 1443 were up-regulated (Figure 1a). Meanwhile, edgeR identified the 1578 down-regulated DEGs and 3121 up-regulated DEGs (Figure 1b...

A detailed protocol of differential expression analysis methods for RNA sequencing was provided: limma, EdgeR, DESeq2.

RNA sequencing (RNA-seq) is one of the most widely used technologies in transcriptomics as it can reveal the relationship between the genetic alteration ...

three-differential-expression-analysis-methods-for-rna-sequencing

Research

JoVE Journal

Cancer Research

This protocol presents a complete and detailed procedure to apply RNA-seq, a powerful next-generation DNA sequencing technology, to profile transcriptomes in human pulmonary microvascular endothelial cells with or without thrombin treatment. This protocol is generalizable to various cells or tissues affected by different reagents or disease states.

RNA-seq Analysis of Transcriptomes in Thrombin-treated and Control Human Pulmonary Microvascular Endothelial Cells

RNA-sequencing and bioinformatics analyses were used to identify significantly and differentially expressed transcription factors in Lin-CD34+ and Lin-CD34- subpopulations of mouse EMLcells. These transcription factors might play important roles in determining the switch between self-renewing Lin-CD34+ and partially differentiated Lin-CD34- cells.

Identification of Key Factors Regulating Self-renewal and Differentiation in EML Hematopoietic Precursor Cells by RNA-sequencing Analysis

We describe a targeted RNA sequencing-based method that includes preparation of indexed cDNA libraries, hybridization and capture with custom probes and data analysis to interrogate selected transcripts for gene expression, mutations, and gene fusions. Targeted RNAseq permits cost-effective, rapid evaluation of selected transcripts on a desktop sequencer.

Targeted RNA Sequencing Assay to Characterize Gene Expression and Genomic Alterations

This protocol presents an approach for whole transcriptome analysis from zebrafish embryos, larvae, or sorted cells. We include isolation of RNA, pathway analysis of RNASeq data, and qRT-PCR-based validation of gene expression changes.

Sample Preparation and Analysis of RNASeq-based Gene Expression Data from Zebrafish

This protocol describes the combinatorial use of ChIP-seq, 4sU-seq, total RNA-seq, and ribosome profiling for cell lines and primary cells. It enables tracking changes in transcription-factor binding, de novo transcription, RNA processing, turnover and translation over time, and displaying the overall course of events in activated and/or rapidly changing cells.

Real-time Analysis of Transcription Factor Binding, Transcription, Translation, and Turnover to Display Global Events During Cellular Activation

Presented here is a stepwise procedure for in vitro differentiation of human primary keratinocytes by contact inhibition followed by characterization at the molecular level by RNA-seq analysis.

Characterization of In Vitro Differentiation of Human Primary Keratinocytes by RNA-Seq Analysis

This method describes the steps to improve the quality and quantity of sequence data that can be obtained from formalin-fixed paraffin-embedded (FFPE) RNA samples. We describe the methodology to more accurately assess the quality of FFPE-RNA samples, prepare sequencing libraries, and analyze the data from FFPE-RNA samples.

Optimization for Sequencing and Analysis of Degraded FFPE-RNA Samples

DiCoExpress is a script-based tool implemented in R to perform an RNA-Seq analysis from quality control to co-expression. DiCoExpress handles complete and unbalanced design up to 2 biological factors. This video tutorial guides the user through the different features of DiCoExpress.

Analyzing Multifactorial RNA-Seq Experiments with DiCoExpress

Alternative splicing (AS) and alternative polyadenylation (APA) expand the diversity of transcript isoforms and their products. Here, we describe bioinformatic protocols to analyze bulk RNA-seq and 3' end sequencing assays to detect and visualize AS and APA varying across experimental conditions.

Identification of Alternative Splicing and Polyadenylation in RNA-seq Data

The current work describes a protocol for running the Pathway2Targets algorithm, an R script that predicts and prioritizes therapeutic targets based on the profile of intracellular signaling pathways generated by comparing case versus control samples from a bulk RNA-sequencing experiment.

Using Human Differentially Expressed Gene Lists to Perform Downstream Pathway Enrichment Analysis and Target Prioritization

This protocol allows initial quality control for RNA-seq experiments for wet-lab biologists with limited bioinformatics experience.

Rup (RNA-seq Usability Assessment Pipeline) - Quality Control for Bulk RNA-seq Experiments in Eukaryotes

Learn differential expression analysis for RNA sequencing using limma, EdgeR, and DESeq2 in RStudio. This protocol compares methods for identifying differentially expressed genes in Cholangiocarcinoma data.

Three Differential Expression Analysis Methods for RNA Sequencing: limma, EdgeR, DESeq2

Explore More Videos