DiBiG
ICBR BioinformaticsPowered by Actor, v1.0

RNAseq - Alignment and differential expression analysis

Title: GE6816
Project: (none)
Started on: 6/24/2023 23:29:06
Hostname: login1.ufhpc
Run directory: /blue/licht/runs/Melanoma-Project/GE6816/GE6816
Configuration GE6816.conf
Table of contents:
  1. Input data
  2. Trimming and quality control
  3. Alignment to transcriptome
  4. Genome coverage
  5. Expression analysis - quantification
  6. Differential expression - protein-coding genes
  7. Differential expression - all genes
  8. Differential expression - isoform level
  9. Differential expression - combined files
  10. Alternative splicing analysis
  11. MultiQC report
  12. UCSC hub
  13. 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).

CategoryData
Summary of input data
Experimental conditions:EV_1205Lu, HDAC8_1205Lu, EV_SKMEL28, HDAC8_SKMEL28
Contrasts:EV_1205Lu vs. HDAC8_1205Lu, EV_SKMEL28 vs. HDAC8_SKMEL28
Number of samples12
Sequencing data data
Total number of reads:425,829,632
Average reads per sample:35,485,802
Table 1. Summary of input data



ConditionSampleNumber of reads% Reads
EV_1205LuEV-1205Lu-134,302,4118.06%
EV-1205Lu-240,601,8509.53%
EV-1205Lu-328,894,7616.79%
HDAC8_1205LuHDAC8-1205Lu-131,065,2727.30%
HDAC8-1205Lu-231,690,3837.44%
HDAC8-1205Lu-330,855,3597.25%
EV_SKMEL28EV-SKMEL28-133,538,2067.88%
EV-SKMEL28-235,915,8588.43%
EV-SKMEL28-338,514,9009.04%
HDAC8_SKMEL28HDAC8-SKMEL28-147,333,34311.12%
HDAC8-SKMEL28-233,214,4877.80%
HDAC8-SKMEL28-339,902,8029.37%
Table 2. Number of reads in each sample.

2. Trimming and quality control
The input sequences were trimmed using trimmomatic. Quality control was performed before and after trimming using FastQC. The following table provides links to the quality control reports before and after trimming, as well as the number of reads in the trimmed files.

SampleReadsetReads before trimQC before trimReads after trimQC after trim% Retained
EV-1205Lu-1EV-1205Lu-1_r134,302,411EV-1205Lu-1_S1_L001_R1_001
EV-1205Lu-1_S1_L001_R2_001
32,558,283EV-1205Lu-1_S1_L001_R1_001.trim.paired
EV-1205Lu-1_S1_L001_R2_001.trim.paired
94.92%
EV-1205Lu-2EV-1205Lu-2_r140,601,850EV-1205Lu-2_S2_L001_R1_001
EV-1205Lu-2_S2_L001_R2_001
38,691,657EV-1205Lu-2_S2_L001_R1_001.trim.paired
EV-1205Lu-2_S2_L001_R2_001.trim.paired
95.30%
EV-1205Lu-3EV-1205Lu-3_r128,894,761EV-1205Lu-3_S3_L001_R1_001
EV-1205Lu-3_S3_L001_R2_001
27,592,976EV-1205Lu-3_S3_L001_R1_001.trim.paired
EV-1205Lu-3_S3_L001_R2_001.trim.paired
95.49%
HDAC8-1205Lu-1HDAC8-1205Lu-1_r131,065,272HDAC8-1205Lu-1_S4_L001_R1_001
HDAC8-1205Lu-1_S4_L001_R2_001
29,494,778HDAC8-1205Lu-1_S4_L001_R1_001.trim.paired
HDAC8-1205Lu-1_S4_L001_R2_001.trim.paired
94.94%
HDAC8-1205Lu-2HDAC8-1205Lu-2_r131,690,383HDAC8-1205Lu-2_S5_L001_R1_001
HDAC8-1205Lu-2_S5_L001_R2_001
29,995,342HDAC8-1205Lu-2_S5_L001_R1_001.trim.paired
HDAC8-1205Lu-2_S5_L001_R2_001.trim.paired
94.65%
HDAC8-1205Lu-3HDAC8-1205Lu-3_r130,855,359HDAC8-1205Lu-3_S6_L001_R1_001
HDAC8-1205Lu-3_S6_L001_R2_001
29,320,456HDAC8-1205Lu-3_S6_L001_R1_001.trim.paired
HDAC8-1205Lu-3_S6_L001_R2_001.trim.paired
95.03%
EV-SKMEL28-1EV-SKMEL28-1_r133,538,206EV-SKMEL28-1_S7_L001_R1_001
EV-SKMEL28-1_S7_L001_R2_001
31,992,820EV-SKMEL28-1_S7_L001_R1_001.trim.paired
EV-SKMEL28-1_S7_L001_R2_001.trim.paired
95.39%
EV-SKMEL28-2EV-SKMEL28-2_r135,915,858EV-SKMEL28-2_S8_L001_R1_001
EV-SKMEL28-2_S8_L001_R2_001
34,224,857EV-SKMEL28-2_S8_L001_R1_001.trim.paired
EV-SKMEL28-2_S8_L001_R2_001.trim.paired
95.29%
EV-SKMEL28-3EV-SKMEL28-3_r138,514,900EV-SKMEL28-3_S9_L001_R1_001
EV-SKMEL28-3_S9_L001_R2_001
36,612,774EV-SKMEL28-3_S9_L001_R1_001.trim.paired
EV-SKMEL28-3_S9_L001_R2_001.trim.paired
95.06%
HDAC8-SKMEL28-1HDAC8-SKMEL28-1_r147,333,343HDAC8-SKMEL28-1_S10_L001_R1_001
HDAC8-SKMEL28-1_S10_L001_R2_001
44,850,606HDAC8-SKMEL28-1_S10_L001_R1_001.trim.paired
HDAC8-SKMEL28-1_S10_L001_R2_001.trim.paired
94.75%
HDAC8-SKMEL28-2HDAC8-SKMEL28-2_r133,214,487HDAC8-SKMEL28-2_S11_L001_R1_001
HDAC8-SKMEL28-2_S11_L001_R2_001
31,629,578HDAC8-SKMEL28-2_S11_L001_R1_001.trim.paired
HDAC8-SKMEL28-2_S11_L001_R2_001.trim.paired
95.23%
HDAC8-SKMEL28-3HDAC8-SKMEL28-3_r139,902,802HDAC8-SKMEL28-3_S12_L001_R1_001
HDAC8-SKMEL28-3_S12_L001_R2_001
37,896,100HDAC8-SKMEL28-3_S12_L001_R1_001.trim.paired
HDAC8-SKMEL28-3_S12_L001_R2_001.trim.paired
94.97%
Table 3. Number of reads in input files and links to QC reports.

The following two tables report the number of reads before and after QC in each sample and in each condition.

SampleReads before QCReads after QC% Retained
EV-1205Lu-134,302,41132,558,28394.92%
EV-1205Lu-240,601,85038,691,65795.30%
EV-1205Lu-328,894,76127,592,97695.49%
HDAC8-1205Lu-131,065,27229,494,77894.94%
HDAC8-1205Lu-231,690,38329,995,34294.65%
HDAC8-1205Lu-330,855,35929,320,45695.03%
EV-SKMEL28-133,538,20631,992,82095.39%
EV-SKMEL28-235,915,85834,224,85795.29%
EV-SKMEL28-338,514,90036,612,77495.06%
HDAC8-SKMEL28-147,333,34344,850,60694.75%
HDAC8-SKMEL28-233,214,48731,629,57895.23%
HDAC8-SKMEL28-339,902,80237,896,10094.97%
Table 4. Number of reads in each sample before and after QC.



ConditionReads before QCReads after QC% Retained
EV_1205Lu103,799,02298,842,91695.23%
HDAC8_1205Lu93,611,01488,810,57694.87%
EV_SKMEL28107,968,964102,830,45195.24%
HDAC8_SKMEL28120,450,632114,376,28494.96%
Table 5. Number of reads in each condition before and after QC.

3. Alignment to transcriptome
The input sequences were aligned to the 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.

SampleInput readsGenome alignmentsGenome alignment rateTranscriptome alignmentsTranscriptome alignment rateAlignment report
EV-1205Lu-132,558,28378,471,3042.4128,159,08386.49%EV-1205Lu-1.star/Log.final.out
EV-1205Lu-238,691,65792,167,6742.3833,839,16887.46%EV-1205Lu-2.star/Log.final.out
EV-1205Lu-327,592,97666,514,2852.4123,949,03286.79%EV-1205Lu-3.star/Log.final.out
HDAC8-1205Lu-129,494,77875,073,6712.5524,569,42383.30%HDAC8-1205Lu-1.star/Log.final.out
HDAC8-1205Lu-229,995,34276,194,3842.5424,972,56383.25%HDAC8-1205Lu-2.star/Log.final.out
HDAC8-1205Lu-329,320,45673,694,6332.5124,623,72683.98%HDAC8-1205Lu-3.star/Log.final.out
EV-SKMEL28-131,992,82078,185,6622.4427,924,35987.28%EV-SKMEL28-1.star/Log.final.out
EV-SKMEL28-234,224,85781,474,9452.3830,396,97188.82%EV-SKMEL28-2.star/Log.final.out
EV-SKMEL28-336,612,77490,221,3602.4631,755,36986.73%EV-SKMEL28-3.star/Log.final.out
HDAC8-SKMEL28-144,850,606106,485,3732.3739,061,18287.09%HDAC8-SKMEL28-1.star/Log.final.out
HDAC8-SKMEL28-231,629,57876,436,8722.4227,226,45486.08%HDAC8-SKMEL28-2.star/Log.final.out
HDAC8-SKMEL28-337,896,10090,989,7022.4032,672,70486.22%HDAC8-SKMEL28-3.star/Log.final.out
Table 6. 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.

NameTotal ntCoverageEffective bpEffective PercEff Coverage
EV-1205Lu-1104,499,387,58134.20768,566,80225.20%135.97
EV-1205Lu-2121,678,809,33439.43790,023,21525.60%154.02
EV-1205Lu-383,709,586,24827.40739,568,96024.20%113.19
HDAC8-1205Lu-177,219,687,23525.03758,446,75824.60%101.81
HDAC8-1205Lu-280,379,304,19826.05762,156,59924.70%105.46
HDAC8-1205Lu-378,937,680,67725.58761,317,85724.70%103.69
EV-SKMEL28-192,336,021,78029.92755,508,59424.50%122.22
EV-SKMEL28-293,685,972,40830.36778,212,72225.20%120.39
EV-SKMEL28-3104,709,176,55033.94778,301,69425.20%134.54
HDAC8-SKMEL28-1129,014,983,07241.80827,557,02326.80%155.90
HDAC8-SKMEL28-290,519,182,16229.33763,257,71724.70%118.60
HDAC8-SKMEL28-3101,176,341,82332.79793,018,62625.70%127.58
Table 7. Genome coverage by sample.

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

NameTotal ntCoverageEffective bpEffective PercEff Coverage
EV_1205Lu297,797,193,13496.49921,746,26329.90%323.08
HDAC8_1205Lu230,328,045,43274.63928,265,16330.10%248.13
EV_SKMEL28282,529,001,84191.54923,926,72329.90%305.79
HDAC8_SKMEL2800.0000.00%0.00
Table 8. Genome coverage by condition

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

File: GE6816.cond.cov.xlsx
Size: 12.58 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: 5.32 MB
Description: Matrix of FPKM values for all genes in all samples.

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

File: genes.xpra.txt
Size: 4.24 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, 153.55 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.

TestControlTotalOverexpressedUnderexpressedTableExpressions
EV_1205LuHDAC8_1205Lu4,9072,6152,292EV_1205Lu.vs.HDAC8_1205Lu.codinggeneDiff.csvEV_1205Lu.vs.HDAC8_1205Lu.gmatrix.csv
EV_SKMEL28HDAC8_SKMEL28592275317EV_SKMEL28.vs.HDAC8_SKMEL28.codinggeneDiff.csvEV_SKMEL28.vs.HDAC8_SKMEL28.gmatrix.csv
Table 9. Results of gene-level differential expression analysis.

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

File: GE6816-allcodingdiff.xlsx
Size: 1.99 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: GE6816.g.deseq2norm.xlsx
Size: 3.10 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, 155.26 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.

TestControlTotalOverexpressedUnderexpressedTableExpressions
EV_1205LuHDAC8_1205Lu6,3363,2123,124EV_1205Lu.vs.HDAC8_1205Lu.geneDiff.csvEV_1205Lu.vs.HDAC8_1205Lu.gmatrix.csv
EV_SKMEL28HDAC8_SKMEL28731326405EV_SKMEL28.vs.HDAC8_SKMEL28.geneDiff.csvEV_SKMEL28.vs.HDAC8_SKMEL28.gmatrix.csv
Table 10. Results of gene-level differential expression analysis (all biotypes).

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

File: GE6816-allgenediff.xlsx
Size: 2.56 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: GE6816-allExpressions.xlsx
Size: 3.59 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.

TestControlTot isoformsOverexpressedUnderexpressedTableExpressions
EV_1205LuHDAC8_1205Lu000EV_1205Lu.vs.HDAC8_1205Lu.isoDiff.csvEV_1205Lu.vs.HDAC8_1205Lu.imatrix.csv
EV_SKMEL28HDAC8_SKMEL28000EV_SKMEL28.vs.HDAC8_SKMEL28.isoDiff.csvEV_SKMEL28.vs.HDAC8_SKMEL28.imatrix.csv
Table 11. Results of isoform-level differential expression analysis.

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

File: GE6816-allisodiff.xlsx
Size: 2.56 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: GE6816-merged.allDiff.xlsx
Size: 597.18 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.

TestControlSERIMXEA3SSA5SSFull
EV_1205LuHDAC8_1205Lu13645413018EV_1205Lu.vs.HDAC8_1205Lu.MATS.xlsx
EV_SKMEL28HDAC8_SKMEL2816375305229EV_SKMEL28.vs.HDAC8_SKMEL28.MATS.xlsx
Table 12. 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/MEL//GE6816/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/MEL//GE6816/hub/hub.json.

13. Methods summary

Short reads were trimmed using trimmomatic (v 0.36) [1], and QC on the original and trimmed reads was performed using FastQC (v 0.11.4) [2] and MultiQC [3].

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

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

Differential expression analysis was performed using DESeq2 [6], 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 [7].


References

  1. Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina Sequence Data. Bioinformatics, btu170.
  2. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
  3. Philip Ewels, Mans Magnusson, Sverker Lundin and Max Kaller (2016). MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics | doi: 10.1093/bioinformatics/btw354 | PubMed: 27312411
  4. 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
  5. 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
  6. 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
  7. 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: 6-24-2023@23:30
© 2023, A. Riva, University of Florida.