<commit_before>---
layout: post
title: "Introduction to RNA-Seq - updated"
modified:
categories: courses
excerpt:
tags: [ngs, rna, introduction, galaxy, qc, alignment, de]
image:
  feature:
date: 2015-10-22
---

<section id="table-of-contents" class="toc">
  <header>
    <h3>Overview</h3>
  </header>
<div id="drawer" markdown="1">
*  Auto generated table of contents
{:toc}
</div>
</section><!-- /#table-of-contents -->

We will be retracing most of the steps required to get from an Illumina FASTQ sequence file received from a sequencing facility as part of an RNA sequencing analysis all the way to transcript assembly and identifying  differentially expressed genes in two different ways. We will focus on the following  steps: initial quality control, read mapping, transcript isoform determination (reference-guided) followed by quantification. This session borrows material from [Jeremy's RNA-seq sample history](https://main.g2.bx.psu.edu/u/jeremy/p/galaxy-rna-seq-analysis-exercise), a [Galaxy workflow](https://usegalaxy.org/u/mejia-guerra/w/basic-rna-seq-analysis---differential-expression-functional-genomics-workshop-2012) and the [UC Davis Bioinformatics Core workshop](http://training.bioinformatics.ucdavis.edu/docs/2012/05/RNA/index.html). 


#### Access to Galaxy
{:.no_toc}

Before tackling this module you should have a [working knowledge of Galaxy](http://bioinformatics.sph.harvard.edu/ngs-workshops/courses/introduction-to-galaxy/), understand how to find relevant tools, annotate your history or get new data sets from the data library.  As a result we will minimize the screenshots and pointers in this document, focusing just on _what_ to do while guiding you through the process in class. As always, do ask questions and most importantly collaborate with your fellow students.

In order for you to be able to access Galaxy on your assigned dedicated machine on the Cloud, you have been given a web or IP address in the form of A.B.C.D where A, B, C and D are numbers separated by dots. You will need it in order to access Galaxy from the web browser on your laptop. 


#### Changes to the default Galaxy installation
{:.no_toc}

The default Galaxy instance relies on [TopHat/Cufflinks](http://tophat.cbcb.umd.edu/manual.html) for RNA-seq analysis, part of a well-tested workflow system. If you are not interested in discovering splice variants -- which can be prone to false positives -- a simple alignment with TopHat followed by analysis of differential expression with [DESeq](http://bioconductor.org/packages/release/bioc/html/DESeq.html) / [DSS](http://www.bioconductor.org/packages/release/bioc/html/DSS.html) (genes) or [DEXSeq](http://bioconductor.org/packages/release/bioc/html/DEXSeq.html) (exons) is likely to be simpler. For this workshop we added a more recent version ([TopHat2](http://genomebiology.com/2013/14/4/r36)) from the [Galaxy ToolShed](http://toolshed.g2.bx.psu.edu/) to our Galaxy instance, and will be running [edgeR](http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) for the differential expression analysis.

> If you want to work through this session outside of this course, take a look at our section on [how to set up your own Galaxy server](http://bioinformatics.sph.harvard.edu/ngs-workshops/posts/running-your-own-galaxy-instance/) or download the [complete introduction dataset](https://dl.dropbox.com/u/407047/Blog/Data/RNA-Seq_Data_files.tgz) and use one of the public Galaxy instances.

## Quality control

The dataset we are using is part of a larger study described in [<NAME> et al, Cell Rep 2014](http://www.ncbi.nlm.nih.gov/pubmed/25464849). The authors are investigating interactions between various genes involved in Fragile X syndrome, a disease in which there is aberrant production of the FMRP protein. FMRP has been linked to the microRNA pathway, as it has been shown to be involved in miRNA mediated translational suppresion. **The authors sought to show that FMRP associates with the RNA helicase MOV10 for regulation of brain mRNAs.**

From this study we are using the [RNA-Seq](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50499) data which is publicly available in the [SRA](http://www.ncbi.nlm.nih.gov/sra). The RNA was extracted from HEK293F cells that were transfected with a MOV10 transgene and normal control cells. Using this data, we will evaluate transcriptional patterns associated with MOV10 overexpression. 

The libraries for this dataset are stranded and were generated using the dUTP method. Sequencing was carried out on the Illumina HiSeq-2500 for 100bp single end reads. The full dataset was sequenced to ~40 million reads per sample, but for this workshop we will be looking at a small subset on chr1 (~300,000 reads/sample) to keep things manageable and allow algorithms to finish within a few minutes. We will switch over to using all samples for the differential expression analysis by retrieving them from the `Data Libraries` in the `Shared Data` menu as needed.For each group we have three replicates as described in the figure below.


![Experimental_design](../../images/screenshots/exp_design.png)


### Exploring the FASTQ files

Open up Galaxy from your machine as before, start up a new blank history and give it a name. In the top menu, move your mouse cursor over `Shared data` and select `Data libraries`. This allows you to import files into Galaxy that have been shared with you.

Browse to the `RNA-Seq` entry, open the`Sequence and reference data` folder and mouse over the arrow to the right of the filename `Mov10oe1.fastq`. In the context menu select `Import this dataset into selected histories`. Move back to the `Analyze Data` menu.

This is the kind of data you would expect to receive from your sequencing core facility: a large number of short nucleotide reads in [FASTQ format](https://dl.dropbox.com/u/407047/Blog/Documents/Literature/QC/Nucleic%20Acids%20Res%202009%20Cock.pdf) (PDF). 

> Take a look at the file contents using Galaxy's preview and file view functionality. What additional information does FASTQ contain over FASTA? How many reads are in your file? Pick the first sequence in your FASTQ file and note the identifier. What is the quality encoding character of it's first and last nucleotide, respectively?

## Quality Controls

The FASTQ file contains output reads from the sequencer that need to be mapped to a reference genome for us to understand where those reads came from on the sequenced genome. However, before we can delve into read mapping, we first need to make sure that our preliminary data is of sufficiently high quality. This involves several steps:

1. Obtaining summary quality statistics on the reads and review diagnostic graphs 
2. Eliminate sequencing artifacts
3. Filter out genetic contaminants (primers, vectors, adaptors)
4. Filter out low-quality reads
5. Recalculate quality statistics and review diagnostic plots on filtered data

Iterate through steps 2-5 until the data is of sufficient quality before proceeding to mapping.

### Obtain Quality Statistics

Under the `NGS:QC and manipulation` tool heading, select the `Compute quality statistics` tool (from the `FASTX toolkit` section) and apply it to your FASTQ file.

> Explore the results. What information does the data matrix contain?

Galaxy contains a number of tools to convert this matrix into graphs (e.g., `Draw Quality Score Boxplot` and `Draw Nucleotides Distribution Chart`), but in most cases you will stick to using a standard QC tool such as [FASTQC](http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/ ) with provides many more diagnostic plots. Find `FASTQC` in the Galaxy tool kit and give it a try.

> Do you observe similar quality statistics as before? What additional information do you get with this tool? What can you say about the per-base GC content? Does this match the human genome? Why is there a shift in the GC distribution from the theoretical distribution? Finally, what is the meaning of the k-mer content tab? This publication from [Schroeder et al](https://dl.dropbox.com/u/407047/Blog/Documents/Literature/QC/PLoS%20ONE%202010%20Schr%C3%B6der.pdf) (PDF) has some background information.

Compare your results with a [pre-generated FastQC plot](https://dl.dropboxusercontent.com/u/74036176/Mov10oe_1-fastqc_report.html) obtained from the full FASTQ data set (i.e., using all reads). 

> Do you note any differences in the output or format? Why is that?

Some of the diagnostic plots from a [Core conference call](http://bioinfo-core.org/index.php/9th_Discussion-28_October_2010) might be useful if you are curious to see what different modes of failure frequently look like.

Depending on your FASTQ results you may need to test for contamination (either by vector or other sourcers) and/or remove adapter sequences from your read prior to alignment. It is good to remove these sequences to increase your mapping efficiency. 

### Screen for contamination

Not all required tools are available through Galaxy just yet. For example, depending on the source of your data you probably want to screen your reads for vector or adapter contamination using tools such as [SeqTrim](https://dl.dropbox.com/u/407047/Blog/Documents/Literature/QC/BMC%20Bioinformatics%202010%20Falgueras.pdf) (PDF). Other tools check for sample cross-contamination ([ContEst](https://dl.dropbox.com/u/407047/Blog/Documents/Literature/QC/Bioinformatics%202011%20Cibulskis.pdf) (PDF)) or contamination with data from other species ([FASTQ Screen](http://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/)).  For now we will just explore how you can spot potential contamination from the FastQC results alone. 

> Where is the most overrepresented sequence in your FastQC results from the h1-hESC sample coming from?

Compare your own FastQC results with a [pre-generated FastQC summary](http://dl.dropbox.com/u/4253254/Training/HSCI_Workshop/D19VHACXX_s1_0_illumina12index_12_SL17357_fastqc/fastqc_report.html) that highlights a potential problem.

> Study the report. Is there any vector contamination present? Is there any adapter contamination present?

As an emergency solution to remove contaminating sequences, Galaxy comes with a way to `Clip adapter sequences`. If you notice a particularly strong contamination with a vector or primer you can use your FASTQ file as the `Library to clip`, keep the default `Minimum sequence length`, and enter the sequence you spotted in the FASTQC report as the `Adapter` to clip. You probably want to keep both clipped and unclipped sequences in this case; keeping only reads that contained an adapter or barcode can be useful if you expect this to be a feature of all valid reads. For command-line use tools such as [cutadapt](https://code.google.com/p/cutadapt/) or [AlienTrimmer](ftp://ftp.pasteur.fr/pub/gensoft/projects/AlienTrimmer/) are more versatile and less labor-intensive to use.
 
Normally, you would take the filtered output and feed it back into Galaxy, but our dataset has very little contamination so we will move on with the original data set. 
<commit_msg>typo correction<commit_after>---
layout: post
title: "Introduction to RNA-Seq - updated"
modified:
categories: courses
excerpt:
tags: [ngs, rna, introduction, galaxy, qc, alignment, de]
image:
  feature:
date: 2015-10-22
---

<section id="table-of-contents" class="toc">
  <header>
    <h3>Overview</h3>
  </header>
<div id="drawer" markdown="1">
*  Auto generated table of contents
{:toc}
</div>
</section><!-- /#table-of-contents -->

We will be retracing most of the steps required to get from an Illumina FASTQ sequence file received from a sequencing facility as part of an RNA sequencing analysis all the way to transcript assembly and identifying  differentially expressed genes in two different ways. We will focus on the following  steps: initial quality control, read mapping, transcript isoform determination (reference-guided) followed by quantification. This session borrows material from [Jeremy's RNA-seq sample history](https://main.g2.bx.psu.edu/u/jeremy/p/galaxy-rna-seq-analysis-exercise), a [Galaxy workflow](https://usegalaxy.org/u/mejia-guerra/w/basic-rna-seq-analysis---differential-expression-functional-genomics-workshop-2012) and the [UC Davis Bioinformatics Core workshop](http://training.bioinformatics.ucdavis.edu/docs/2012/05/RNA/index.html). 


#### Access to Galaxy
{:.no_toc}

Before tackling this module you should have a [working knowledge of Galaxy](http://bioinformatics.sph.harvard.edu/ngs-workshops/courses/introduction-to-galaxy/), understand how to find relevant tools, annotate your history or get new data sets from the data library.  As a result we will minimize the screenshots and pointers in this document, focusing just on _what_ to do while guiding you through the process in class. As always, do ask questions and most importantly collaborate with your fellow students.

In order for you to be able to access Galaxy on your assigned dedicated machine on the Cloud, you have been given a web or IP address in the form of A.B.C.D where A, B, C and D are numbers separated by dots. You will need it in order to access Galaxy from the web browser on your laptop. 


#### Changes to the default Galaxy installation
{:.no_toc}

The default Galaxy instance relies on [TopHat/Cufflinks](http://tophat.cbcb.umd.edu/manual.html) for RNA-seq analysis, part of a well-tested workflow system. If you are not interested in discovering splice variants -- which can be prone to false positives -- a simple alignment with TopHat followed by analysis of differential expression with [DESeq](http://bioconductor.org/packages/release/bioc/html/DESeq.html) / [DSS](http://www.bioconductor.org/packages/release/bioc/html/DSS.html) (genes) or [DEXSeq](http://bioconductor.org/packages/release/bioc/html/DEXSeq.html) (exons) is likely to be simpler. For this workshop we added a more recent version ([TopHat2](http://genomebiology.com/2013/14/4/r36)) from the [Galaxy ToolShed](http://toolshed.g2.bx.psu.edu/) to our Galaxy instance, and will be running [edgeR](http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) for the differential expression analysis.

> If you want to work through this session outside of this course, take a look at our section on [how to set up your own Galaxy server](http://bioinformatics.sph.harvard.edu/ngs-workshops/posts/running-your-own-galaxy-instance/) or download the [complete introduction dataset](https://dl.dropbox.com/u/407047/Blog/Data/RNA-Seq_Data_files.tgz) and use one of the public Galaxy instances.

## Quality control

The dataset we are using is part of a larger study described in [Kenny PJ et al, Cell Rep 2014](http://www.ncbi.nlm.nih.gov/pubmed/25464849). The authors are investigating interactions between various genes involved in Fragile X syndrome, a disease in which there is aberrant production of the FMRP protein. FMRP has been linked to the microRNA pathway, as it has been shown to be involved in miRNA mediated translational suppresion. **The authors sought to show that FMRP associates with the RNA helicase MOV10 for regulation of brain mRNAs.**

From this study we are using the [RNA-Seq](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50499) data which is publicly available in the [SRA](http://www.ncbi.nlm.nih.gov/sra). The RNA was extracted from HEK293F cells that were transfected with a MOV10 transgene and normal control cells. Using this data, we will evaluate transcriptional patterns associated with MOV10 overexpression. 

The libraries for this dataset are stranded and were generated using the dUTP method. Sequencing was carried out on the Illumina HiSeq-2500 for 100bp single end reads. The full dataset was sequenced to ~40 million reads per sample, but for this workshop we will be looking at a small subset on chr1 (~300,000 reads/sample) to keep things manageable and allow algorithms to finish within a few minutes. We will switch over to using all samples for the differential expression analysis by retrieving them from the `Data Libraries` in the `Shared Data` menu as needed.For each group we have three replicates as described in the figure below.


![Experimental_design](../../images/screenshots/exp_design.png)


### Exploring the FASTQ files

Open up Galaxy from your machine as before, start up a new blank history and give it a name. In the top menu, move your mouse cursor over `Shared data` and select `Data libraries`. This allows you to import files into Galaxy that have been shared with you.

Browse to the `RNA-Seq` entry, open the`Sequence and reference data` folder and mouse over the arrow to the right of the filename `Mov10oe1.fastq`. In the context menu select `Import this dataset into selected histories`. Move back to the `Analyze Data` menu.

This is the kind of data you would expect to receive from your sequencing core facility: a large number of short nucleotide reads in [FASTQ format](https://dl.dropbox.com/u/407047/Blog/Documents/Literature/QC/Nucleic%20Acids%20Res%202009%20Cock.pdf) (PDF). 

> Take a look at the file contents using Galaxy's preview and file view functionality. What additional information does FASTQ contain over FASTA? How many reads are in your file? Pick the first sequence in your FASTQ file and note the identifier. What is the quality encoding character of it's first and last nucleotide, respectively?

## Quality Controls

The FASTQ file contains output reads from the sequencer that need to be mapped to a reference genome for us to understand where those reads came from on the sequenced genome. However, before we can delve into read mapping, we first need to make sure that our preliminary data is of sufficiently high quality. This involves several steps:

1. Obtaining summary quality statistics on the reads and review diagnostic graphs 
2. Eliminate sequencing artifacts
3. Filter out genetic contaminants (primers, vectors, adaptors)
4. Filter out low-quality reads
5. Recalculate quality statistics and review diagnostic plots on filtered data

Iterate through steps 2-5 until the data is of sufficient quality before proceeding to mapping.

### Obtain Quality Statistics

Under the `NGS:QC and manipulation` tool heading, select the `Compute quality statistics` tool (from the `FASTX toolkit` section) and apply it to your FASTQ file.

> Explore the results. What information does the data matrix contain?

Galaxy contains a number of tools to convert this matrix into graphs (e.g., `Draw Quality Score Boxplot` and `Draw Nucleotides Distribution Chart`), but in most cases you will stick to using a standard QC tool such as [FASTQC](http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/ ) with provides many more diagnostic plots. Find `FASTQC` in the Galaxy tool kit and give it a try.

> Do you observe similar quality statistics as before? What additional information do you get with this tool? What can you say about the per-base GC content? Does this match the human genome? Why is there a shift in the GC distribution from the theoretical distribution? Finally, what is the meaning of the k-mer content tab? This publication from [Schroeder et al](https://dl.dropbox.com/u/407047/Blog/Documents/Literature/QC/PLoS%20ONE%202010%20Schr%C3%B6der.pdf) (PDF) has some background information.

Compare your results with a [pre-generated FastQC plot](https://dl.dropboxusercontent.com/u/74036176/Mov10oe_1-fastqc_report.html) obtained from the full FASTQ data set (i.e., using all reads). 

> Do you note any differences in the output or format? Why is that?

Some of the diagnostic plots from a [Core conference call](http://bioinfo-core.org/index.php/9th_Discussion-28_October_2010) might be useful if you are curious to see what different modes of failure frequently look like.

Depending on your FASTQ results you may need to test for contamination (either by vector or other sourcers) and/or remove adapter sequences from your read prior to alignment. It is good to remove these sequences to increase your mapping efficiency. 

### Screen for contamination

Not all required tools are available through Galaxy just yet. For example, depending on the source of your data you probably want to screen your reads for vector or adapter contamination using tools such as [SeqTrim](https://dl.dropbox.com/u/407047/Blog/Documents/Literature/QC/BMC%20Bioinformatics%202010%20Falgueras.pdf) (PDF). Other tools check for sample cross-contamination ([ContEst](https://dl.dropbox.com/u/407047/Blog/Documents/Literature/QC/Bioinformatics%202011%20Cibulskis.pdf) (PDF)) or contamination with data from other species ([FASTQ Screen](http://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/)).  For now we will just explore how you can spot potential contamination from the FastQC results alone. 

> Where is the most overrepresented sequence in your FastQC results from the Mov10oe1 sample coming from?

Compare your own FastQC results with a [pre-generated FastQC summary](http://dl.dropbox.com/u/4253254/Training/HSCI_Workshop/D19VHACXX_s1_0_illumina12index_12_SL17357_fastqc/fastqc_report.html) that highlights a potential problem.

> Study the report. Is there any vector contamination present? Is there any adapter contamination present?

As an emergency solution to remove contaminating sequences, Galaxy comes with a way to `Clip adapter sequences`. If you notice a particularly strong contamination with a vector or primer you can use your FASTQ file as the `Library to clip`, keep the default `Minimum sequence length`, and enter the sequence you spotted in the FASTQC report as the `Adapter` to clip. You probably want to keep both clipped and unclipped sequences in this case; keeping only reads that contained an adapter or barcode can be useful if you expect this to be a feature of all valid reads. For command-line use tools such as [cutadapt](https://code.google.com/p/cutadapt/) or [AlienTrimmer](ftp://ftp.pasteur.fr/pub/gensoft/projects/AlienTrimmer/) are more versatile and less labor-intensive to use.
 
Normally, you would take the filtered output and feed it back into Galaxy, but our dataset has very little contamination so we will move on with the original data set. 

{'fraction_non_alphanumeric': 0.07106484610878284, 'fraction_numerical': 0.017754945698706684, 'mean_word_length': 4.743271221532091, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 36, 'https://': 18, 'lorem ipsum': 0, 'www.': 14, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '935181', 'n_tokens_mistral': 6390, 'n_tokens_neox': 5723, 'n_words': 2880}