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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Differential+Expression+Analysis+in+RNA-seq+Data+Using+a+Geometric+Approach.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=RNA-Seq:+a+revolutionary+tool+for+transcriptomics.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Count-based+differential+expression+analysis+of+RNA+sequencing+data+using+R+and+Bioconductor.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Interpretation+of+differential+gene+expression+results+of+RNA-seq+data:+review+and+integration.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=RNA-Seq+differential+expression+analysis:+An+extended+review+and+a+software+tool.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=RNA-seq+analysis+is+easy+as+1-2-3+with+limma,+Glimma+and+edgeR.">RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR.</a> F1000Research. 5, (2016).</li><li>Varet, H., Brillet-Gu&#233;guen, L., Copp&#233;e, J. Y., Dillies, M. A. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=SARTools:+A+DESeq2-+and+EdgeR-Based+R+Pipeline+for+Comprehensive+Differential+Analysis+of+RNA-Seq+Data.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=limma+powers+differential+expression+analyses+for+RNA-sequencing+and+microarray+studies.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=edgeR:+a+Bioconductor+package+for+differential+expression+analysis+of+digital+gene+expression+data.">edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.</a> Bioinformatics. 26, 139-140 (2010).</li><li>Love, M. I., Huber, W., Anders, S. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Moderated+estimation+of+fold+change+and+dispersion+for+RNA-seq+data+with+DESeq2.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Bioconductor:+open+software+development+for+computational+biology+and+bioinformatics.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=voom:+Precision+weights+unlock+linear+model+analysis+tools+for+RNA-seq+read+counts.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Linear+models+and+empirical+bayes+methods+for+assessing+differential+expression+in+microarray+experiments.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Detecting+differential+expression+in+RNA-sequence+data+using+quasi-likelihood+with+shrunken+dispersion+estimates.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Evaluating+statistical+analysis+models+for+RNA+sequencing+experiments.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Excess+False+Positive+Rates+in+Methods+for+Differential+Gene+Expression+Analysis+using+RNA-Seq+Data.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=Comparison+and+calibration+of+transcriptome+data+from+RNA-Seq+and+tiling+arrays.">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/entrez/query.fcgi?db=PubMed&cmd=Search&doptcmdl=Citation&defaultField=Title+Word&term=EBSeq:+an+empirical+Bayes+hierarchical+model+for+inference+in+RNA-seq+experiments.">EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments.</a> Bioinformatics. 29, 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

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|&#8805;2 and adj.P.Val &#60;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</...

Watch this Scientific Journal Video about Three Differential Expression Analysis Methods for RNA Sequencing: limma, EdgeR, DESeq2 at JoVE.com

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

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

Universidad San Sebastian

Research

JoVE Journal

Cancer Research

36.4K Views.  Renmin Hospital of Wuhan University. A detailed protocol of differential expression analysis methods for RNA sequencing was provided: limma, EdgeR, DESeq2.

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

Using RNA-sequencing to Detect Novel Splice Variants Related to Drug Resistance in In Vitro Cancer Models

Drug resistance remains a major problem in the treatment of cancer for both hematological malignancies and solid tumors. Intrinsic or acquired resistance can be caused by a range of mechanisms, including increased drug elimination, decreased drug uptake, drug inactivation and alterations of drug targets. Recent data showed that other than by well-known genetic (mutation, amplification) and epigenetic (DNA hypermethylation, histone post-translational modification) modifications, drug resistance mechanisms might also be regulated by splicing aberrations. This is a rapidly growing field of investigation that deserves future attention in order to plan more effective therapeutic approaches. The protocol described in this paper is aimed at investigating the impact of aberrant splicing on drug resistance in solid tumors and hematological malignancies. To this goal, we analyzed the transcriptomic profiles of several in vitro models through RNA-seq and established a qRT-PCR based method to validate candidate genes. In particular, we evaluated the differential splicing of DDX5 and PKM transcripts. The aberrant splicing detected by the computational tool MATS was validated in leukemic cells, showing that different DDX5 splice variants are expressed in the parental vs. resistant cells. In these cells, we also observed a higher PKM2/PKM1 ratio, which was not detected in the Panc-1 gemcitabine-resistant counterpart compared to parental Panc-1 cells, suggesting a different mechanism of drug-resistance induced by gemcitabine exposure.

Using RNA-sequencing to Detect Novel Splice Variants Related to Drug Resistance in In Vitro Cancer Models

Here we describe a protocol aimed at investigating the impact of aberrant splicing on drug resistance in solid tumors and hematological malignancies. To this goal, we analyzed the transcriptomic profiles of parental and resistant in vitro models through RNA-seq and established a qRT-PCR based method to validate candidate genes.

Drug resistance remains a major problem in the treatment of cancer for both hematological malignancies and solid tumors. Intrinsic or acquired resistance ...

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

The analysis of global gene expression changes is a valuable tool for identifying novel pathways underlying observed phenotypes. The zebrafish is an excellent model for rapid assessment of whole transcriptome from whole animal or individual cell populations due to the ease of isolation of RNA from large numbers of animals. Here a protocol for global gene expression analysis in zebrafish embryos using RNA sequencing (RNASeq) is presented. We describe preparation of RNA from whole embryos or from cell populations obtained using cell sorting in transgenic animals. We also describe an approach for analysis of RNASeq data to identify enriched pathways and Gene Ontology (GO) terms in global gene expression data sets. Finally, we provide a protocol for validation of gene expression changes using quantitative reverse transcriptase PCR (qRT-PCR). These protocols can be used for comparative analysis of control and experimental sets of zebrafish to identify novel gene expression changes, and provide molecular insight into phenotypes of interest.

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.

The analysis of global gene expression changes is a valuable tool for identifying novel pathways underlying observed phenotypes. The zebrafish is an ...

Genetics

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

Upon activation, cells rapidly change their functional programs and, thereby, their gene expression profile. Massive changes in gene expression occur, for example, during cellular differentiation, morphogenesis, and functional stimulation (such as activation of immune cells), or after exposure to drugs and other factors from the local environment. Depending on the stimulus and cell type, these changes occur rapidly and at any possible level of gene regulation. Displaying all molecular processes of a responding cell to a certain type of stimulus/drug is one of the hardest tasks in molecular biology. Here, we describe a protocol that enables the simultaneous analysis of multiple layers of gene regulation. We compare, in particular, transcription factor binding (Chromatin-immunoprecipitation-sequencing (ChIP-seq)), de novo transcription (4-thiouridine-sequencing (4sU-seq)), mRNA processing, and turnover as well as translation (ribosome profiling). By combining these methods, it is possible to display a detailed and genome-wide course of action.
Sequencing newly transcribed RNA is especially recommended when analyzing rapidly adapting or changing systems, since this depicts the transcriptional activity of all genes during the time of 4sU exposure (irrespective of whether they are up- or downregulated). The combinatorial use of total RNA-seq and ribosome profiling additionally allows the calculation of RNA turnover and translation rates. Bioinformatic analysis of high-throughput sequencing results allows for many means for analysis and interpretation of the data. The generated data also enables tracking co-transcriptional and alternative splicing, just to mention a few possible outcomes.
The combined approach described here can be applied for different model organisms or cell types, including primary cells. Furthermore, we provide detailed protocols for each method used, including quality controls, and discuss potential problems and pitfalls.

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.

Upon activation, cells rapidly change their functional programs and, thereby, their gene expression profile. Massive changes in gene expression occur, for ...

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

Human primary keratinocytes are often used as in vitro models for studies on epidermal differentiation and related diseases. Methods have been reported for in vitro differentiation of keratinocytes cultured in two-dimensional (2D) submerged manners using various induction conditions. Described here is a procedure for 2D in vitro keratinocyte differentiation method by contact inhibition and subsequent molecular characterization by RNA-seq. In brief, keratinocytes are grown in defined keratinocyte medium supplemented with growth factors until they are fully confluent. Differentiation is induced by close contacts between the keratinocytes and further stimulated by excluding growth factors in the medium. Using RNA-seq analyses, it is shown that both 1) differentiated keratinocytes exhibit distinct molecular signatures during differentiation and 2) the dynamic gene expression pattern largely resembles cells during epidermal stratification. As for comparison to normal keratinocyte differentiation, keratinocytes carrying mutations of the transcription factor p63 exhibit altered morphology and molecular signatures, consistent with their differentiation defects. In conclusion, this protocol details the steps for 2D in vitro keratinocyte differentiation and its molecular characterization, with an emphasis on bioinformatic analysis of RNA-seq data. Because RNA extraction and RNA-seq procedures have been well-documented, it is not the focus of this protocol. The experimental procedure of in vitro keratinocyte differentiation and bioinformatic analysis pipeline can be used to study molecular events during epidermal differentiation in healthy and diseased keratinocytes.

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.

Human primary keratinocytes are often used as in vitro models for studies on epidermal differentiation and related diseases. Methods have been reported ...

Optimization for Sequencing and Analysis of Degraded FFPE-RNA Samples

Gene expression analysis by RNA sequencing (RNA-seq) enables unique insights into clinical samples that can potentially lead to mechanistic understanding of the basis of various diseases as well as resistance and/or susceptibility mechanisms. However, FFPE tissues, which represent the most common method for preserving tissue morphology in clinical specimens, are not the best sources for gene expression profiling analysis. The RNA obtained from such samples is often degraded, fragmented, and chemically modified, which leads to suboptimal sequencing libraries. In turn, these generate poor quality sequence data that may not be reliable for gene expression analysis and mutation discovery. In order to make the most of FFPE samples and obtain the best possible data from low quality samples, it is important to take certain precautions while planning experimental design, preparing sequencing libraries, and during data analysis. This includes the use of appropriate metrics for precise sample quality control (QC), identifying the best methods for various steps during the sequencing library generation, and careful library QC. In addition, applying correct software tools and parameters for sequence data analysis is critical in order to identify artifacts in RNA-seq data, filter out contamination and low quality reads, assess uniformity of gene coverage, and measure the reproducibility of gene expression profiles among biological replicates. These steps can ensure high accuracy and reproducibility for profiling of very heterogeneous RNA samples. Here we describe the various steps for sample QC, library preparation and QC, sequencing, and data analysis that can help to increase the amount of useful data obtained from low quality RNA, such as that obtained from FFPE-RNA tissues.

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.

Gene expression analysis by RNA sequencing (RNA-seq) enables unique insights into clinical samples that can potentially lead to mechanistic understanding ...

Analyzing Multifactorial RNA-Seq Experiments with DiCoExpress

The proper use of statistical modeling in NGS data analysis requires an advanced level of expertise. There has recently been a growing consensus on using generalized linear models for differential analysis of RNA-Seq data and the advantage of mixture models to perform co-expression analysis. To offer a managed setting to use these modeling approaches, we developed DiCoExpress that provides a standardized R pipeline to perform an RNA-Seq analysis. Without any particular knowledge in statistics or R programming, beginners can perform a complete RNA-Seq analysis from quality controls to co-expression through differential analysis based on contrasts inside a generalized linear model. An enrichment analysis is proposed both on the lists of differentially expressed genes, and the co-expressed gene clusters. This video tutorial is conceived as a step-by-step protocol to help users take full advantage of DiCoExpress and its potential in empowering the biological interpretation of an RNA-Seq experiment.

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.

The proper use of statistical modeling in NGS data analysis requires an advanced level of expertise. There has recently been a growing consensus on using ...

Engineering

RNA-Seq Analysis of Differential Gene Expression in Electroporated Chick Embryonic Spinal Cord

In ovo electroporation of the chick neural tube is a fast and inexpensive method for identification of gene function during neural development. Genome wide analysis of differentially expressed transcripts after such an experimental manipulation has the potential to uncover an almost complete picture of the downstream effects caused by the transfected construct. This work describes a simple method for comparing transcriptomes from samples of transfected embryonic spinal cords comprising all steps between electroporation and identification of differentially expressed transcripts. The first stage consists of guidelines for electroporation and instructions for dissection of transfected spinal cord halves from HH23 embryos in ribonuclease-free environment and extraction of high-quality RNA samples suitable for transcriptome sequencing. The next stage is that of bioinformatic analysis with general guidelines for filtering and comparison of RNA-Seq datasets in the Galaxy public server, which eliminates the need of a local computational structure for small to medium scale experiments. The representative results show that the dissection methods generate high quality RNA samples and that the transcriptomes obtained from two control samples are essentially the same, an important requirement for detection of differential expression genes in experimental samples. Furthermore, one example is provided where experimental overexpression of a DNA construct can be visually verified after comparison with control samples. The application of this method may be a powerful tool to facilitate new discoveries on the function of neural factors involved in spinal cord early development.

This video shows the basic steps for performing whole transcriptome analysis on dissected chick embryonic spinal cord samples after transfection with in ovo electroporation.

In ovo electroporation of the chick neural tube is a fast and inexpensive method for identification of gene function during neural development. Genome ...

Neuroscience

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

Hematopoietic stem cells (HSCs) are used clinically for transplantation treatment to rebuild a patient&#39;s hematopoietic system in many diseases such as leukemia and lymphoma. Elucidating the mechanisms controlling HSCs self-renewal and differentiation is important for application of HSCs for research and clinical uses. However, it is not possible to obtain large quantity of HSCs due to their inability to proliferate in vitro. To overcome this hurdle, we used a mouse bone marrow derived cell line, the EML (Erythroid, Myeloid, and Lymphocytic) cell line, as a model system for this study.
RNA-sequencing (RNA-Seq) has been increasingly used to replace microarray for gene expression studies. We report here a detailed method of using RNA-Seq technology to investigate the potential key factors in regulation of EML cell self-renewal and differentiation. The protocol provided in this paper is divided into three parts. The first part explains how to culture EML cells and separate Lin-CD34+ and Lin-CD34- cells. The second part of the protocol offers detailed procedures for total RNA preparation and the subsequent library construction for high-throughput sequencing. The last part describes the method for RNA-Seq data analysis and explains how to use the data to identify differentially expressed transcription factors between Lin-CD34+ and Lin-CD34- cells. The most significantly differentially expressed transcription factors were identified to be the potential key regulators controlling EML cell self-renewal and differentiation. In the discussion section of this paper, we highlight the key steps for successful performance of this experiment.
In summary, this paper offers a method of using RNA-Seq technology to identify potential regulators of self-renewal and differentiation in EML cells. The key factors identified are subjected to downstream functional analysis in vitro and in vivo.

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.

Hematopoietic stem cells (HSCs) are used clinically for transplantation treatment to rebuild a patient&#39;s hematopoietic system in many diseases such as ...

Biology

Targeted RNA Sequencing Assay to Characterize Gene Expression and Genomic Alterations

RNA sequencing (RNAseq) is a versatile method that can be utilized to detect and characterize gene expression, mutations, gene fusions, and noncoding RNAs. Standard RNAseq requires 30 - 100 million sequencing reads and can include multiple RNA products such as mRNA and noncoding RNAs. We demonstrate how targeted RNAseq (capture) permits a focused study on selected RNA products using a desktop sequencer. RNAseq capture can characterize unannotated, low, or transiently expressed transcripts that may otherwise be missed using traditional RNAseq methods. Here we describe the extraction of RNA from cell lines, ribosomal RNA depletion, cDNA synthesis, preparation of barcoded libraries, hybridization and capture of targeted transcripts and multiplex sequencing on a desktop sequencer. We also outline the computational analysis pipeline, which includes quality control assessment, alignment, fusion detection, gene expression quantification and identification of single nucleotide variants. This assay allows for targeted transcript sequencing to characterize gene expression, gene fusions, and mutations.

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.

RNA sequencing (RNAseq) is a versatile method that can be utilized to detect and characterize gene expression, mutations, gene fusions, and noncoding ...

Identification of Alternative Splicing and Polyadenylation in RNA-seq Data

As well as the typical analysis of RNA-Seq to measure differential gene expression (DGE) across experimental/biological conditions, RNA-seq data can also be utilized to explore other complex regulatory mechanisms at the exon level. Alternative splicing and polyadenylation play a crucial role in the functional diversity of a gene by generating different isoforms to regulate gene expression at the post-transcriptional level, and limiting analyses to the whole gene level can miss this important regulatory layer. Here, we demonstrate detailed step by step analyses for identification and visualization of differential exon and polyadenylation site usage across conditions, using Bioconductor and other packages and functions, including DEXSeq, diffSplice from the Limma package, and rMATS.

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.

As well as the typical analysis of RNA-Seq to measure differential gene expression (DGE) across experimental/biological conditions, RNA-seq data can also ...

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

Summary

Explore More Videos

Chapters in this video

Using RNA-sequencing to Detect Novel Splice Variants Related to Drug Resistance in In Vitro Cancer Models

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

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

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

Optimization for Sequencing and Analysis of Degraded FFPE-RNA Samples

Analyzing Multifactorial RNA-Seq Experiments with DiCoExpress

RNA-Seq Analysis of Differential Gene Expression in Electroporated Chick Embryonic Spinal Cord

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

Targeted RNA Sequencing Assay to Characterize Gene Expression and Genomic Alterations

Identification of Alternative Splicing and Polyadenylation in RNA-seq Data