1 Introduction

SeqNjoy is a stand-alone application to generate and analyze differential gene expression results from RNA-Seq experiments. Analyses can be performed from the very initial step using the files containing short reads (also known as FASTQ files), or from an intermediate stage using alignments files (BAM/SAM files). The result of the SeqNjoy analyses are files containing annotated list of the genes (or other genetic features) differentially expressed between groups of samples in different conditions, and interactive HTML pages that allow the user to graphically analyze, filter and download selected genesets of interest.

This tutorial is intended to guide you through the process of analyzing your RNA-Seq data with SeqNjoy. For each step, we describe the input type of files required, the available parameters and the outputs the step yields. In addition, to make you more aware of what you are doing, we explain how are these files constituted, the meaning and consequences of each parameter and some key concepts you should be clear about when analyzing RNA-Seq data. The goal is that you become able to perform a complete RNA-Seq gene expression analysis on your own.

SeqNjoy has a series of minimal requirements. Before trying to install it, check if your computer meets these requirements.

2 Computer requirements and installation

      tar xzvf SeqNjoy_X.X.X-Linux.tar.gz

      cd SeqNjoy_X.X.X-Linux

      ./runSeqNjoy.sh

3 RNA-Seq in short

RNA-Seq is a high-throughput technique of sequencing that allows to detect and quantify the presence of RNA in biological samples, which is a proxy for detecting and quantifying gene expression, among other applications.

The “wet lab” part of the RNA-Seq process (Figure 3.1, top) starts with the “library preparation”, a process that comprises some steps: RNA extraction from the samples, filtering to keep only mRNA molecules and deplete other RNA species (ribosomal RNA, for example), reverse transcription from mRNA molecules to more stable cDNA molecules, fragmentation and size-selection of the fragments to get smaller chunks of nucleic acids (more convenient to deal with) and the sequencing of short reads obtained from the fragments.

The sequencer is the machine in charge of “reading” the nucleotide sequences corresponding to the mRNA transcripts, generating millions of “short reads” that are stored in a number of “FASTQ files”, one per sample (or biological replica) of your experiment. These FASTQ files are the input files of the RNA-Seq data analysis process. And you can perform this analysis by yourself using SeqNjoy!

Figure 3.1 displays an schematic view of the standard RNA-Seq pipeline to analyze mRNA expression in biological samples. The upper part describes the experimental part of the protocol. The yellowish area (bottom) contains the three main steps for the analysis of the RNA-Seq data of the sequenced samples.

Figure 3.1: Schematic view of the standard RNA-Seq pipeline to analyze mRNA expression

In brief, the bioinformatic analysis of a RNA-Seq experiment comprises three consecutive steps:

1. Alignment to reference genome: The short reads of the input FASTQ files are aligned to a FASTA file containing the genome sequence selected to align to. As a result, files with the aligned reads in BAM format are generated (one BAM file per each FASTQ file or pair of FASTQ, depending of the RNA-Seq type of sequencing: single-end or paired-ends sequencing).

2. Gene quantification: The aligned reads of each BAM file are mapped to a file containing the start-to-end positions of each gene (GFF file). For each gene, all reads mapped to it are counted. As a results, a “counts file” is generated per BAM file. A counts file contains a list of all the genes of the genome and the number of reads (“counts”) mapped to each gene.

3. Samples comparison / Differential expression calculation: The counts files of the samples are grouped into two sets: one set of counts files is assigned to the “Control” samples group and other set of counts files are assigned to the “Case” samples group. The aim of this step is to statistically compare gene by gene the expression of that gene in each group and determine if there is a significant differential expression between both groups and the magnitude of this difference. The final result of this step is a file containing a table with the results of the differential expression test per gene.

4 SeqNjoy framework in short

SeqNjoy utilizes user-friendly interfaces that allow you to perform a complete RNA-Seq expression gene data analysis in a short number of steps. If you have never performed this kind of analysis, we recommend reading RNA-Seq in short section as a starting point to get familiarized with the overall process.

Flowchart below (Figure 4.1) summarizes how SeqNjoy addresses the different parts of the analysis:

Figure 4.1: Flowchart of SeqNjoy

The Homepage is the entry point of the application (upper-left part of the flowchart). In this initial page is where you create your projects. A project is a container for one single study and comprises all the files used for a single RNA-Seq data analysis. The files of a project are a) the ones loaded as input (short reads files, gene feature files, genome reference file, functional annotations) plus b) all the files generated in the different steps of the analysis (alignments files, counts files and lists of differentially expressed genes as a result of comparing samples).

From the list of created projects registered in the Homepage, you can select any of them and access the Project’s Main Page, which is the central hub of the application. All the functionalities are launched from this page. From the Project’s Main Page you can:

• Load the input files needed for the analysis.
• Access the menus to perform the successive steps of the analysis.
• Browse and operate the list of files obtained in the successive analysis steps.
• Access some specific utilities that allow to inspect/analyze the content of the files.

As in a recipe you need to add some initial ingredients to start your dish, you need to load certain input files to start the analysis of your project. You should load one of the following two combinations of files for starting the processing in SeqNjoy:

• A number of short reads files (FASTQ format) + 1 Genome Sequence file (FASTA format) + 1 Genome Feature file (GFF3 or GTF format). This combination is the common input for a complete RNA-Seq data analysis.

  or,

• A number of alignment files (in SAM or BAM format) + 1 Genome Feature file (GFF3 or GTF format).

In the first case, when you load FASTQ plus GFF plus FASTA files, the analysis process (“the recipe”) comprises three successive steps: Align Reads, Quantify Genes and Compare Results, as described in the Figure 4.2:

Figure 4.2: Analysis steps in SeqNjoy starting from Short Reads (FASTQ) files

In the second case, when you load alignment files (BAM/SAM files) the workflow is similar. The only difference is that the initial “Align Reads” step is skipped (as you have already aligned reads in your BAM/SAM files…), see Figure 4.3:

Figure 4.3: Analysis steps in SeqNjoy starting from Alignments (BAM/SAM) files

As seen in Figure 4.1, there are some Utilities that you can use to inspect or analyze the content of files of certain types (Figure 4.4).

Figure 4.4: Utilities in SeqNjoy for file inspection/analysis

Figure 4.5 displays the whole picture of the functionalities of SeqNjoy to perform a RNA-Seq data analysis. All of them are launched from the “Project’s Main Page”. The RNA-Seq data analysis proceeds from top to bottom. Observe how the files being the output of some steps are the input of the following:

Figure 4.5: Functionalities in SeqNjoy.

In the next sections, we describe each step of the Analysis, each Utility, as well as the characteristics and use of the involved files. For FIESTA 2.0 graphical viewer you have available an specific tutorial accessible here.

5 First steps in SeqNjoy: creating and managing projects

Once you have launched SeqNjoy (clicking the SeqNjoy icon (Windows) or typing ./runSeqNjoy.sh in the console (Linux/MacOS)) you land in the Homepage (Figure 5.1). In the upper part, the page welcomes you and gives you a summary and important information relative to the application.

Figure 5.1: Homepage

The lower panel contains the list of Available Projects. As no project has been created yet, this list is displayed empty at first.

First of all, let’s define what a project is. A project is a container for one single study and comprises all the files used for a single RNA-Seq data analysis. These files are the ones uploaded as input plus all the files resulting from the analysis: alignments, counts, lists of differentially expressed genes.

Each project is identified by a project name that you provide in the text box located at the bottom of the page, as seen in Figure 5.2.

Figure 5.2: Creating your first project

After clicking the Start New Project button, you are directed to the Project’s Main Page to start working with the project (see the Main Page section).

Once you have one or more projects created, the Available Projects area of the Homepage contains a list with the name and date of creation of your projects (Figure 5.3).

Figure 5.3: List of started projects

At this point, you can perform two simple operations for a particular project by clicking the icons located at both extremes of the project name. If you click the “play” icon the project is opened and you are directed to the Main Page to start or continue your data analysis. On the other side, by clicking the “basket” icon the project is deleted from your system. IMPORTANT: if you delete a project, all its associated files (the input and the results) are deleted so the project can NOT be recovered.

6 Project’s Main Page: Loading files

Once the project has been created, you can start uploading your input files. First, open the project by clicking the “play” icon in the Homepage and access the Project’s Main Page of the application (Figure 6.1).

Figure 6.1: Project’s Main page

Now, click the Load Files button and select in your computer one or several files of the allowed file types. SeqNjoy supports the following file types and extensions:

If you are not familiarized with any of these file types, go to the section File Types to know more about the structure of the format and the function and details of each file type.

Figure 6.2 shows the application loading three files that contain a genome reference (FASTA) file, a gene features (GFF) file and an annotations (ANN) file:

Figure 6.2: Uploading files in SeqNjoy

The pop-up window reports the status of file uploading and post-processing by means of a progress bar and informative messages describing the step being performed. When all is OK, you should see a window with an “all in green” aspect like this (Figure 6.3):

Figure 6.3: Uploading process successful

Take into account that in the case of some large and/or complex files the upload may take a considerable amount of time.

7 Project’s Main Page: File List

After you have uploaded your input files, your Project’s Main Page shows the list of the loaded files. The list offers different information and a number of options to operate on them. From the Project’s Main Page you can also perform the successive steps of your RNA-Seq gene expression data analysis project.

This is the aspect of the Project’s Main Page after a standard RNA-Seq data analysis project has been performed (Figure 7.1):

Figure 7.1: Project’s Main Page

The Project’s Main Page is structured as follows.

In the top area the following functionalities are available:

• Load Files: Click here to load your input files for the analysis. See Loading Files section for a description of the process of loading and File Types in SeqNjoy for a full description of the allowed file types.
• Project Name: This text box informs you of the current project you are working with. You can change the project name at any moment by typing the new name in the text box.
• Align Reads: Click here to access the first step of the analysis (short reads alignment). See the corresponding section for more information.
• Quantify Genes: Click here to access the second step of the analysis (quantification of genes). See the corresponding section for more information.
• Compare Samples: Click here to access the third step of the analysis (differential expression between sample groups). See the corresponding section for more information.
• About: Click here to see Credits or obtain Help.

The rest of the Project’s Main Page is dedicated to the File List. The File List area contains information about all the files loaded by the user or generated by the processing with SeqNjoy, and provides a number of operations that the user can run on one or more selected files. The File List comprises the following buttons and columns, described from the left to the right:

• Select/Unselect: A column of check boxes. When you check one or more boxes, the corresponding file(s) is/are marked as selected. If you click a checked box, the file(s) returns to unselected status. In case you want to select all the files of the list, you can do in a more convenient way by clicking the uppermost check box (it also serves to “unselect all”).

  There are two possible operations you can execute with the selected files:

  • Download the selected file(s): you are asked the destination path of the selected files in your system. In the case of results files, a file with the script in R that generated this result is also provided. In the case of diffexp results, an additional third file is provided: an interactive FIESTA 2.0 html page that you can open with a web browser.
  • Delete the selected file(s): The files are deleted from the project. Note that files that are being used as input of other files can not be deleted. Be aware that deleted files can not be recovered`.

  Both operations, Download and Delete, are also available on a file-by-file basis in the two rightmost columns of the File List.

• Show buttons: You can filter the File List to hide or show the files having a particular file type(s). For example, if you click a colored “bam” button , all bam files are hidden from the File List and the bam button becomes grayish . Click again the grayish icon to return your bam files to the list. You can click multiple Show buttons to show/hide several file types at one time. Note that only file types with no members selected can be hidden.

  The numbers accompanying the file type on the button indicate the total number of existing files of that type and how many of them have been selected via the checkboxes (see Select/Unselect above).

• Type: It informs you of the type of the file. Refer the File Types in SeqNjoy section to know more about the file types and their characteristics.

• Details: This column contains diverse information about the file and provides access to additional utilities for some file types.

• File name: the background color of the file name indicates the file type (e.g red for FASTQ, green for FASTA, etc.)
• File content: e.g. “short reads”, “chromosomes”, “genes”, “features”, “alignments”… In addition, it displays also the number of items (e.g. “1000000 reads”). For the alignment files, it displays the result of the duplication removal, if applied.
• File origin: For some generated files, the files whose processing` gives rise to those files are displayed. (e.g. the file pH5-Rep1.bam is a product of processing the FASTQs 1M-SRR9336468_1.fq and 1M-SRR9336468_2.fq.
• File glimpse: Displays a list of the first identifiers contained in the file (e.g. chromosome identifiers, gene ids)
• Introns: Displays information calculated from the GFF file about the introns in the genome. Some of this information (minimum and maximum intron length) is useful to take full advantage of the splice-aware feature of hisat2 aligner.
• For some file types, the Details column displays a specific button that launches an specific utility. Each button is colored according to the file type of the selected file. These are the supported utilities:

  File Types	Operation	Tool	Icon
  fasta, bam	Genomic Browser	IGV
  gff	Genome features structure	Rgff
  diffexp	Visualize, plot and filter gene expression results	FIESTA 2.0

• Bytes: Size of the file in bytes.
• Date: Day and time of creation of the file.
• (Download file icon): Downloads the file in your computer (you are asked where).
• (Delete file icon): Deletes the file (Note: this action cannot be reverted).

8 Align Reads

Click the “Align Reads” button of Main Page to access this window. For performing this task, you inform a genome sequence (FASTA) file and the corresponding short reads files (FASTQ files) and obtain files with the alignments (BAM formatted) of the short reads contained in the FASTQ files to the reference genome (FASTA file).

Figure 8.1: Align Reads: selecting alignment method (aligner)

After choosing any of them you access a window where you have to fill the input FASTQ files and alignment parameters that are necessary to complete the task.

Figure 8.2: Window for the reads alignment. Input files are the genome reference (FASTA file) and the short reads files (FASTQ files). The shaded areas contains bowtie-specific or hisat2-specific parameters and therefore are differents according to the selected aligner.

The Align Reads to Genomes window has the following parts (see Figure 8.2):

• Select Genome sequence: select here one FASTA file among the listed in the drop-down menu. Usually, you select the genome reference of the species/variety the sequenced samples come from (or a close one), but you can select any uploaded genome reference that could be more convenient in your particular study.

• Select FASTQ files: Below the Genome Sequence you have a grid of drop-down menus arranged in rows and columns. To fill this part in, you must know that there are two modalities of RNA-Seq named single-end and paired-ends sequencing. Refer to the corresponding section to know about them if you need it.

  • In “Fastq 1” column you select: a) all your single-end FASTQ files, if your FASTQ file are of the single-end type or b) the files corresponding to the first member of the pair in the case of paired-ends FASTQ files (normally, the filename of these “read 1” files contains the suffix “_1” or similar).

  • “Fastq 2” column only applies in the case of paired-end sequencing. In that case, you have to put here the second member of the pair (normally, files with the suffix “_2” or similar in their filesnames). This second member of the pair must be selected side-by-side with the first ("_1“) member that was selected in”Fastq 1" column (as in Figure 8.2). Do not select any file in this column if all the FASTQ files in your experiment are of “single-end” type.

  • In addition, you have to fill in a Sample label for each Fastq 1 (single-end) or Fastq 1 plus Fastq 2 (paired-end) item you fill in. This Sample label will name the resulting alignment (BAM) file.

  At the end, the aspect of this area will be similar to one of the displayed in the Figure 8.2.

• Output files suffix (optional): The text you type in this box will be added to the name of all the output files (BAM files) of this step, immediately after the sample label and before the .bam extension, that is, “Sample_name.suffix.bam” (Example: “Sample1.aligned_with_bowtie2.bam”, “Sample2.aligned_with_bowtie2.bam”, etc.). If no “output files suffix” is typed, the BAM output files will be named simply with the sample name followed by the “.bam” extension (i.e. [Sample_name].bam).

• Alignment parameters: The shaded areas of the window contain the parameters that are used by the aligner of your choice, bowtie2 or hisat2. Some options are the present irrespective of the aligner you choose:

  • Trim reads from 5’ end / Trim reads from 3’ end: Previously to align reads, usually the FASTQ files are submitted to programs to assess the quality (QA) of the short reads. FASTQC (Andrews et al., 2012) is probably one of the most popular. In Figure 8.3 we can see an example of one of the plots that FASTQC provides and that represents the “quality per base sequence” (see here how to interpret this plot), i.e. an overview of the range of quality values across all bases at each position in the FASTQ file. It is not unusual that sequencing decreases the quality of the positions at the 3’ end (“right” end). In Figure 8.3 the quality drop is so marked from position 120 (approx.) onward that a trimming of this chunk on the 3’ extreme of the reads might be advisable.
    

    Figure 8.3: FASTQC Plot representing the quality per base sequence on a FASTQ file with quality drop in the 3’ extreme.

    Other reason to trim nucleotides from the reads ends may be the presence of adapters that for some reason have not been removed after the sequencing. The aim of the trimming is to remove these read chunks of very low quality or that do not correspond to the sequenced sample.

    With the Trim reads from 5’ end parameter you can set the number of nucleotides to be removed from the reads of the file at the 5’ end (i.e. “on the left”). With Trim reads from 5’ end you can set the number of nucleotides to be removed from the reads of the file at the 3’ end (i.e. “on the right”).

  • Max CPU Cores: bowtie2 and hisat2 support parallel processing, that is, it can use multiple CPU cores. Multi-core is enabled by default to use (up to) 4 cores. You can modify this value in the text box with a numeric value greater or equal to 1. If your computer has less CPU cores than the number typed here, the aligner program will use the maximum your computer has.

  • Remove duplicated alignments: In RNA-Seq protocols (especially if a PCR amplification step is included), there is a certain probability of introducing artificial generation of duplicated reads, that is, reads of the same size aligning exactly in the same chunk of the genome (same start and end positions in the genome, Figure 8.4). The presence of duplicated reads above a certain number may be an indication of the existence of such unwanted technical artifacts during the sequencing process. The application provides this parameter to remove reads that are duplicated a number of times equal of greater to a value (50 by default).

    

    Figure 8.4: Duplicated reads (in red).

• Specific options for bowtie2.

  We use the R package Rbowtie2 (Wei et al., 2018) to align with bowtie2 (Langmead & Salzberg, 2012).

  • Performance: Alignment methods search a reasonable trade-off between speed and sensitivity/accuracy in their results in order to reduce false detected alignments without missing too many true alignments. In bowtie2 you can adjust this balance by setting one of the following values the “performance” option: very-sensitive, sensitive (default value), fast or very-fast. The default value, sensitive, keeps a good balance between sensitivity and precision. In the case you need more precision (at the cost of a slower performance) choose very-sensitive. Conversely, choose fast or very-fast to prioritize speed performance (at the cost of a lesser precision).

  • Mode: By default, bowtie2 aligns “end-to-end” mode, meaning that it attempts to align the entire read without omitting characters at either extreme. The alternative, local mode, does not require that reads align from one extreme to the other, that is, it can yield alignments that do not include one (or both) of their extremes. Observe in the example of Figure 8.5 that in "local" mode the extremes are not aligned while they are both aligned in "end-to-end" mode.

    

    Figure 8.5: “end-to-end” and “local” modes.

• Specific options for hisat2.

  We use the R package Rhisat2 (Soneson, 2021) to align with hisat2 (Kim et al., 2015).

  All the functionalities of these hisat2-specific parameters are linked to the fact that hisat2 is a “splice-aware” aligner. In practical terms, that means that it uses the information deduced and/or provided about the genomic locations of the splice sites and introns in order to improve the precision of the predicted alignments. Of course, this not applies for intronless genomes. By informing these parameters the alignment may be slightly more accurate (if you do not inform them, hisat2 will do his best anyway.)

  • Min(inimum)/Max(imum) intron length: By default, 20 nucleotides and 100000 nucleotides, respectively. Set here the (estimated) length of the shortest/longest intron in your genome. If possible, this information is provided for you by SeqNjoy in the “Details” column (File List) of the corresponding GFF file.
  • GFF file (optional). You can select here one of the GFF files you loaded corresponding to this sequencing. From this GFF file SeqNjoy will try to extract the known splice sites of the genome, that will be used by hisat2 to improve the alignments involving splice sites.

Once you have filled the form in with the required files and set your parameter values, the “Launch!” button turns from disabled (grey) to enabled (red) so you can click it and launch the alignment(s). A pop-up window appears informing you of the progress of the alignment task. When the process ends OK a message “done” is displayed in the pop-up window and a “Close” button is enabled. By clicking “Close” you return to the “Align Reads to Genome” pane. Now, you could make a new Align Reads to Genome round or click the “X” icon for closing the “Align Reads to Genome” window. You will then return to the “Project’s Main Page”, which will have been updated to contain the new alignment (.bam) files in its File List.

Figure 8.6: Align Reads to Genome dialog ready to launch the alignments processes.

9 Quantify Genes

Click the “Quantify Genes” button in the “Project’s Main Page” to access this window (Figure 9.1).

Figure 9.1: Quantification window

In summary, you inform here a genome features file (GFF) file and a number of alignment files in BAM. Normally, these alignments have been generated previously in the “Align Reads” step. Alternatively, the BAM files were generated elsewhere and the user uploaded to SeqNjoy.

The output of the quantification is one “counts” file per BAM alignment file. Each “counts” file contains a list with the “counts” , i.e. the number of “reads” or “fragments” that have been mapped to each gene (or other feature type that you select) of the genome.

• Genome Features: Select one GFF3 or GTF file from the list of loaded GFF files displayed in the drop-down menu.

• Output file suffix (optional): you can type in this box a text chunk that will be added to the name of all the output files of this step (“counts” files), immediately before the .counts extension, that is, “bam_name.suffix.counts” (Example: “Sample1.exons_by_genes.bam”, “Sample2.exons_by_genes.bam”, etc.). If no output file suffix is typed, the resulting files will be named similarly to the corresponding BAM (“sample_name.counts”).

• Alignments: In this column of selectors you select the samples (their corresponding BAM alignment files) from which you want to obtain the quantification of counts for one (or several) genomic feature(s). At least one BAM file must be selected.

The following parameters are used by the tool that performs the quantification task, the function featureCounts of the R package Rsubread (Liao et al., 2019):

• Strand specificity: Select the type of protocol that was used in the library preparation: “direct”, “reverse” or “unstranded”.
  • With unstranded, a read will be counted for a gene irrespective to which strand it maps to.
  • With direct and having single end data, the read has to map to the same strand as the gene. For paired end data, the first read of a pair has to map to the same strand as the gene, and the second read has to map to the opposite strand.
  • With reverse and having paired end data, the second read has to map to the same strand as the gene, and the first read has to map to the opposite strand.

  If you have any doubt, there is a specific section on strand specificity in RNA-Seq in this Tutorial. No default value is set, as it depends on how the sequencing of your reads has been made. If you finally don’t know which strandness protocol was followed in the sequencing, select unstranded as it is the safest option.

• Primary only: A read may map ambiguously to multiple locations, e.g. repeated or closely similar regions in the genome. In that cases, only one of the multiple read alignments is considered “primary”, that is, “the best alignment” based on diverse criteria. All other alignments are marked as “secondary”. By default (box checked) the quantification will count only the alignments marked as “primary alignments” and discard the “secondary”. If not checked, all alignments will be counted. In Figure 9.2, if "primary only" is checked, only the alignment tagged as primary will be counted in feature A for read R001. If not, all the three alignments will be counted: one count for feature A, one for feature B and one for feature C.

  

  Figure 9.2: Primary and secondary alignments.

• Ignore Marked Dup[licate]s: This option applies only to alignment (BAM) files uploaded by the user and which have duplicate marking (with Picard Tools, for example). As said in a similar parameter of the “Align Reads to Genomes” dialog, duplicated reads are sometimes generated due to technical issues (Figure 8.4). By default, duplicated reads are not considered in the quantification (checked box). If you uncheck the box, the duplicated reads will be counted.
• Max CPU Cores: featureCounts supports parallel processing, that is, it can use multiple CPU cores. Multi-core is enabled by default to use (up to) 4 cores. You can modify this value in the text box with a numeric value greater or equal to 1. If your computer has less cores than the number typed here, the quantification program will use the maximum your computer has.

• Genomic features and blocks: The GFF file selected has a number of feature types that you can quantify by. For example, you may wish to get the list of reads that align (counts) per each gene. For that, in the drop-down menu “Select Feature” you should select the “feature type” gene among the existing in your GFF file. Some features, like gene, are constructs that have constituent sub-features underneath. One key sub-feature related to gene expression is exon. In principle, it is more accurately (specially if your GFF is well annotated) to take into account for quantification not all the reads aligning the gene span but only those reads aligning within the (block of) exons of the gene (Figure 9.3).

  

  Figure 9.3: Genomic features and blocks for computing counts in ‘GeneX’. When only the feature ‘gene’ is selected, all reads aligning in the GeneX span are counted (13 reads). When you add the block ‘exon’ to the feature ‘gene’ you are imposing the restriction that only the reads inside exons are counted, so the counts for GeneX are 10.

• Add more button: Click to set additional feature types to quantify. For example, if the GFF file contains “genes” and also “pseudogenes” you may want quantify both.

Once you have filled the form with the required files and set your parameter values, the “Launch!” button turns from disabled (grey) to enabled (red) so you can click it and launch the quantification. In Figure 9.4, the quantification process will be launched to count genes for nine samples, with reverse strandness mode. Only primary alignments will be counted. The output of the quantification is a number of counts files equal to the number of input BAM files, that is, nine in this case.

Figure 9.4: Quantification dialog ready to be launched

A pop-up window appears informing you of the progress of the quantification task. When the process ends OK a message “done” is displayed in the pop-up window and a “Close” button is enabled. By clicking “Close” you return to the Quantify Genomic Features pane. Now, you could make a new Quantify Genomic Features round or click the “X” icon for closing the window. You will then return to the “Project’s Main Page”, that it will have been updated to contain the new “counts” files.

10 Compare Samples (Differential Expression)

In this step you can perform an statistical analysis of the differences in the expression for a given set of features between two groups of samples. You can use it repeatedly to make more that one binary comparison (e.g. “SampleGroupDrugA versus SampleGroupDrugB”, and then “SampleGroupDrugA versus SampleGroupPlacebo” and finally “SampleGroupDrugB versus SampleGroupPlacebo”).

Figure 10.1: Compare Samples: selecting the method for calculating differential expression

After selecting the desired option, the Calculate Differential Expression dialog is displayed (Figure 10.2). Most of this window is the same regardless of the selected software. Only some specific parameters located in the shaded area of the window are different and specific of the selected method (DESeq2 or edgeR) to perform the differential expression calculations.

Figure 10.2: Calculate Differential Expression (Compare Samples) dialog as displayed initially.

The following fields are present in the Calculate Differential Expression window irrespective of the selected method:

• Counts Group: Select one “counts group” from the listed in the drop-down menu. A “counts group” is a set of counts files linked a) to a particular GFF AND b) to a particular “feature” or “block-by-feature” counting decided for that GFF in the “Quantify Genes” step. This is expressed in the options displayed in the “Counts Group” drop-down menu with the syntax “gff_file.gff feature (block)” (or “gff_file.gff feature” if no block was selected in “Quantify Genes” step).
• Case Label and Control Label: Put here a denomination for the two groups of samples (or biological replicas) to compare, one for the “case”(or “test”, “mutated”, etc.) group and one for the “control” group. These labels are parts of the name of the results file, that has the following syntax: “[CaseLabel]vs[ControlLabel].[Output_suffix_if_typed].xlsx”.
• Select Replicate files: Here you have two columns with drop-down menus. Select from the counts group the counts files corresponding to the samples (or biological replicas) you wish to include in the “case” condition (left column) and do the same operation with the samples included in the “control” condition (put them in the right column). To perform the “differential expression calculation”, it is required that you select at least two samples (or biological replicas) in one of the conditions (but take into account that this minimum valid number of samples is extremely poor in statistical terms as input for the differential expression calculation).
• Output files suffix (optional): You can type in this box a text label that will be added to the name of the results file, immediately after the “[CaseLabel]vs[ControlLabel]” tag and before the “.xlsx” extension specific for differential expression files, that is, “[CaseLabel]vs[ControlLabel].[output_suffix].xlsx”. If no output files suffix is typed, the resulting file will be named simply “[CaseLabel]vs[ControlLabel].xlsx”.
• Annotations (optional): You can add here a file with annotations of your genes (or the specific feature type of the counts files being analyzed). As you can see in the corresponding section, it must match these requirements:
  • Be a tab-separated text file, with extension .ann, .ann.txt, .ann.tsv, or an Excel file, with extension .ann.xls or .ann.xlsx.
  • Have a first column containing (gene) identifiers, that must be of the same class as the identifiers present in the counts file (for example, if the gene IDs are ENSEMBL IDs in the counts files, the gene IDs must be ENSEMBL IDs in the first column of the user-provided annotations file).
  • Have a first column containing (gene) unique identifiers, that must be the same as the identifiers present in the counts file. Regarding your annotations file, take into account that different databases can assign different identifiers to the same feature (for example, ENSEMBL gene IDs are the format ENSG…. and Genbank gene accessions with the format NM…..).
  • The rest of columns contain a number of data relative to the respective (gene) identifier (e.g. name, description, genomic location, functional annotations, other IDs, etc), one data per column.
    
    

• Specific options for DESeq (Love et al., 2014):

  • Independent Filtering: The default (checked box) is to perform independent filtering. The purpose of independent filtering is to run a process, subsequent to and independent of the particular statistical tests performed by DESeq2, in order to remove from the results those that have no, or little chance of showing significant evidence, what typically increases the statistical/detection power. The simplest criterion to do the filtering is using the mean of all the normalized counts for the evaluated gene. In the box is unchecked, this additional filtering is disabled. More information in this link

  • Use beta prior The default (checked box) is to use beta prior, which is equivalent to perform a “Log fold change (LFC) shrinkage”, that is, to shrink the LFC estimates of genes with little information or high dispersion toward more likely (lower) LFC estimates. Figure 10.3 illustrates the effect of using (right) or not (left) beta prior (PC) on the shrinkage of LFC. Observe how using PC a shrinkage on the LFC (“Log2Ratio”, in the Y-axis) is performed. This shrinkage affects those genes that have low expression (low “Log2BaseMean”, X-axis). Red and green dots represent genes upregulated and downregulated, respectively, with p-adjusted (FDR) < 0.05.

    

    Figure 10.3: Effect of beta prior on the shrinkage of LFC.

    Use beta prior procedure when your data have a high number of low-expressed genes (low number of counts) and/or a high number of genes showing high difference of expression within a condition (high dispersion in the counts values). You can find additional explanations in this link.

• Specific options for edgeR (Robinson et al., 2009):

  • Prior Count (optional): This argument is the counterpart in edgeR to the Use beta prior argument in DESeq2. When checked, you can type in the text box a positive number that will be used in the calculations of differential expression as a factor to shrink Log Fold-Changes (LFC). The default is 0.125, a value that produce a minimal LFC shrinkage. Larger values of prior count will produce more shrinkage. Figure 10.4 shows the effect of different values of Prior Count (PC) on the shrinkage of LFC. Observe how the higher the value of PC the larger the shrinkage of LFC (“Log2Ratio”, in the Y-axis). Observe also how the shrinkage affects only the those genes with low expression (low “Log2BaseMean”, X-axis). Red and green dots represent genes upregulated and downregulated, respectively, with p-adjusted (FDR) < 0.05.

    

    Figure 10.4: Effect of different values of Prior Count on the shrinkage of LFC.

Once you have filled the form with the required files and set your parameters values, the Launch! button turns from disabled (grey) to enabled (red) so you can press it and launch the differential expression calculation. In the figure 10.5 differential expression calculation between the Case ("pH5_CO2") group of samples and the Control ("pH5") group of samples is displayed ready to be launched (Launch! button enabled). In the example, an annotations file has been provided so the information contained in this file will be incorporated to the diffexp file with the differential expression calculations.

Figure 10.5: Calculate Differential Expression dialog ready to be launched.

After launching the process, a pop-up window informs you the progress of the “differential expression” task. When the process ends OK the message “done” is displayed in the pop-up window and a close button is enabled. By clicking the Close button you return to the Calculate Differential Expression pane. Now, you can fill in again the form with a new comparison between groups of samples to calculate the differential expression between them or click the “X” icon for closing the window. You will then return to the “Project’s Main Page”, that it will have been updated to contain the newly created diffexp files in Excel format (extension .xlsx) containing the results of the calculated differential expressions of the comparisons performed.

The “diffexp” output file (syntax: [CaseLabel]vs[ControlLabel].[OutputSuffix_if_typed].xlsx) contains the list of differentially expressed genes (DEG) detected between the “Case” and the “Control” groups of samples. Check in the File Types section the information content of the diffexp files. You can graphically inspect and filter the differentially expressed genes results using the FIESTA 2.0 analytic tool.

11 Utilities

In addition to the pipeline to perform the RNA-Seq gene expression data analysis, SeqNjoy provides three utilities for inspecting the content of some of the types of files involved in the analysis. These utilities are launched from buttons located in the File List.

11.1 GFF feature structure hierarchical tree

To access this utility from a GFF (GFF3 or GTF formatted) of the File List click the icon located rightmost in the Details column.

When launched, SeqNjoy will run the plot_features function of the Rgff R package (Garcia-Martin et al., 2022) to generate a picture containing a tree representing the hierarchical structure (dependencies) of the feature types included in the selected GFF file (Figure 11.1). This tree representation is useful for a global look of what is annotated in the genome and it can help you decide which feature(s) and blocks to use in the Quantify Genes step of the analysis.

Figure 11.1: A hierarchical tree of a genome feature file (GFF) produced by Rgff

As it is shown in the example tree above, the highest feature type in this GFF is gene and there are 7127 items annotated as genes. gene has a child node, the transcript feature. Included into transcript there are five child features, that are siblings from each other: CDS, exon, five_prime_UTR, start_codon and stop_codon.

With this information, we know that we could count genes, exons-by-gene, transcript, exons-by-transcript, exons, etc. in our analysis but not “exon-by-CDS” for example, as exon is not a descendant of the CDS feature.

11.2 Visual exploration of genomic files with IGV

You can access this utility from a FASTA file containing the genomic sequence or from an alignment file (BAM), by clicking the icon located on the right of the Details column of the corresponding file. The option is not available for user-provided BAMs.

After clicking the button, the widely used “Integrative Genomic Viewer” (IGV) (Robinson et al., 2011, 2022) tool will be launched and will load the corresponding FASTA or BAM file. If you have loaded a FASTA, a screen similar to the one in Figure 11.2 will be displayed.

Figure 11.2: IGV launched with a FASTA file as input

If you use a BAM file with IGV a similar window is displayed. In addition to the BAM file, the associated genomic sequence (FASTA file) is automatically incorporated to IGV in the same session. See in Figure 11.3 the initial aspect displayed in this case and a minimal “use case” for IGV to visually compare the alignment files of two samples (i.e. two BAM files) for a particular gene.

Figure 11.3: A simple usage example of IGV to compare two BAM files in a certain genomic region.

IGV has a plethora of functionalities and possibilities to interact with your alignment files. Please refer the IGV user general guide, to take full advantage of this versatile tool. In addition, you can inspect your files with the web/on-line version that IGV provides at https://igv.org/app/, documented at https://igvteam.github.io/igv-webapp/.

11.3 Graphical analysis of differential expression results with FIESTA 2.0

FIESTA 2.0 (Oliveros, 2022) is an interactive viewer for differential expression results files that is integrated in SeqNjoy. It provides multiple ways of representing your data as tunable interactive plots that can be downloaded as high resolution images. In addition, FIESTA 2.0 allows you sorting and filtering genes by one or several columns and saving the relevant ones into files.

Figure 11.4 illustrates a usage example for FIESTA 2.0. When a diffexp file is loaded into FIESTA 2.0, a default representation of the expression data as a MA-plot is displayed. From this point, you can combine filters (top) to select gene subsets. A very common (and useful) filter combination is to select and highlight with user-selected colours genes (center) that show to be significantly (P-adjusted < 0.05) under- (negative Log2Ratio) or overexpressed (positive Log2Ratio). Finally, user can select any pair of numerical columns present in the diffexp file to be represented (bottom), so other popular representations like Volcano plots are also available. User-generated plots and subsets of relevant genes are both downloadable (in the case of the plots, with high resolution).

Figure 11.4: FIESTA 2.0 usage example.

To access FIESTA 2.0 viewer from a file containing the differential expression results (“diffexp” files in .xlsx format) click the icon located rightmost in the Details column.

More details about FIESTA 2.0 here.

12 File Types in SeqNjoy

SeqNjoy uses up to seven types of files in the analyses. The files that you provide have these types: FASTA, FASTQ, GFF3/GTF and ANN. Other types are only held by files generated through the SeqNjoy analysis: COUNT, DIFFEXP. Finally, the format BAM/SAM is that of the alignment files, which can be either internally generated by SeqNjoy or provided (loaded) by the user.

To follow, we describe the aspect, source, uses, format requirements and content of each file type. You can find a good complement to know more about them in this link by Andrew Severin and also in this page.

12.1 FASTA file

• File type and extensions: plain text file, uncompressed (fa, fasta, fna) or gzip-compressed (fa.gz, fasta.gz, fna.gz), e.g. my_fasta_file.fa, my_fasta_file.fasta.gz
• Source: User-provided. In case you do not have a FASTA file of the sequenced species/variety, you have good repositories of FASTA files in Ensembl, EnsemblGenomes, USCS and NCBI. Very often, you will find or have “.gz” FASTA files. “gz” files are compressed to save space. SeqNjoy will process them normally, you don’t need to uncompress them before loading.
• Use: As input file In the “Align Reads” task. It is used as genomic reference sequence where the short reads (FASTQ files) are aligned to.
• Content: FASTA files have a very simple format to contain the sequence of nucleotides of a genome, composed by one or several genomic sequences. Each genomic sequence item (a chromosome, for example) is preceded by a header line that starts by a “>” character and contains the identifier of that genomic sequence (Figure 12.1)

Figure 12.1: FASTA file format

12.3 Gene Feature File (GFF): GFF3 and GTF formats

• File type and extensions: plain text files, uncompressed (.gff3, .gtf) or gzip-compressed (.gff3.gz, .gtf.gz). e.g.: my_gene_features.gff3, my_gene_features.gff3.gz.
• Source: User-provided. GFF3 or GFF files are obtained very often in the same place that the file with the genome (FASTA file). This makes sense as a particular “Gene Features File” (GFF) is linked to a specific “FASTA reference file”. If you don’t have a GFF file for your sequenced species/variety you can try to find one in repositories like Ensembl, EnsemblGenomes, USCS or NCBI. Very often, the GFF3/GTF files are compressed (“.gz” extension) to save space. SeqNjoy will process them normally, you don’t need to uncompress them before loading.
• Use: As input file for the “Quantify Genes” task. In a simplistic explanation, in Quantify Genes step the aligned reads contained in each alignment (BAM) file are mapped to the GFF3/GTF file to count the number of reads that map within the limits of a each gene (or other selected feature type) contained in the GFF3/GFT file.

• Content: A GFF3 or GTF file contains information for the so-called genomic features. A feature is a genomic region with some involvement in a biological process. For example, the mRNA “EDEN.1” (mRNA00001) of the gene “EDEN” (gene00001) in the Figure 12.3 is a feature. In the same figure, that specific feature is described in Line 3 of the exemplified GFF file:

  

  Figure 12.3: Canonical representation of the representation of a single gene in GFF3. Adapted from the GFF3 specifications

  The Figure 12.3 depicts how a gene and its features composition (mRNAs, exons, CDSs) are represented in a GFF3. After some header lines (starting with “#”), each feature is written in a line. You can see that each feature has information about its location (start, end, strand) and feature type (gene, mRNA, exon, etc.). A key point is the last column, where relations of dependency between the features are coded. For example, the feature “mRNA00001” (with feature type “mRNA” ) in Line 3 depends on (note its “Parent” tag) another feature named “gene00001” (aka EDEN) in Line 2, which has the feature type “gene”. SeqNjoy has a specific utility for displaying the feature structure of your GFF3/GTF file in a graphical way.

More formally, each line of a GFF3 (or GTF) file has 9 tab-separated columns:

• Column 1: seqID (e.g. chromosome-id: “chr1” or “1”, “chr2”, etc..). IMPORTANT: These ids must be the same as the used in the FASTA file. See in the example below (Figure 12.4) how only the GFF3 in green has matching chromosomes ids (chr1, chr2…) with those in the FASTA file displayed.

  

  Figure 12.4: GFF and FASTA chromosome ids must match.

  • Column 2: Source (where the feature comes from: a program, a database…) Not important for us…
    
  • Column 3: Feature type (gene, mRNA, CDS, exon, etc.)
    
  • Column 4: Start position of feature in the genome (i.e. in the FASTA file)
  • Column 5: End position of feature in the genome (i.e. in the FASTA file)
    
  • Column 6: Score (confidence value). Often this field contains only a dot ‘.’ (uninformed). Anyway, not important for us.
    
  • Column 7: Strand (+,-,.). A dot means uninformed.
    
  • Column 8: Phase (0,1,2 or not informed (.) ).
  • Column 9: List of attributes, with a format tag=value per attribute. Multiple attributes are separated by “;”. It holds (among other data) the dependency relations of the feature.

GFF3 is currently the preferred format to represent Gene Features, but it is not uncommon to find gene feature files in GTF format. GFF3 provides more flexibility to represent the hierarchy of the features in a genome but GTF has some advantages in certain cases. SeqNjoy can process files in GFF3 or GTF formats. In both cases, columns 1 to 8 are the same. The 9th column has a different syntax in GFF3 and GTF (Figure 12.5 )

Figure 12.5: Comparison of GTF and GFF3 files for the representation of the features of the same gene

Important: For SeqNjoy, GFF files must contain at least one “gene” in the feature type column.

For more detailed information about GFF formats refer to this, this or this link.

12.5 Counts file

• File type and extension: tab-separated plain text file (.counts extension).
• Source: Generated by SeqNjoy in the “Quantify Reads” step of the analysis.
  
• Use: count files are used as input files in the “Compare samples” step of the analysis. When comparing two groups of samples in two different conditions, the counts files of the samples in “Case condition” are compared with the counts files of the samples in “Control condition”.
• Content: The format of a counts file is very simple. It has only two columns: the Feature_id and the number of raw counts for that feature id. The feature type of the ids in the first column corresponds to that selected for quantifying in the “Quantify Genomic Features” dialog, for example, gene ids if “gene” was selected as “feature”. The second column, raw counts, contains the number of reads that map in the reference genome within the feature span.

Figure 12.6: Counts file example

12.7 Differential Expression file

• File type and extension: Differential expression files (diffexp) are Excel files with extension .xlsx.
• Source: Generated by SeqNjoy in the Compare samples step of the analysis.
  
• Use: Differential expression (diffexp) file is the final result of the analysis in SeqNjoy. Additionally, you can carry out further graphical analysis of the results contained in diffexp files by accessing the FIESTA 2.0 utility from the Project’s Main Page.
• Content: diffexp files contain information generated by the program selected for the detection of Differential Expressed Genes (DEGs) between two groups of samples that you determine to be in different conditions (Figure 12.8). In addition, diffexp files contain the raw counts and the normalized counts of all the samples involved in the comparison and the location of the genes in the genome (information taken from the GFF file). Finally, if an annotations file (ANN) is provided in the Calculate Differential Expression dialog, the diffexp file incorporates the biological information contained in the ANN file for each matching gene of the diffexp file.

Figure 12.8: diffexp file format.

This is the meaning of the columns present in the diffexp files (Figure 12.8):

• ID: Feature ID. Identifier of the feature type that was counted. In the case of “gene” this column contains gene IDs.
  • Type: Feature Type. e.g. gene.
  • Location: Chromosome, initial position and end position of the feature in the genome as annotated in the GFF file, in the format Chromosome:IniPos-EndPos.
  • Strand: Strand of the feature (as stated in the GFF file).
  • Log2BaseMean: Average of the normalized counts taken over all samples in logarithm to base 2. This value is calculated by the program that performs the differential expression (DESeq2 or edgeR).
  • FoldChange(CaseCondition/ControlCondition): The fold change or “ratio” is a measure that describes how much the expression of the feature changes from the Control to the Case condition. For example, given a particular gene with a expression value (averaged between the biological replica of this condition and normalized) for Control equal to 30.5 and for Case is equal to 61.0, the ratio would be 2. That is, “Case” expression is two-fold the “Control” expression for that gene.
  • Log2Ratio(CaseCondition/ControlCondition): Logarithm to base 2 of the ratio (fold change) between the counts for Case condition and the counts for the Control condition. Also known as “log Fold Change”. This value is calculated by the program that performs the differential expression (DESeq2 or edgeR). The value in this column is the logarithm to base 2 of the ratio, which is used commonly for analysis and visualization purposes as it is easy to interpret, e.g. a doubling value of the ratio (ratio=2) is a logRatio of 1; if it is a quadrupling value (ratio=4) then the logRatio is 2 and so on. Moreover, the logRatio is symmetric: when a change decreases by an equivalent amount e.g. a reduction by half (ratio=0.5) is equal to a logRatio=-1, a reduction by quarter is equal to a logRatio=-2, etc.
  • StdErr(Log2Ratio): The standard error of the Log2Ratio estimate. Only in case you have selected DESeq2 for differential expression detection.
  • P-value: Statistical significance (p-value). The p-value is a measure of how likely you are to get the particular LogRatio of that gene (or a more extreme value) if no real difference existed. Therefore, a small p-value indicates that there is a small chance of getting this LogRatio if no real difference existed, so logically you decide that the difference of expression levels between both conditions is significant. By “small” often a threshold of p-value=0.05 is taken, but note that it is a convention.
  • P-adjusted: This field contains an “adjusted p-value” obtained with the “False Discovery Rate” (FDR) method by Benjamini and Hochberg (Benjamini & Hochberg, 1995). FDR is a method to cope with the “multiple comparisons problem”. In this link you will find a very good explanation on this problem and why it is a good idea to take the FDR as “the” measure for evaluating statistical significance of your results.
  • Annotations columns: Additional columns will be displayed only if you provided an annotations file in the “Calculate Differential Expression” dialog. The annotations contained for each gene/feature in the annotations file will be added to the differential expression information of that gene/feature.
  • Raw counts for the CONTROL group of samples (or biological replicas): Several columns. Each column contains the raw counts detected in the Quantify Genomic Features menu task for each gene/feature in a “Control” sample (or biological replica), i.e. the counts contained in the corresponding counts files “as-is”.
  • Raw counts for the CASE group of samples (or biological replicas): Idem as the previous “raw counts” but for “Case” samples (or biological replicas).
  • Norm counts columns for CONTROL group of samples (or biological replicas): Several columns. Each column contains the normalized counts for each gene/feature in each sample (or biological replica). The normalization method depends on the program selected to perform the Compare samples step. In the normalization step of the differential expression analysis, raw counts are adjusted to account for factors that prevent direct comparison of expression measures. Considered factors in the normalization calculations are the sequencing depth and the RNA composition. See a very intuitive explanation at this link.
  • Norm counts columns for CASE group of samples (or biological replicas): Idem as the previous “norm counts” but for “Case” samples (or biological replicas).