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 border="1" fo:keep-together.within-page="1" fo:keep-with-next.within-page="always">
	<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

Research

JoVE Journal

Cancer Research

35.2K 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 ...

Two-dimensional Gel Electrophoresis Coupled with Mass Spectrometry Methods for an Analysis of Human Pituitary Adenoma Tissue Proteome

Human pituitary adenoma (PA) is a common tumor that occurs in the human pituitary gland in the hypothalamus-pituitary-targeted organ axis systems, and may be classified as either clinically functional or nonfunctional PA (FPA and NFPA). NFPA is difficult for early stage diagnosis and therapy due to barely elevating hormones in the blood compared to FPA. Our long-term goal is to use proteomics methods to discover reliable biomarkers for clarification of PA molecular mechanisms and recognition of effective diagnostic, prognostic markers and therapeutic targets. Effective two-dimensional gel electrophoresis (2DE) coupled with mass spectrometry (MS) methods were presented here to analyze human PA proteomes, including preparation of samples, 2D gel electrophoresis, protein visualization, image analysis, in-gel trypsin digestion, peptide mass fingerprint (PMF), and tandem mass spectrometry (MS/MS). 2-Dimensional gel electrophoresis matrix-assisted laser desorption/ionization mass spectrometry PMF (2DE-MALDI MS PMF), 2DE-MALDI MS/MS, and 2DE-liquid chromatography (LC) MS/MS procedures have been successfully applied in an analysis of NFPA proteome. With the use of a high-sensitivity mass spectrometer, many proteins were identified with the 2DE-LC-MS/MS method in each 2D gel spot in an analysis of complex PA tissue to maximize the coverage of human PA proteome.

Here, we present a two-dimensional gel electrophoresis (2DE) coupled with mass spectrometry (MS) to separate and identify human pituitary adenoma tissue proteome, which presents a good and reproducible 2DE pattern. Many proteins are observed in each 2DE spot when analyzing complex cancer proteome with the use of high-sensitivity MS.

Human pituitary adenoma (PA) is a common tumor that occurs in the human pituitary gland in the hypothalamus-pituitary-targeted organ axis systems, and may ...

Sequencing Small Non-coding RNA from Formalin-fixed Tissues and Serum-derived Exosomes from Castration-resistant Prostate Cancer Patients

Ablation of androgen receptor (AR) signaling by androgen deprivation is the goal of the first line of therapy for prostate cancer that initially results in cancer regression. However, in a significant number of cases, the disease progresses to advanced, castration-resistant prostate cancer (CRPC), which has limited therapeutic options and is often aggressive. Distant metastasis is mostly observed at this stage of the aggressive disease. CRPC is treated by a second generation of AR pathway inhibitors that improve survival initially, followed by the emergence of therapy resistance. Neuroendocrine prostate cancer (NEPC) is a rare variant of prostate cancer (PCa) that often develops as a result of therapy resistance via a transdifferentiation process known as neuroendocrine differentiation (NED), wherein PCa cells undergo a lineage switch from adenocarcinomas and show increased expression of neuroendocrine (NE) lineage markers. In addition to the genomic alterations that drive the progression and transdifferentiation to NEPC, epigenetic factors and microenvironmental cues are considered essential players in driving disease progression. This manuscript provides a detailed protocol to identify the epigenetic drivers (i.e., small non-coding RNAs) that are associated with advanced PCa. Using purified microRNAs from formalin-fixed paraffin-embedded (FFPE) metastatic tissues and corresponding serum-derived extracellular vesicles (EVs), the protocol describes how to prepare libraries with appropriate quality control for sequencing microRNAs from these sample sources. Isolating RNA from both FFPE and EVs is often challenging because most of it is either degraded or is limited in quantity. This protocol will elaborate on different methods to optimize the RNA inputs and cDNA libraries to yield most specific reads and high-quality data upon sequencing.

Therapy resistance often develops in patients with advanced prostate cancer, and in some cases, cancer progresses to a lethal subtype called neuroendocrine prostate cancer. Assessing the small non-coding RNA-mediated molecular changes that facilitate this transition would allow better disease stratification and identification of causal mechanisms that lead to development of neuroendocrine prostate cancer.

Ablation of androgen receptor (AR) signaling by androgen deprivation is the goal of the first line of therapy for prostate cancer that initially results ...

Evaluating Cell Death Using Cell-Free Supernatant of Probiotics in Three-Dimensional Spheroid Cultures of Colorectal Cancer Cells

This manuscript describes a protocol to evaluate cancer cell deaths in three dimensional (3D) spheroids of multicellular types of cancer cells using supernatants from Lactobacillus fermentum cell culture, considered as probiotics cultures. The use of 3D cultures to test Lactobacillus cell-free supernatant (LCFS) are a better option than testing in 2D monolayers, especially as L. fermentum can produce anti-cancer effects within the gut. L. fermentum supernatant was identified to possess increased anti-proliferative effects against several colorectal cancer (CRC) cells in 3D culture conditions. Interestingly, these effects were strongly related to the culture model, demonstrating the notable ability of L. fermentum to induce cancer cell death. Stable spheroids were generated from diverse CRCs (colorectal cancer cells) using the protocol presented below. This protocol of generating 3D spheroid is time saving and cost effective. This system was developed to easily investigate the anti-cancer effects of LCFS in multiple types of CRC spheroids. As expected, CRC spheroids treated with LCFS strongly induced cell death during the experiment and expressed specific apoptosis molecular markers as analyzed by qRT-PCR, western blotting, and FACS analysis. Therefore, this method is valuable for exploring cell viability and evaluating the efficacy of anti-cancer drugs.

Here methods are presented to understand anti-cancer effects of Lactobacillus cell-free supernatant (LCFS). Colorectal cancer cell lines show cell deaths when treated with LCFS in 3D cultures. The process of generating spheroids can be optimized depending on the scaffold and the analysis methods presented are useful for evaluating the involved signaling pathways.

This manuscript describes a protocol to evaluate cancer cell deaths in three dimensional (3D) spheroids of multicellular types of cancer cells using ...

An Automated Differential Nuclear Staining Assay for Accurate Determination of Mitocan Cytotoxicity

The contribution of mitochondria to oncogenic transformation is a subject of wide interest and active study. As the field of cancer metabolism becomes more complex, the goal of targeting mitochondria using various compounds that inflict mitochondrial damage (so-called mitocans) is becoming quite popular. Unfortunately, many existing cytotoxicity assays, such as those based on tetrazolium salts or resazurin require functional mitochondrial enzymes for their performance. The damage inflicted by compounds that target mitochondria often compromises the accuracy of these assays. Here, we describe a modified protocol based on differential staining with two fluorescent dyes, one of which is cell-permeant (Hoechst 33342) and the other of which is not (propidium iodide). The difference in staining allows living and dead cells to be discriminated. The assay is amenable to automated microscopy and image analysis, which increases throughput and reduces bias. This also allows the assay to be used in high-throughput fashion using 96-well plates, making it a viable option for drug discovery efforts, particularly when the drugs in question have some level of mitotoxicity. Importantly, results obtained by Hoechst/PI staining assay show increased consistency, both with trypan blue exclusion results and between biological replicates when the assay is compared to other methods.

The protocol describes a rapid, high-throughput, reliable, inexpensive, and unbiased assay for efficiently determining cellular viability. This assay is particularly useful when cells&#39; mitochondria have been damaged, which interferes with other assays. The assay uses automated counting of cells stained with two nuclear dyes &#8211; Hoechst 33342 and propidium iodide.

The contribution of mitochondria to oncogenic transformation is a subject of wide interest and active study. As the field of cancer metabolism becomes ...

Spatial Profiling of Protein and RNA Expression in Tissue: An Approach to Fine-Tune Virtual Microdissection

Multiplexing enables the assessment of several markers on the same tissue while providing spatial context. Spatial Omics technologies allow both protein and RNA multiplexing by leveraging photo-cleavable oligo-tagged antibodies and probes, respectively. Oligos are cleaved and quantified from specific regions across the tissue to elucidate the underlying biology. Here, the study demonstrates that automated custom antibody visualization protocols can be utilized to guide ROI selection in conjunction with spatial proteomics assays. This specific method did not show acceptable performance with spatial transcriptomics assays. The protocol describes the development of a 3-plex immunofluorescent (IF) assay for marker visualization on an automated platform, using tyramide signal amplification (TSA) to amplify the fluorescent signal from a given protein target and increase the antibody pool to choose from. The visualization protocol was automated using a thoroughly validated 3-plex assay to ensure quality and reproducibility. In addition, the exchange of DAPI for SYTO dyes was evaluated to allow imaging of TSA-based IF assays on the spatial profiling platform. Additionally, we tested the ability of selecting small ROIs using the spatial transcriptomics assay to allow the investigation of highly-specific areas of interest (e.g., areas enriched for a given cell type). ROIs of 50 &#181;m and 300 &#181;m diameter were collected, which corresponds to approximately 15 cells and 100 cells, respectively. Samples were made into libraries and sequenced to investigate the capability to detect signals from small ROIs and profile-specific regions of the tissue. We determined that spatial proteomics technologies highly benefit from automated, standardized protocols to guide ROI selection. While this automated visualization protocol was not compatible with spatial transcriptomics assays, we were able to test and confirm that specific cell populations can successfully be detected even in small ROIs with the standard manual visualization protocol.

Here, we describe a protocol for fine-tuning regions of interest (ROIs) for Spatial Omics technologies to better characterize the tumor microenvironment and identify specific cell populations. For proteomics assays, automated customized protocols can guide ROI selection, while transcriptomics assays can be fine-tuned utilizing ROIs as small as 50 &#181;m.

Multiplexing enables the assessment of several markers on the same tissue while providing spatial context. Spatial Omics technologies allow both protein ...

Transmitochondrial Cybrid Generation From Suspension-Growing Cancer Cells

Transmitochondrial Cybrid Generation Using Cancer Cell Lines

In recent years, the number of studies dedicated to ascertaining the connection between mitochondria and cancer has significantly risen. However, more efforts are still needed to fully understand the link involving alterations in mitochondria and tumorigenesis, as well as to identify tumor-associated mitochondrial phenotypes. For instance, to evaluate the contribution of mitochondria in tumorigenesis and metastasis processes, it is essential to understand the influence of mitochondria from tumor cells in different nuclear environments. For this purpose, one possible approach consists of transferring mitochondria into a different nuclear background to obtain the so-called cybrid cells. In the traditional cybridization techniques, a cell line lacking mtDNA (&#961;0, nuclear donor cell) is repopulated with mitochondria derived from either enucleated cells or platelets. However, the enucleation process requires good cell adhesion to the culture plate, a feature that is partially or completely lost in many cases in invasive cells. In addition, another difficulty found in the traditional methods is achieving complete removal of the endogenous mtDNA from the mitochondrial-recipient cell line to obtain pure nuclear and mitochondrial DNA backgrounds, avoiding the presence of two different mtDNA species in the generated cybrid. In this work, we present a mitochondrial exchange protocol applied to suspension-growing cancer cells based on the repopulation of rhodamine 6G-pretreated cells with isolated mitochondria. This methodology allows us to overcome the limitations of the traditional approaches, and thus can be used as a tool to expand the comprehension of the mitochondrial role in cancer progression and metastasis.

Author Spotlight: Transmitochondrial Cybrid Generation Using Cancer Cell Lines  

This protocol describes a technique for cybrid generation from suspension-growing cancer cells as a tool to study the role of mitochondria in the tumorigenic process.

In recent years, the number of studies dedicated to ascertaining the connection between mitochondria and cancer has significantly risen. However, more ...

Isolation of Nuclei for Multiome Sequencing in Solid Tumor Specimen

Optimized Nuclei Isolation from Fresh and Frozen Solid Tumor Specimens for Multiome Sequencing

Multiome sequencing, which provides same-cell/paired single-cell RNA- and the assay for transposase-accessible chromatin with sequencing (ATAC-sequencing) data, represents a breakthrough in our ability to discern tumor cell heterogeneity-a primary focus of translational cancer research at this time. However, the quality of sequencing data acquired using this advanced modality is highly dependent on the quality of the input material. 
Digestion conditions need to be optimized to maximize cell yield without sacrificing quality. This is particularly challenging in the context of solid tumors with dense desmoplastic matrices that must be gently broken down for cell release. Freshly isolated cells from solid tumor tissue are more fragile than those isolated from cell lines. Additionally, as the cell types isolated are heterogeneous, conditions should be selected to support the total cell population. 
Finally, nuclear isolation conditions must be optimized based on these qualities in terms of lysis times and reagent types/ratios. In this article, we describe our experience with nuclear isolation for the 10x Genomics multiome sequencing platform from solid tumor specimens. We provide recommendations for tissue digestion, storage of single-cell suspensions (if desired), and nuclear isolation and assessment.

Author Spotlight: Enhancing Nuclei Isolation for Multiome Sequencing in Challenging Tumor Microenvironments

The protocol provides a reliable and optimized approach to the isolation of nuclei from solid tumor specimens for multiome sequencing using the 10x Genomics platform, including recommendations for tissue dissociation conditions, cryopreservation of single-cell suspensions, and assessment of isolated nuclei.

Multiome sequencing, which provides same-cell/paired single-cell RNA- and the assay for transposase-accessible chromatin with sequencing (ATAC-sequencing) ...

Diagnosis of Vitreoretinal Lymphoma Using Cell-Free DNA

Cell-Free DNA Extraction of Vitreous and Aqueous Humor Specimens for Diagnosis and Monitoring of Vitreoretinal Lymphoma

Vitreoretinal lymphoma (VRL) represents an aggressive lymphoma, often categorized as primary central nervous system diffuse large B-cell lymphoma. To diagnose VRL, specimens such as vitreous humor and, more recently, aqueous humor are collected. Diagnostic testing for VRL on these specimens includes cytology, flow cytometry, and molecular testing. However, both cytopathology and flow cytometry, along with molecular testing using cellular DNA, necessitate intact whole cells. The challenge lies in the fact that vitreous and aqueous humor typically have low cellularity, and many cells get destroyed during collection, storage, and processing. Moreover, these specimens pose additional difficulties for molecular testing due to the high viscosity of vitreous humor and the low volume of both vitreous and aqueous humor. This study proposes a method for extracting cell-free DNA from vitreous and aqueous specimens. This approach complements the extraction of cellular DNA or allows the cellular component of these specimens to be utilized for other diagnostic methods, including cytology and flow cytometry.

Author Spotlight: Advancing VRL Diagnosis Using Cell-Free DNA Extraction from Vitreous Humor 

A procedure to extract cell-free DNA from vitreous and aqueous humor to perform molecular studies for diagnosing vitreoretinal lymphoma is established here. The method offers the ability to concurrently extract DNA from the cellular component of the sample or to reserve it for ancillary testing.

Vitreoretinal lymphoma (VRL) represents an aggressive lymphoma, often categorized as primary central nervous system diffuse large B-cell lymphoma. To ...