Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
Loading...
This tutorial gives an overview of RNA-Seq analysis with Partek Flow. It will guide you through creating an RNA-Seq analysis pipeline. The goals of the analysis are to create a list of differentially expressed genes, visualize these gene expression signatures by hierarchical clustering, and interpret the gene lists using gene ontology (GO) enrichment.
This tutorial will illustrate:
This tutorial uses a subset of the data set published in Xu et al. 2013 (PMID: 23902433). In the experiment, mRNA was isolated from HT29 colon cancer cells treated with the drug 5-aza-deoxy-cytidine (5-aza) at three different doses: 0μM (control), 5μM, or 10μM. The mRNA was sequenced using Illumina HiSeq (paired end reads). The goal of the experiment was to identify differentially expressed genes between the different treatment groups.
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
With our reads trimmed, we now have high-quality reads for each sample. The next step is to align the reads to a reference genome. Alignment matches each of the short sequencing reads to a location in the reference genome.
Click Trimmed reads
Click Aligners in the task menu to display available aligners (Figure 1)
Partek Flow offers a variety of different aligners. Mouse over any option for a short description. For this tutorial, we will use STAR, a fast and accurate aligner commonly used for RNA-Seq data. For more information about STAR and the other aligners, please consult the Aligners user guide.
Click STAR
The STAR aligner options allow us to select the genome build (assembly) and index. For this tutorial, our data set contains only reads that map to chromosome 22 to minimize the time required for resource-intensive tasks, such as alignment.
Click Finish to run with hg19 selected for Assembly and Whole genome for the Aligner index (Figure 2)
Alignment is a resource-intensive task and may take around 20 minutes to complete, even when mapping only reads from a single chromosome. Task and data nodes that have been queued, but not completed, are shown in a ligher color than completed tasks (Figure 3).
The Align reads task generates an Aligned reads data node once complete. You can wait for the alignment task to finish or you can continue building the pipeline while the results of alignment are pending; additional tasks can be added to the pipeline and queued before the current task has completed.
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
With attributes added, we can begin building our pipeline.
Click the Analyses tab
In the Analysis tab, data are represented as circles, termed data nodes. One data node, mRNA, should be visible in the Analysis tab (Figure 1).
Click the mRNA node
Clicking a data node brings up the context-sensitive task menu with tasks that can be performed on the data node (Figure 2).
Pre-alignment QA/QC assesses the quality of the unaligned reads and will help us determine whether trimming or filtering is necessary.
Click Pre-alignment QA/QC in the QA/QC section of the task menu
Click Finish to run the task with default settings
Running a task creates a task node, e.g. the blue rectangle labeled Pre-alignment QA/QC (Figure 3), which contains details on the task and a report.. While tasks have been queued or are in progress they have a lighter color. Any output nodes that the task will generate are also displayed in a lighter color until the task completes. Once the task begins running, a progress bar is displayed on the task node.
Click the Pre-alignment QA/QC node
The context-sensitive task menu (Figure 4) shows the option to view the Task report and the Task details. You can also access a task report by double-clicking on a task node.
Click Task report
Pre-aligment QA/QC provides information about the sequencing quality of unaligned reads (Figure 5). Both project level summaries and sample-level summaries are provided.
Click sample SSR592573 in the data table of the report to open its sample-level report
The Average base quality score per position graph in the upper right-hand panel (Figure 6) gives the average Phred score for each position in the reads.
A Phred score is a measure of base call accuracy with a higher score indicating greater accuracy.
By convention, a score above 20 is considered adequate. As you can see, the standard error bars in the graph show that some reads have quality scores below 20 for some of their base pair calls near the 3' end.
Based on the results of Pre-alignment QA/QC, while most of the reads are high quality, we will need to perform read trimming and filtering. For more information about the information included in the task report, please see the Pre-alignment QA/QC user guide.
Click RNA-Seq 5-AZA to return to the Analyses tab
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
Attributes describe samples. Examples of sample attributes include treatment group, age, sex, and time point. Attributes can be added individually in the Metadata tab or in bulk using a text file. In this tutorial, we will add one attribute, 5-AZA Dose, manually.
Click the Metadata tab
Click Manage under Sample attributes (Figure 1)
Click Add new attribute (Figure 2)
To configure a new attribute, at Name, type in 5-AZA Dose as the name of the attribute
Click Add to add 5-AZA Dose as a categorical, project-specific attribute (Figure 3)
Name the first New category 0uM
Click the green plus icon to add category (Figure 4)
Repeat for two additional categories, 5uM and 10uM (Figure 5)
Click Back to metadata tab
The data table now includes an Attribute column for 5-AZA Dose (Figure 6). Next, we need to assign samples attribute categories for 5-AZA Dose.
Select Assign values
The option to edit the 5-AZA Dose field for each sample will appear as a drop-down menu (Figure 7).
Select the 5-AZA Dose text box for a sample to bring up a drop-down menu with the 5-AZA Dose attribute categories (0uM, 5uM, 10uM)
Use the drop-down menus to add a treatment group for each sample
The first three samples (SRR592573-5) should be 0uM, the next three samples (SRR592576-8) should be 5uM, and the final three samples (SRR592579-81) should be 10uM (Figure 8).
Click Apply changes
The data table will now show each a 5-AZA Dose attribute for each sample.
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
Based on pre-alignment QA/QC, we need to trim low quality bases from the 3' end of reads.
Click the Unaligned reads data node
Click Pre-alignment tools in the task menu
Click Trim bases (Figure 1)
By default, Trim bases removes bases starting at the 3' end and continuing until it finds a base pair call with a Phred score of equal to or greater than 35 (Figure 2).
Click Finish to run Trim bases with default settings
The Trim bases task will generate a new data node, Trimmed reads (Figure 3). We can view the task report for Trim bases by double-clicking either the Trim bases task node or the Trimmed reads data node or choosing Task report from the task menu.
Double-click the Trimmed reads data node to open the task report
The report shows the percentage of trimmed reads and reads removed in a spreadsheet and a two graphs (Figure 4).
The results are fairly consistent across samples with ~2% of reads untrimmed, ~86% trimmed, and ~12% removed for each. The average quality score for each sample is increased with higher average quality scores at the 3' ends.
Click RNA-Seq 5-AZA to return to the Analyses tab
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
| Phred Quality Score | Probability of incorrect base call | Base call accuracy |
|---|---|---|
20
1 in 100
99%
30
1 in 1000
99.9%
40
1 in 10,000
9.99%
50
1 in 100,000
99.999%
After normalizing the data, we can perform differential analysis to identify genes that are differentially expressed based on treatment.
Click the Normalized counts node
Click Statistics in the task menu
Click Differential analysis in the task menu (Figure 1)
Check 5-AZA Dose and click Add factors to add the attribute to the statistical model.
Select Next to continue with 5-AZA Dose as the selected attribute
The Comparisons page will open (Figure 4).
It is easiest to think about comparisons as the questions we are asking. In this case, we want to know what are the differentially expressed genes between untreated and treated cells. We can ask this for each dose individually and for both collectively.
The upper box will be the numerator and the lower box will be the denominator in the comparison calculation so we will select the 0μM control in the lower box.
Drag 5μM to the upper box
Drag 0μM to the lower box
Click Add comparison to add 5μM vs. 0μM to the comparison table (Figure 5)
Repeat to create comparisons for 10μM vs. 0μM and 10μM,5μM vs. 0μM (Figure 6)
Click Finish to perform DESeq2 as configured
A DESeq2 task node and a DESeq2 data node will be added to the pipeline (Figure 7).
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
The tutorial data set includes 9 samples equally divided into 3 treatment groups. Sequencing was performed by an Illumina HiSeq (paired-end reads), but the workflow can be easily adapted for data generated by other sequencers. Each sample has 2 fastq files for a total of 18 fastq files.
You can obtain the tutorial data set through Partek Flow.
Click your avatar
Click Settings in the drop-down menu (Figure 1)
At the top of the System information page, there is a section labeled Download tutorial data (Figure 2).
Click RNA-Seq 5-AZA to download the tutorial data set
A new project will be created and you will be directed to the Analyses Tab. The data will be downloaded automatically (Figure 3) and imported into your project. Because this is a tutorial project, there is no need to click on Add data as it will be done automatically.
At first the project is empty, but the file download will start automatically in the background. You can wait a few minutes then refresh your browser or you can monitor the download progress using the Queue.
Click Queue
Click View Queued Tasks in the drop-down menu
The Queued tasks page will open (Figure 4).
Click Projects
Click RNA-Seq 5-AZA in the drop-down menu
The Analyses tab will open (Figure 5). If you download has completed, you will see a blue circle titled mRNA.
Once the download completes, the sample table will appear in the Metadata tab.
Click the Metadata tab The Metadata tab includes the sample table with the names of each imported sample (Figure 6).
In the next section of the tutorial, we will add a sample attribute that indicates the treatment group of each sample.
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
Low expression genes may be indistinguishable from noise and will decrease the sensitivity of differential expression analysis.
Click the Gene counts node
Click Filtering in the task menu
Click Filter features (Figure 1)
Click Noise reduction filter
Set the filter to maximum <= 10
Click Finish (Figure 2)
A new Filtered counts node will be created (Figure 3).
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
After alignment has completed, we can view the quality of alignment by performing post-alignment QA/QC.
Click the Aligned reads data node
Click QA/QC inthe task menu
Click Post-alignment QA/QC from the QA/QC section of the task menu (Figure 1)
A Post-alignment QA/QC task node will be generated (Figure 2).
Double-click the Post-alignment QA/QC task node to view the task report
Similar to the Pre-alignment QA/QC task report, general quality information about the whole data set is displayed and sample-level reports can be opened by clicking a sample name in the table.
The top two graphs in the data set view (Figure 3) show the alignment breakdown and coverage.
From these graphs, we can see that more than 95% of reads were aligned, but the total number of reads for each sample varies. Normalizing for the variability in total read counts will be addressed in a later section of the tutorial.
For more information about the graphs and information presented in the Post-alignment QA/QC task report, see the Post-alignment QA/QC user guide.
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
Because different samples have different total numbers of reads, it would be misleading to calculate differential expression by comparing read count numbers for genes across samples without normalizing for the total number of reads.
Click the Filtered counts data node
Click Normalization and scaling in the task menu
Click Normalization (Figure 1)
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
The Count normalization menu will open (Figure 2).
Normalization can be performed by sample or by feature. By sample is selected by default; this is appropriate for the tutorial data set.
Available normalization methods are listed in the left-hand panel. For more information about these options, please see the Normalize counts user guide.
For this tutorial, we will use the recommended default normalization settings.
This adds the Median ratio normalization method, which is suitable for performing differential expression analysis using DESeq2 (Figure 3).
Click Finish to perform normalization
A Normalize counts task node and a Normalized counts data node are added to the pipeline (Figure 4)
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
RNA-Seq uses the number of sequencing reads per gene or transcript to quantify gene expression. Once reads are aligned to a reference genome, we need to assign each read to a known transcript or gene to give a read-count per transcript or gene.
Click the Aligned reads data node
Click Quantification in the task menu
We will use Partek E/M to quantify reads to an annotation model in this tutorial. For more information about the other quantification options, please see the Quantification user guide.
Click Quantify to an annotation model (Partek E/M) (Figure 1)
Choose the latest RefSeq Transcripts 95 annotation from the Gene/feature annotation drop-down menu (you may need to download it first, via Library File Management)
Click Finish (Figure 2)
The Quantify to annotation model task node outputs two data nodes, Gene counts and Transcript counts (Figure 3).
To view the results of quantification, we can select either data node output.
Double-click the Gene counts data node to view the task report
The task report details the number of reads within exons, introns, and intergenic regions. For detailed information about the quantification results, see the Quantify to annotation model (Partek E/M) user guide.
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
The principal components analysis (PCA) scatter plot allows us to visualize similarities and differences between the samples in a data set.
Click the Normalized counts data node
Click Exploratory analysis in the task menu
Click PCA
Click Finish to run PCA with the default options
The PCA task node will be added to the pipeline (Figure 1)
Double click the PCA data node to open the PCA scatter plot (Figure 2)
In the Data Viewer, click Style under Configure and set the Color by drop-down to 5-AZA Dose. The scatter plot shows each sample as a sphere, colored by treatment group, in a three dimensional plot. The x, y, and z axes are the first three principal components. The percentage of total variance explained by each is listed next to the axis label. The size of each axis is determined by the variance along that axis. The plot is fully interactive; it can be rotated and points selected.
Here, we can see that samples separate based on treatment, but there is noticeable separation within treatment groups, particularly the 0μM and 10μM treatment groups.
For more detailed information about the PCA scatter plot, please see the PCA user guide.
If you need additional assistance, please visit our support page to submit a help ticket or find phone numbers for regional support.
To check how well our list of differentially expressed genes distinguishes one treatment group from another, we can perform hierarchical clustering based on the gene list. Clustering can also be used to discover novel groups within your data set, identify gene expression signatures distinguishing groups of samples, and to identify genes with similar patterns of gene expression across samples.
Click the Feature list data node
Click Exploratory analysis in the task menu
Click Hierarchical clustering (Figure 1)
The Hierarchical clustering menu will open (Figure 2). Hierarchical clustering can be performed with a heatmap or bubble map plot. Cluster must be selected under Ordering for both Feature order and Cell order if both the features (columns) and cells (rows) are to be clustered.
Click Finish to run with default settings
A Hierarchical clustering task node will be added to the pipeline (Figure 3).
Double-click the Hierarchical clustering / heatmap task node to launch the heatmap
The Dendrogram view will open showing a heatmap with the hierarchical clustering results (Figure 4).
Samples are shown on rows and genes on columns. Clustering for samples and genes is shown through the dendrogram trees. More similar samples/genes are separated by fewer branch points of the dendrogram tree.
The heatmap displays standardized expression values with a mean of zero and standard deviation of one.
The heatmap can be customized to improve data visualization using the menu on the Configuration panel on the left.
Expand the Annotations > Data card.
In the dialog, click on the Gene counts node
Now set the Row annot to 5-AZA Dose
Samples are now labeled with their 5-AZA Dose group (Figure 5).
Samples cluster based on treatment group and the 5μM and 10μM groups are more similar to each other than to the 0μM group.
We can save the heatmap as a publication-quality image.
Choose size and resolution using the Save as SVG dialog (Figure 6)
Select Save
The heatmap will be saved as a .svg file and downloaded in your web browser.
Once we have performed DESeq2 to identify differentially expressed genes, we can create a list of significantly differentially expressed genes using cutoff thresholds.
Double click the Feature list data node to open the task report
The task report spreadsheet will open showing genes on rows and the results of the DESeq2 on columns (Figure 1).
To get a sense of what filtering thresholds to set, we can view a volcano plot for a comparison.
A volcano plot will open showing p-value on the y-axis and fold-change on the x-axis (Figure 2). If the gene labels are on (not shown), click on the plot to turn them off.
Thresholds for the cutoff lines are set using the Statistics card (Configuration panel > Configure > Statistics). The default thresholds are |2| for the X axis and 0.05 for the Y axis.
Switch to the browser tab showing the DESeq2 report
Click FDR step up
Click the triangle next to FDR step up to open the FDR step up options
Leave All contrasts selected
Set the cutoff value to 0.05. Hit Enter.
This will include genes that have a FDR step up value of less than or equal to 0.05 for all three contrasts, 5μM vs. 0μM, 10μM vs. 0μM and 5μM:10μM vs. 0μM. FDR step up is the false discovery rate adjusted p-value used by convention in microarray and next generation sequencing data sets in place of unadjusted p-value.
Click Fold-change
Click the triangle next to Fold-change to open the Fold-change options
Leave All contrasts selected
Set to From -2 to 2 with Exclude range selected. Hit Enter.
Note that the number of genes that pass the filter is listed at the top of the filter menu next to Results: and will update to reflect any changes to the filter. Here, 27 genes pass the filter (Figure 3). Depending on your settings, the number may be slightly different.
This creates a Filter list task node and a Filtered feature list data node (Figure 4).
To learn more about the biology underlying gene expression changes, we can use gene ontology (GO) or pathway enrichment analysis. Enrichment analysis identifies over-represented GO terms or pathways in a filtered list of genes.
Click the filtered Filtered feature list data node
Click Biological interpretation in the task menu
Click Gene set enrichment then select Gene set database to perform GO enrichment analysis (Figure 1)
Select the latest gene set from geneontology.org from the Gene set database drop-down menu
Click Finish
A GO e_nrichment_ task node will be added to the pipeline (Figure 2).
Double-click the Gene set enrichment task node to open the task report (Figure 3)
The GO e_nrichment_ task report spreadsheet lists GO terms by ascending p-value with the most significant GO term at the top of the list. Also included are the enrichment score, the number of genes from that GO term in the list, and the number of genes from that GO term that are not in the list.
KEGG enrichment analysis identifies pathways that are over-represented in a gene list data node.
Click the filtered Filtered feature list data node
Click Biological interpretation in the task menu
Click Gene set enrichment then select KEGG database
Click Finish in the configuration dialog to run KEGG analysis with the Homo sapiens KEGG database
A Pathway e_nrichment_ task node will be added to the pipeline (Figure 4).
Double-click the Pathway enrichment task node to open the task report
The Pathway enrichment task report is similar to the Enrichment analysis task report (Figure 5).
To view an interactive KEGG pathway map, click the pathway ID (Gene set column).
By following the steps in this tutorial, you have built a pipeline. You can save this pipeline for future use.
Select Create new pipeline near the bottom left-hand side of the browser window
Select the Pre-alignment QA/QC, Tr****im bases, Align reads, Post-alignment QA/QC, Quantify to annotation model, Filter features, Normalize counts, PCA, and GSA task nodes to include them in the pipeline
Name the pipeline; we have chosen RNA-Seq basic analysis
Give a description for the pipeline; we have noted trim <20, STAR, normalize with total count and add 0.001, GSA
Select Create pipeline (Figure 1)
To access this pipeline in the future, select an unaligned reads data node and choose Pipelines from the task menu. Available saved pipelines will be available to choose from the Pipelines section of the task menu (Figure 2).
In addition to the volcano plot showing all genes, we can view expression levels of each gene on a dot plot.
Double-click the **Filtered feature list **data node to open the task report
Click the FDR step up header in the 5uM vs. 0uM section to sort by ascending FDR step up
In the task report table, there is a column labeled View with three icons in each row.
Select to open a dot plot for the gene SELENOM
The dot plot for SELENOM (Figure 1) shows each sample as a point with normalized reads on the y-axis. Samples are separated and colored by treatment group.
From the DESeq2 task report, we can browse to any gene in the Chromosome view.
Click in the SELENOM row to open Chromosome view (Figure 1)
A new tab will open showing SELENOM in the Chromosome view (Figure 2).
Chromosome View shows reference genome, annotation, and data set information together aligned at genomic coordinates.
The top track shows average number of total count normalized reads for each of the three treatment groups in a stacked histogram. The second track shows the RefSeq annotation.
We can add tracks from any data node using Select Tracks.
Click Select tracks
A pop-up dialog showing the pipeline allows us to choose which data to display as tracks in Chromosome view (Figure 3).
Click Reads pileup under Aligned reads on the left-hand side of the dialog
Click Display tracks to make the change
The reads pileup track is now included (Figure 4).
Select the appropriate differential analysis method (Figure 2). In this tutorial we are going to use DESeq2, but Partek Flow offers a number of alternatives. Hover the mouse over the symbol for more information on each differential analysis method, or see our Differential Analysis user guide for a more in-depth look.
Select
We will use the default options for quantification. To learn more about the different options, please see the Quantify to annotation model (Partek E/M) user guide or mouse over the next to each option.
Click on the gray button on the right hand side of Row annot None available
Click the Export image icon in the top right corner of the plot
For more information about hierarchical clustering and the Dendrogram view, please see the user guide.
If you need additional assistance, please visit to submit a help ticket or find phone numbers for regional support.
Click next to the 5uM vs. 0uM comparison
Click to create a data node with only the genes that pass the filter
If you need additional assistance, please visit to submit a help ticket or find phone numbers for regional support.
To view the genes associated with each GO term, select to open the extra details page. To view additional information about a GO term, click the blue gene set ID to open the linked geneontology.org entry in a new tab.
For more information about GO enrichment analysis, please see the user guide.
If you need additional assistance, please visit to submit a help ticket or find phone numbers for regional support.
After selecting the pipeline, you will be prompted to choose the reference genome for alignment, the annotation for quantification, and the contrasts for GSA. After selecting these options, the pipeline will automatically run. See for more information.
If you need additional assistance, please visit to submit a help ticket or find phone numbers for regional support.
For more information about the dot plot, please see the user guide. To return to the DESeq2 report, switch to the browser table with filtered feature list.
In the DESeq2 report, you can select to view additional information about the statistical results for a gene or select to view the region in Chromosome View. Chromosome View is discussed in the next section of the tutorial.
If you need additional assistance, please visit to submit a help ticket or find phone numbers for regional support.
Each track has Configure track and Move track buttons that can be used to modify each track.
To learn more about Chromosome view, please consult the user guide.
If you need additional assistance, please visit to submit a help ticket or find phone numbers for regional support.
