	DiBiG
ICBR Bioinformatics	Powered by Actor, v1.0

RNAseq - Alignment and differential expression analysis

Title: GE7218
Project: (none)
Started on: 1/5/2024 11:56:54
Hostname: login7.ufhpc
Run directory: /blue/licht/runs/H2B-E77K-Project/GE7218/GE7218
Configuration GE7218.conf

Table of contents:

Input data
Trimming and quality control
Alignment to transcriptome
Genome coverage
Expression analysis - quantification
Differential expression - protein-coding genes
Differential expression - all genes
Differential expression - isoform level
Differential expression - combined files
Alternative splicing analysis
MultiQC report
UCSC hub
Methods summary

1. Input data
The following table summarizes the samples, conditions, and contrasts in this analysis. A readset is either a single fastq file or a pair of fastq files (for paired-end sequencing).

Category	Data
*Summary of input data*
Reference genome:	hg38
Experimental conditions:	WT, HET, MUT
Contrasts:	HET vs. WT, MUT vs. WT
Number of samples	9
*Sequencing data data*
Total number of reads:	403,987,071
Average reads per sample:	44,887,452

Table 1. Summary of input data

Condition	Sample	Number of reads	% Reads
WT	WT1	42,810,999	10.60%
	WT2	42,591,089	10.54%
	WT3	48,038,535	11.89%
HET	HET1	42,992,044	10.64%
	HET2	51,532,469	12.76%
	HET3	40,280,567	9.97%
MUT	MUT1	41,032,996	10.16%
	MUT2	47,955,844	11.87%
	MUT3	46,752,528	11.57%

Table 2. Number of reads in each sample.

2. Trimming and quality control
The input sequences were trimmed using fastp (version 0.22.0). The following table provides links to the quality control reports after trimming, as well as the number of reads in the trimmed files.

Sample	Readset	Reads before trim	Reads after trim	QC after trim	% Retained
WT1	WT1_r1	42,810,999	35,172,826	WT1_r1	82.16%
WT2	WT2_r1	42,591,089	31,307,401	WT2_r1	73.51%
WT3	WT3_r1	48,038,535	39,831,210	WT3_r1	82.92%
HET1	HET1_r1	42,992,044	35,713,353	HET1_r1	83.07%
HET2	HET2_r1	51,532,469	43,027,255	HET2_r1	83.50%
HET3	HET3_r1	40,280,567	32,404,304	HET3_r1	80.45%
MUT1	MUT1_r1	41,032,996	31,882,547	MUT1_r1	77.70%
MUT2	MUT2_r1	47,955,844	37,180,855	MUT2_r1	77.53%
MUT3	MUT3_r1	46,752,528	38,696,989	MUT3_r1	82.77%

Table 3. Number of reads in input files and links to QC reports.

3. Alignment to transcriptome
The input sequences were aligned to the hg38 transcriptome using 2.7.9a. The following table reports the number of alignments to the genome and the transcriptome for each sample. Please note that the number of alignments will in general be higher than the number of reads because the same read may align to multiple isoforms of the same gene. The WIG files can be uploaded to the UCSC Genome Browser as custom tracks.

Sample	Input reads	Genome alignments	Genome alignment rate	Transcriptome alignments	Transcriptome alignment rate	Alignment report
WT1	35,172,826	98,379,023	2.80	27,954,386	79.48%	WT1.star/Log.final.out
WT2	31,307,401	96,938,332	3.10	21,397,769	68.35%	WT2.star/Log.final.out
WT3	39,831,210	89,852,146	2.26	36,210,567	90.91%	WT3.star/Log.final.out
HET1	35,713,353	82,528,094	2.31	32,140,232	89.99%	HET1.star/Log.final.out
HET2	43,027,255	98,725,795	2.29	39,009,598	90.66%	HET2.star/Log.final.out
HET3	32,404,304	78,901,854	2.43	27,305,815	84.27%	HET3.star/Log.final.out
MUT1	31,882,547	101,174,535	3.17	21,636,097	67.86%	MUT1.star/Log.final.out
MUT2	37,180,855	125,750,994	3.38	18,698,944	50.29%	MUT2.star/Log.final.out
MUT3	38,696,989	106,340,044	2.75	30,806,465	79.61%	MUT3.star/Log.final.out

Table 4. Number of alignments to genome and transcriptome.

4. Genome coverage
The following table reports the overall and effective genome coverage in each sample. The Total nt column reports the total number of nucleotides sequenced, i.e. the number of aligned reads times the length of each read. Coverage is this number divided by the size of the genome. Effective bp reports the number of bases in the genome having coverage greater than 5, and the Effective Perc column shows what percentage this is of the genome size. Note that, especially in the case of RNA-seq, the effective genome size may be much smaller than the full size. Eff Coverage is the average coverage over the effectively covered fraction of the genome.

Name	Total nt	Coverage	Effective bp	Effective Perc	Eff Coverage
WT1	91,572,436,618	29.97	684,762,924	22.40%	133.73
WT2	47,631,381,297	15.43	622,793,317	20.20%	76.48
WT3	120,049,534,001	38.90	715,037,579	23.20%	167.89
HET1	110,121,431,347	35.69	710,918,761	23.00%	154.90
HET2	125,725,940,172	40.75	728,198,692	23.60%	172.65
HET3	80,867,465,551	26.20	685,053,628	22.20%	118.05
MUT1	57,467,007,347	18.63	565,599,303	18.30%	101.60
MUT2	45,635,720,548	14.79	604,648,569	19.60%	75.47
MUT3	91,206,744,992	29.55	743,906,817	24.10%	122.61

Table 5. Genome coverage by sample.

The following table reports the overall and effective genome coverage in each condition.

Name	Total nt	Coverage	Effective bp	Effective Perc	Eff Coverage
WT	250,709,316,560	81.22	828,833,093	26.90%	302.48
HET	305,756,307,569	99.06	872,726,461	28.30%	350.35
MUT	186,744,183,109	60.50	921,187,673	29.80%	202.72

Table 6. Genome coverage by condition

File: GE7218.sample.cov.xlsx
Size: 25.47 kB
Description: Per-chromosome coverage data, by sample.

File: GE7218.cond.cov.xlsx
Size: 11.80 kB
Description: Per-chromosome coverage data, by condition.

5. Expression analysis - quantification
Gene and transcript expression values were quantified using RSEM v1.3.1. The following files contain the raw FPKM values for all genes/transcripts in all samples. NOTE: these values are not normalized yet, please apply the appropriate normalization before using them in analysis.

File: genes.rawmatrix.csv
Size: 4.52 MB
Description: Matrix of FPKM values for all genes in all samples.

File: transcripts.rawmatrix.csv
Size: 16.73 MB
Description: Matrix of FPKM values for all transcripts in all samples.

File: genes.xpra.txt
Size: 3.41 MB
Description: Counts table suitable for ExpressAnalyst.

The following scatterplots show the level of similarity between replicates of the same condition.

Principal Component Analysis on raw (un-normalized) expression data. Click on the thumbnail to display the full-size image.

(png format, 82.32 kB)

The following image displays the Multi-Dimensional Scaling (MDS) plot for the raw (un-normalized) expression data. Click on the thumbnail to display the full-size image.

6. Differential expression - protein-coding genes
Differential gene expression was analyzed using DESeq2. The following table reports the number of differentially expressed genes in each contrast with abs(log2(FC)) >= 1.0 and FDR-corrected P-value <= 0.05. The files under the Table heading contain the log2(FC) and P-value of all significant genes, while the files under the Expressions heading contain normalized expression values for the significant genes in all replicates of the two conditions being compared. The lists of differentially expressed genes for all contrasts can also be downloaded as a single Excel file using the link below.

Test	Control	Total	Overexpressed	Underexpressed	Table	Expressions
HET	WT	130	52	78	HET.vs.WT.codinggeneDiff.csv	HET.vs.WT.gmatrix.csv
MUT	WT	1,503	868	635	MUT.vs.WT.codinggeneDiff.csv	MUT.vs.WT.gmatrix.csv

Table 7. Results of gene-level differential expression analysis.

File: GE7218-codingdiff.xlsx
Size: 130.15 kB
Description: Excel file containing differentially expressed genes for all contrasts (one sheet per contrast). Only includes protein-coding genes.

File: GE7218-allcodingdiff.xlsx
Size: 1.75 MB
Description: Excel file containing differential expression values for all tested genes in all contrasts (one sheet per contrast). Only includes protein-coding genes. Note: genes with very low average expression in all conditions were removed.

File: GE7218.g.deseq2norm.xlsx
Size: 2.18 MB
Description: Excel file containing normalized (DESeq2) expression values for all protein-coding genes in all conditions. Note: genes with very low average expression in all conditions were removed.

Principal Component Analysis on normalized expression data. Click on the thumbnail to display the full-size image.

(png format, 91.08 kB)

The following image displays the Multi-Dimensional Scaling (MDS) plot for the normalized expression data. In this plot, relative distances between samples reflect the similarity of their gene expression profiles. Ideally, replicates of the same condition should be close together, and well separated from other conditions.

Volcano plots for all contrasts. Use the menu to select a contrast.

7. Differential expression - all genes
The following table reports results from the same differential analysis as above, but includes all biotypes instead of coding genes only.

Test	Control	Total	Overexpressed	Underexpressed	Table	Expressions
HET	WT	136	54	82	HET.vs.WT.geneDiff.csv	HET.vs.WT.gmatrix.csv
MUT	WT	1,687	1,025	662	MUT.vs.WT.geneDiff.csv	MUT.vs.WT.gmatrix.csv

Table 8. Results of gene-level differential expression analysis (all biotypes).

File: GE7218-genediff.xlsx
Size: 146.29 kB
Description: Excel file containing differentially expressed genes for all contrasts (one sheet per contrast). Includes all genes and pseudo-genes.

File: GE7218-allgenediff.xlsx
Size: 2.13 MB
Description: Excel file containing differential expression values for all genes in all contrasts (one sheet per contrast). Includes all genes and pseudo-genes.

File: GE7218-allExpressions.xlsx
Size: 2.98 MB
Description: Excel file containing normalized (RSEM) expression values for all genes in all conditions.

8. Differential expression - isoform level
The following table reports the number of differentially expressed isoforms in each contrast with abs(log2(FC)) >= 1.0 and FDR-corrected P-value <= 0.05. The lists of differentially expressed isoforms for all contrasts can also be downloaded as a single Excel file using the link below.

Test	Control	Tot isoforms	Overexpressed	Underexpressed	Table	Expressions
HET	WT	363	123	240	HET.vs.WT.isoDiff.csv	HET.vs.WT.imatrix.csv
MUT	WT	2,032	890	1,142	MUT.vs.WT.isoDiff.csv	MUT.vs.WT.imatrix.csv

Table 9. Results of isoform-level differential expression analysis.

File: GE7218-isodiff.xlsx
Size: 208.41 kB
Description: Excel file containing differentially expressed isoforms for all contrasts (one sheet per contrast).

File: GE7218-allisodiff.xlsx
Size: 7.63 MB
Description: Excel file containing differential expression values for all isoforms in all contrasts (one sheet per contrast).

9. Differential expression - combined files
The following file contains merged differential expression data. The first sheet contains fold changes for all genes that were found to be differentially expressed in at least one contrast. The second and third sheets contain the same information for coding genes only, and all transcripts.

File: GE7218-merged.allDiff.xlsx
Size: 306.20 kB
Description: Merged fold changes for all differentially expressed genes, coding genes, and transcripts respectively.

10. Alternative splicing analysis
Alternative splicing analysis was performed using rMATS version v4.1.0. The following table reports the number of events in each class for each contrast. The link in the last column allows you to download an Excel file containing full results.

Test	Control	SE	RI	MXE	A3SS	A5SS	Full
HET	WT	24	9	9	9	5	HET.vs.WT.MATS.xlsx
MUT	WT	8	8	6	3	3	MUT.vs.WT.MATS.xlsx

Table 10. Number of alternative splicing events, by class, for each contrast. Classes are: SE=exon skipping; RI=intron retention; MXE=mutually exclusive exons; A3SS=alternative 3' splice site; A5SS=alternative 5' splice site.

11. MultiQC report
MultiQC is a general Quality Control tool for a large number of bioinformatics pipelines. The report on this analysis (generated using MultiQC version 1.12) is available here:

MultiQC report

12. UCSC hub

UCSC Genome Browser: use the previous link to display the data tracks automatically, or copy the the URL https://lichtlab.cancer.ufl.edu/reports/H2B//GE7218/hub/hub.txt and paste it into the "My Hubs" form in this page.

WashU EpiGenome Browser: use the previous link to display the data tracks automatically, or copy the following URL into the "Datahub by URL Link" field: https://lichtlab.cancer.ufl.edu/reports/H2B//GE7218/hub/hub.json.

13. Methods summary

Trimming and QC on short reads were performed by fastp (v 0.22.0) [1].

The reads were aligned to the transcriptome using STAR version 2.7.9a [2].

Transcript abundance was quantified using RSEM (RSEM v1.3.1) [3].

Differential expression analysis was performed using DESeq2 [4], with an FDR-corrected P-value threshold of 0.05. The output files were further filtered to extract transcripts showing a 2.0-fold change in either direction. Results were reported for protein-coding genes only, and for all transcript types.

Alternative splicing analysis was performed using rMATS version v4.1.0 [5].

References

Shifu Chen, Yanqing Zhou, Yaru Chen, Jia Gu; fastp: an ultra-fast all-in-one FASTQ preprocessor, Bioinformatics, Volume 34, Issue 17, 1 September 2018, Pages i884–i890, https://doi.org/10.1093/bioinformatics/bty560
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 29(1):15-21 | doi: 10.1093/bioinformatics/bts635
Li B and Dewey CN (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12:323 | doi: 10.1186/1471-2105-12-323
Love MI, huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15,550 (2014). | doi: 10.1186/s13059-014-0550-8
Shen S., Park JW., Lu ZX., Lin L., Henry MD., Wu YN., Zhou Q., Xing Y. rMATS: Robust and Flexible Detection of Differential Alternative Splicing from Replicate RNA-Seq Data. PNAS, 111(51):E5593-601 | doi: 10.1073/pnas.1419161111

Completed: 1-5-2024@12:03