0 Description of Workflow

Previous studies suggested that de novo assembly is beneficial for significantly differentially expressed genes (DEG) even when a reference genome is available. [1] Moreover, it is now widely recognized that performing transcript quantification and afterwards obtaining the gene expression by adding together the expression from the individual transcripts will result in improved gene-level analysis. [2, 3] Therefore, the current analysis for this bulk RNA-seq data set is based on the transcript-level with de novo transcriptome assembly. The current analysis protocol was modified from a guideline attached to a peer-reviewed bioinformatic software called IsoformSwitchAnalyzeR [4] (Please Click HERE to check the origianl protocol). Briefly, the current protocol can be described as a 6-step process:

  1. After QC and potential quality trimming, the reads were mapped to the genome with the STAR aligner.

  1. Then the de novo transcriptome assembly was performed with StringTie software and was guided by existing annoation from ENSEMBL database.

  1. The third step is to merge the assembled transcripts from all individual StringTie runs into one combined reference transcriptome with StringTie software. The merged GTF file was further annotated with SQANTI3 software, which is designed to annotate the data generated from oxford nanopore sequencing. Meanwhile, the Circular RNA was detected with CIRCexplorer2 software [5] from the mapping results created in Step 1 based on the merged GTF file.

  1. To make the analysis of the individual samples comparable, all samples were re-quantified using the merged GTF file created in Step 3 with a quasi-mapping aligner called Salmon.

  1. Applying a “rescue” algorithm embbed in the IsoformSwitchAnalyzeR [4] software to the merged GTF file, to ensures that all novel isoforms are assigned to genes and that “reference gene_id” and “reference gene_names” are used instead of “Stringtie gene_is” for all already annotated genes.The “rescued” GTF file was annotated by SQANTI3 software again and the still remained novel genes (“Stringtie genes” not overlapping any “Refrence genes”) was removed for the downstream analysis.

  1. Finally, after correction of batch effect by RUVSeq software, the significant differentially-expressed genes (DEGs) were detected by DESeq2 software; the significant differential exon usage (DEU) were detected by DEXseq software; the analysis of isoform switches with functional consequences and the associated alternative splicing was performed using IsoformSwitchAnalyzeR software; the significant differentially-expressed circular RNAs were detected by circRNAprofiler software.

The source code is avaible on the Supplementary section to reproduce this analysis.

In addition, other data set was generated in order to use the leafcutter software to check the alternative splicing at nucleotide level; to use the Ularcirc software to visulize, to predict ORF, and to estimate the potenial biological function of detected circular RNAs.

The above workflow is developed and mantained by Bioinformatics Study Group in Okayama University (BSGOU). BSGOU is an international academic community committed to advancing the digital transformation of biological and biomedical research. We unite students, researchers, clinicians, and engineers to collaboratively explore how high-throughput data and integrative computation can drive new theories, models, and discoveries in life sciences. For more information in details, please visit the homepage (https://labonom.github.io/) of BSGOU.

[1] Wang S, Gribskov M. Comprehensive evaluation of de novo transcriptome assembly programs and their effects on differential gene expression analysis. Bioinformatics. 2017 Feb 1;33(3):327-33.

[2] Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research. 2015;4.

[3] Yi L, Pimentel H, Bray NL, Pachter L. Gene-level differential analysis at transcript-level resolution. Genome biology. 2018 Dec;19(1):1-1.

[4] Vitting-Seerup K, Sandelin A. IsoformSwitchAnalyzeR: analysis of changes in genome-wide patterns of alternative splicing and its functional consequences. Bioinformatics. 2019 Nov 1;35(21):4469-71.

[5] Zhang XO, Dong R, Zhang Y, Zhang JL, Luo Z, Zhang J, Chen LL, Yang L. Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. Genome research. 2016 Sep 1;26(9):1277-87.

1 Sample Sequencing Statistics

1.1 Reference

Annotation Files

Help

The de novo assembly transcriptome and mapping results could be visualized using IGV software. The detials could be checked below:

1.2 Quality Control Report

Click HERE to check the FastQC results.

Click HERE to inspect the de novo transcriptome assembly results.

1.3 Quasi-Mapping by Salmon software

*Mapping Rate (%): the number of reads mapped on genome features (i.e. a genome region contains genes)

Automatically detected most likely library type as U.

Click HERE to learn more detials about the fragment library types.

Fragment Library Types

1.4 Correction of Batch Effect by RUVSeq

Before Correction

RLE plot

RLE plot: Relative log expression plots.

PCA plot

After Correction

RLE plot

RLE plot: Relative log expression plots.

PCA plot

2 Analysis

2.1 Basic Information

Experiment Design

Name Group Numbering Description
No-inf_S17_L001 No-inf No-inf_S17_L001 No-inf_S17_L001
No-inf_S17_L002 No-inf No-inf_S17_L002 No-inf_S17_L002
No-inf_S17_L003 No-inf No-inf_S17_L003 No-inf_S17_L003
No-inf_S17_L004 No-inf No-inf_S17_L004 No-inf_S17_L004
pgKDN-inf_S19_L001 pgKDN-inf pgKDN-inf_S19_L001 pgKDN-inf_S19_L001
pgKDN-inf_S19_L002 pgKDN-inf pgKDN-inf_S19_L002 pgKDN-inf_S19_L002
pgKDN-inf_S19_L003 pgKDN-inf pgKDN-inf_S19_L003 pgKDN-inf_S19_L003
pgKDN-inf_S19_L004 pgKDN-inf pgKDN-inf_S19_L004 pgKDN-inf_S19_L004
pgwt-inf_S18_L001 pgwt-inf pgwt-inf_S18_L001 pgwt-inf_S18_L001
pgwt-inf_S18_L002 pgwt-inf pgwt-inf_S18_L002 pgwt-inf_S18_L002
pgwt-inf_S18_L003 pgwt-inf pgwt-inf_S18_L003 pgwt-inf_S18_L003
pgwt-inf_S18_L004 pgwt-inf pgwt-inf_S18_L004 pgwt-inf_S18_L004

Parameter Setting

  • The genes with read number > 10 in at least 3 samples were defined as expressed genes.
  • Differentially Expressed Genes (DEGs) were defined as FDR < 0.05 and |log2FC| > 1.
  • DEGs testing method in DESeq2 software was “LRT”.
  • Hierarchical cluster analysis for PCA and Heatmap plots was “ward.D2” based on “pearson” distance.
  • Significantly enriched GO terms in GSEA analysis were defined as P-value < 0.05.
  • Significantly enriched KEGG terms in GSEA analysis were defined as P-value < 0.05.
  • Significantly enriched GO terms in Over Representation Analysis were defined as P-value < 0.05.
  • Significantly enriched KEGG terms in Over Representation Analysis analysis were defined as P-value < 0.05.
  • Transcripts with isoform fraction (IF) > 5% were defined as expressed transcripts (only validate when using IsoformSwitchAnalyzeR).
  • Significant isoform switching were defined as FDR < 0.05 and dIf > 0.05.
  • Differentially Expressed Exons were defined as FDR < 0.1.
  • Differential intron excision in leafcutter software were defined as FDR < 0.05.

2.2 Differential Gene Expression Analysis

Click HERE to check the read counts, TPM, and feature length of all mapped genes across all samples in a Microsoft .excel file. Click HERE to check all parameters of DESeq2 model for all detected genes of all samples in a Microsoft .excel file. Click HERE to check the DESeq2 model estimated mean gene expression level of each condition (the intercept of a linear model) and its standard errors. For more detials of the DESeq2 model, please Click HERE.

Overview of Differentially Expressed Genes (DEGs)
    Sig. AS
    Sig. DEGs        in sig. DEGs
All
Detected
Genes		 Non Sig.
DEGs		     Up-
DEGs		 Down-
DEGs		     Genes with
Sig. AS		 in Non Sig.
DEGs		     in Up-
DEGs		 in Down-
DEGs
pgwt-inf vs No-inf 14951 11084     1671 2196     466 323     81 62
pgKDN-inf vs No-inf 14760 11062     1331 2367     2071 1622     211 238
pgKDN-inf vs pgwt-inf 13988 13558     136 294     382 363     7 12
† TPM, Transcripts Per Kilobase Million; DEGs, Differentially-Expressed Genes; Sig., Significant; Up-DEGs, Significantly Up-Regulated Genes; Down-DEGs, Significantly Down-Regulated Genes; AS, Alternative Splicing

2.3 Functional Enrichment Analysis

Number of Significantlly-Enriched Terms of Functional Enrichment Analysis (FEA)
GSEA    ORA
    GO-BP    GO-MF    GO-CC    KEGG
BP MF CC KEGG     From
All Sig.
DEGs		 From
Up-
DEGs		 From
Down-
DEGs		     From
All Sig.
DEGs		 From
Up-
DEGs		 From
Down-
DEGs		     From
All Sig.
DEGs		 From
Up-
DEGs		 From
Down-
DEGs		     From
All Sig.
DEGs		 From
Up-
DEGs		 From
Down-
DEGs
pgwt-inf vs No-inf 429 80 72 57     2236 2008 847     239 206 104     146 136 72     104 98 54
pgKDN-inf vs No-inf 430 94 74 63     1960 1797 763     217 182 104     129 114 67     101 87 49
pgKDN-inf vs pgwt-inf 123 26 13 18     800 328 722     62 33 61     34 31 26     48 0 47
† GSEA, Gene Set Enrichment Analysis; ORA, Over Representation Analysis; DEGs, Differentially-Expressed Genes; GO-BP, Gene Ontology Biological Processes; GO-MF, Gene Ontology Molecular Function; GO-CC, Gene Ontology Cellular Component; KEGG, Kyoto Encyclopedia of Genes and Genomes; Sig., Significant; Up-DEGs, Significantly Up-Regulated Genes; Down-DEGs, Significantly Down-Regulated Genes

2.4 Analysis of Isoform Switches

Click HERE to download all results of the significantly-switched isoforms analysis (q < 0.05 and dIF > 0.05).

Click HERE to download the alternative splicing analysis results for visualization by LeafCutter software.

Alternative Splicing Events

Amount of Event