Mapping to Reference Exercise - Revision history

Rf at 14:02, 14 May 2017

2017-05-14T14:02:04Z

Rf at 17:44, 11 May 2017

2017-05-11T17:44:59Z

Rf at 17:44, 11 May 2017

2017-05-11T17:44:13Z

Rf at 17:36, 11 May 2017

2017-05-11T17:36:31Z

Rf at 16:21, 7 May 2017

2017-05-07T16:21:19Z

Rf at 16:18, 7 May 2017

2017-05-07T16:18:55Z

Rf at 14:52, 7 May 2017

2017-05-07T14:52:59Z

Rf at 15:48, 6 May 2017

2017-05-06T15:48:15Z

Rf at 15:40, 6 May 2017

2017-05-06T15:40:43Z

Rf: Created page with "= Aims = Mapping to a reference genome is a vital step to generate counts and do differential gene expression thereafter. For RNA-Seq data it is important to choose an aligne..."

2017-05-06T15:17:37Z

Created page with "= Aims = Mapping to a reference genome is a vital step to generate counts and do differential gene expression thereafter. For RNA-Seq data it is important to choose an aligne..."

New page

= Aims =

Mapping to a reference genome is a vital step to generate counts and do differential gene expression thereafter. For RNA-Seq data it is important to choose an aligner which is splice aware.

In this part you will learn to:
* align RNA-Seq reads to a reference genome
* calculate the mapping rate

You will use the following software:
* TopHat2 v2.0.11: http://ccb.jhu.edu/software/tophat/index.shtml
* Bowtie2 v2.2.0: http://bowtie-bio.sourceforge.net/bowtie2/index.shtml

The data set you'll be using is downloaded from ENA (http://www.ebi.ac.uk/ena/data/view/SRP019027). The reads belong to sample SRR769316. The data set is tailored with respect to the time allocated for the exercise.

Change directory:

cd /home/training/Data/04_Mapping_to_a_reference_genome

Index the reference genome using Bowtie2:

cd Reference
bowtie2-build mm10_chr19-1-20000000.fasta mm10_chr19-1-20000000

Run the alignment using TopHat2:

cd ..
tophat2 -o tophat2 --no-mixed \
--rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina \
-G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 \
/home/training/Data/03_Quality_control_and_data_preprocessing/Read_1.fastq \
/home/training/Data/03_Quality_control_and_data_preprocessing/Read_2.fastq

where:
* --no-mixed: For paired reads, only report read alignments if both reads in a pair can be mapped
* --rg-id: Read group ID
* --rg-sample: Sample ID
* --rg-center: Sequencing Centre name
* --rg-platform: Sequencing platform descriptor
* -G: Supply TopHat with a set of gene model annotations and/or known transcripts, as a GTF 2.2 or GFF3 formatted file.
* -o: Output directory

Check the output of TopHat2:

cd tophat2
ls

Get the mapping rate:

cat align_summary.txt

Get the number of reads mapped.
Run the alignment of filtered data using TopHat2

cd ..
tophat2 -o tophat2_with_filtered_data --no-mixed \
--rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina \
-G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 \
/home/training/Data/03_Quality_control_and_data_preprocessing/Read_1_q30l50.fastq \
/home/training/Data/03_Quality_control_and_data_preprocessing/Read_2_q30l50.fastq

Check the output of TopHat2:

cd tophat2_with_filtered_data
ls

Get the mapping rate:

cat align_summary.txt

Get the number of reads mapped.
* What difference does using the filtered data make?

@@ Line 11: / Line 11: @@
 You will use the following software:
 * TopHat v2.1.1: http://ccb.jhu.edu/software/tophat/index.shtml. (This is the latest version)
 The data set you'll be using is downloaded from ENA (http://www.ebi.ac.uk/ena/data/view/SRP019027). The reads belong to sample SRR769316. The data set is tailored with respect to the time allocated for the exercise.
 = Indexing =

@@ Line 85: / Line 85: @@
   watch ls -altR
-You'll need <code>ctrl+c</code> to get out of the watch command.
+You'll need <code>ctrl+c</code> to get out of the <code>watch</code> command.
 If it has finished, we can go on to check the output of TopHat:

@@ Line 82: / Line 82: @@
 * Note the final three arguments, they do not have a switch: the index name, and our two reads.
-This will now take a while to run, some 10 minutes.
+This will now take a while to run, some 10 minutes. If you feel curious you can switch to another screen session, and see if it's actually writing any files. The <code>watch</code> command can be used like so:
 If it has finished, we can go on to check the output of TopHat:
@@ Line 100: / Line 103: @@
   ../01_Quality_Control_and_Preprocessing
-You may have used some different names for our filtered files, build up your symbolic link command line as follows
+You may have used some different names for your filtered files, build up your symbolic link command line as follows
   ln -s ../01<TAB>
 You should see suggestions. Keep hitting the TAB key until you get the file you want, then add a space and <code>.</code>. We have done this before, but please ask if you are still unclear on this command.

@@ Line 20: / Line 20: @@
 Go to the appropriate folder/directory:
-  cd $HOME/i2rda_data/Mapping_to_Reference
+  cd $HOME/i2rda_data/02_Mapping_to_Reference
-We keep the reference in a separte subfolder, we index the reference genome using one of the Bowtie2 tools: <code>bowtie2-build</code>
+We keep the reference in a separate subfolder, we index the reference genome using one of the Bowtie2 tools: <code>bowtie2-build</code>
-  cd Reference
+  cd Reference_files
   bowtie2-build mm10_chr19-1-20000000.fasta mm10_chr19-1-20000000
@@ Line 31: / Line 31: @@
 = Running the splice-aware aligner =
-Now the index is build, we can go on and verify we have tophat
+Now the index is built, we can go on and verify we have tophat
   module list
@@ Line 54: / Line 54: @@
 If in <code>screen</code>'s copy mode, type escape to get out. There is a reference to the website (http://ccb.jhu.edu/software/tophat/manual.shtml), where there is a certain comment:
-	Please note that it is highly recommended that a FASTA file with the sequence(s) the genome being indexed be present in the same directory with the Bowtie index files and having the name <genome_index_base>.fa. If not present, TopHat will automatically rebuild this FASTA file from the Bowtie index files.
+	"Please note that it is highly recommended that a reference FASTA file with the sequence(s) the genome being indexed be present in the same directory with the Bowtie index files and having the name <genome_index_base>.fa. If not present, TopHat will automatically rebuild this FASTA file from the Bowtie index files."
-Rather idiosyncratically, the <code>fa</code> extension is required for our reference genome name. But our's ends in fasta. What we can do is go back in and create a symlink in the following way:
+Rather idiosyncratically, the <code>fa</code> extension is required for our reference genome name. But our own ends in <code>fasta</code>. What we can do is go back in and create a symboic link. First make sure we are in the <code>Reference_files</code> subfolder:
   pwd
-  <user_home>/i2rda_data/Mapping_to_Reference/Reference
+  ~/i2rda_data/02_Mapping_to_Reference/Reference_files
 then
@@ Line 70: / Line 68: @@
   cd ..
-Run the alignment using TopHat2:
+Perhaps when you looked in this folder, you will have noticed that the read files themselves are symbolic links to the folder we used in the last session. SYmbolic links are very useful, they avoid copying large files, and help (somewhat) to administer files. Now let's starts running the alignment using TopHat:
-  tophat -o tophat2 --no-mixed --rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina -G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 $HOME/i2rda_data/Mapping_to_Reference/Read_1.fastq.gz $HOME/i2rda_data/Mapping_to_Reference/Read_2.fastq.gz
+  tophat -o tophat2 --no-mixed --rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina -G Reference_files/mm10_chr19-1-20000000_Ensembl.gtf Reference_files/mm10_chr19-1-20000000 Read_1.fastq.gz Read_2.fastq.gz
 <ins>where</ins>:
@@ Line 99: / Line 97: @@
 = Same procedure with filtered data =
-  cd ..
+You will have to symbolically link your filtered files from:
-  tophat -o tophat2_filtered --no-mixed --rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina -G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 $HOME/i2rda_data/Mapping_to_Reference/Read_1_q30l50.fastq.gz $HOME/i2rda_data/Mapping_to_Reference/Read_2_q30l50.fastq.gz
+  ../01_Quality_Control_and_Preprocessing
-After it's finished, check the output of TopHat2:
+After it's finished, check the output of TopHat:
   cd tophat2_filtered
@@ Line 113: / Line 126: @@
 <code>Question</code>:
 * What difference does using the filtered data make?

@@ Line 45: / Line 45: @@
   tophat --help
-That's rather long. Yo can use <code>ctrl+l,Escape</code> to enter screen scrolling mode, also called copy-mode, and then use <code>PgUp</code> but the template for launching it is here:
+That's rather long. Yo can use <code>ctrl+l,Escape</code> to enter <code>screen</code>'s scrolling mode, also called copy-mode, and then use <code>PgUp</code> but the template for launching it is here:
   tophat [options] <bowtie_index> <reads1[,reads2,...]> [reads1[,reads2,...]] [quals1,[quals2,...]] [quals1[,quals2,...]]
@@ Line 52: / Line 52: @@
 So, typically we will have some options and then <bowtie_index> and <reads1... These are not surrounded by square brackets, which means they are obligatory. Also, because they do not have a switch, such as -x or something like that, they need to be in order. OK, now we can look at the full command:
-If in screen's copy mode, type escape to get out. There is a reference to the website (http://ccb.jhu.edu/software/tophat/manual.shtml), where there is a certain comment:
+If in <code>screen</code>'s copy mode, type escape to get out. There is a reference to the website (http://ccb.jhu.edu/software/tophat/manual.shtml), where there is a certain comment:
 	Please note that it is highly recommended that a FASTA file with the sequence(s) the genome being indexed be present in the same directory with the Bowtie index files and having the name <genome_index_base>.fa. If not present, TopHat will automatically rebuild this FASTA file from the Bowtie index files.
@@ Line 58: / Line 58: @@
 Rather idiosyncratically, the <code>fa</code> extension is required for our reference genome name. But our's ends in fasta. What we can do is go back in and create a symlink in the following way:
-First make sure we are in the reference subfolder
+First make sure we are in the <code>Reference</code> subfolder:
   pwd
@@ Line 86: / Line 86: @@
 This will now take a while to run, some 10 minutes.
-If it has finished, we can go on to check the output of TopHat2:
+If it has finished, we can go on to check the output of TopHat:
   cd tophat2

@@ Line 10: / Line 10: @@
 You will use the following software:
-* TopHat2 v2.0.11: http://ccb.jhu.edu/software/tophat/index.shtml
+* TopHat v2.1.1: http://ccb.jhu.edu/software/tophat/index.shtml. (This is the latest version)
-* Bowtie2 v2.2.0: http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
+:- loaded via <code>module load tophat</code>
 The data set you'll be using is downloaded from ENA (http://www.ebi.ac.uk/ena/data/view/SRP019027). The reads belong to sample SRR769316. The data set is tailored with respect to the time allocated for the exercise.
@@ Line 17: / Line 18: @@
 = Indexing =
-Go to the right folder/directory:
+Go to the appropriate folder/directory:
   cd $HOME/i2rda_data/Mapping_to_Reference
-Index the reference genome using one of the Bowtie2:
+We keep the reference in a separte subfolder, we index the reference genome using one of the Bowtie2 tools: <code>bowtie2-build</code>
   cd Reference
-  bowtie2-build mm10_chr19-1-20000000.fa mm10_chr19-1-20000000
+  bowtie2-build mm10_chr19-1-20000000.fasta mm10_chr19-1-20000000
-Note the "fa" extension for the reference this is due to a preference of tophat which we'll be using below.
+Note the final argument is the index name we give it
 Run the alignment using TopHat2:
   tophat -o tophat2 --no-mixed --rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina -G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 $HOME/i2rda_data/Mapping_to_Reference/Read_1.fastq.gz $HOME/i2rda_data/Mapping_to_Reference/Read_2.fastq.gz
-where:
+<ins>where</ins>:
-* --no-mixed: For paired reads, only report read alignments if both reads in a pair can be mapped
+* --no-mixed: Important: for paired reads, only report read alignments if both reads in a pair can be mapped
-* --rg-id: Read group ID
+* --rg-id: Read group ID, for now we use an invented one.
-* --rg-sample: Sample ID
+* --rg-sample: Sample ID, for now we use an invented one.
-* --rg-center: Sequencing Centre name
+* --rg-center: Sequencing Centre name, for now we use an invented one.
 * --rg-platform: Sequencing platform descriptor
 * -G: Supply TopHat with a set of gene model annotations and/or known transcripts, as a GTF 2.2 or GFF3 formatted file.
 * -o: Output directory
-Check the output of TopHat2:
+If it has finished, we can go on to check the output of TopHat2:
   cd tophat2
   ls
-Get the mapping rate:
+Make a note of the mapping rate:
   cat align_summary.txt
 Get the number of reads mapped.
-Run the alignment of filtered data using TopHat2
   cd ..
-  tophat2 -o tophat2_with_filtered_data --no-mixed \
+  tophat -o tophat2_filtered --no-mixed --rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina -G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 $HOME/i2rda_data/Mapping_to_Reference/Read_1_q30l50.fastq.gz $HOME/i2rda_data/Mapping_to_Reference/Read_2_q30l50.fastq.gz
-Check the output of TopHat2:
+After it's finished, check the output of TopHat2:
-  cd tophat2_with_filtered_data
+  cd tophat2_filtered
   ls
-Get the mapping rate:
+Again have a look at the mapping rate & the number of reads mapped
   cat align_summary.txt
-Get the number of reads mapped.
+<code>Question</code>:
 * What difference does using the filtered data make?

@@ Line 32: / Line 32: @@
   --rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina \
   -G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 \
-  $HOME/i2rda_data/Mapping_to_Reference/Read_1.fastq \
+  $HOME/i2rda_data/Mapping_to_Reference/Read_1.fastq.gz \
-  $HOME/i2rda_data/Mapping_to_Reference/Read_2.fastq
+  $HOME/i2rda_data/Mapping_to_Reference/Read_2.fastq.gz
 where:

@@ Line 1: / Line 1: @@
 = Aims =
 In this part you will learn to:
@@ Line 13: / Line 15: @@
 The data set you'll be using is downloaded from ENA (http://www.ebi.ac.uk/ena/data/view/SRP019027). The reads belong to sample SRR769316. The data set is tailored with respect to the time allocated for the exercise.
-Change directory:
+= Indexing =
-  cd /home/training/Data/04_Mapping_to_a_reference_genome
+  cd $HOME/i2rda_data/Mapping_to_Reference
-Index the reference genome using Bowtie2:
+Index the reference genome using one of the Bowtie2:
   cd Reference
@@ Line 28: / Line 32: @@
   --rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina \
   -G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 \
-  /home/training/Data/03_Quality_control_and_data_preprocessing/Read_1.fastq \
+  $HOME/i2rda_data/Mapping_to_Reference/Read_1.fastq \
-  /home/training/Data/03_Quality_control_and_data_preprocessing/Read_2.fastq
+  $HOME/i2rda_data/Mapping_to_Reference/Read_2.fastq
 where:
@@ Line 56: / Line 60: @@
   --rg-id Lane-1 --rg-sample sample1 --rg-center XYZ --rg-platform Illumina \
   -G Reference/mm10_chr19-1-20000000_Ensembl.gtf Reference/mm10_chr19-1-20000000 \
-  /home/training/Data/03_Quality_control_and_data_preprocessing/Read_1_q30l50.fastq \
+  $HOME/i2rda_data/Mapping_to_Reference/Read_1_q30l50.fastq \
-  /home/training/Data/03_Quality_control_and_data_preprocessing/Read_2_q30l50.fastq
+  $HOME/i2rda_data/Mapping_to_Reference/Read_2_q30l50.fastq
 Check the output of TopHat2: