Quality of Mapping Exercise - Revision history

Rf at 14:09, 14 May 2017

2017-05-14T14:09:51Z

Rf at 13:49, 14 May 2017

2017-05-14T13:49:53Z

Rf at 13:46, 14 May 2017

2017-05-14T13:46:17Z

Rf at 20:16, 11 May 2017

2017-05-11T20:16:37Z

Rf at 18:18, 11 May 2017

2017-05-11T18:18:02Z

Rf at 18:24, 8 May 2017

2017-05-08T18:24:27Z

Rf at 18:17, 8 May 2017

2017-05-08T18:17:52Z

Rf at 00:16, 8 May 2017

2017-05-08T00:16:40Z

Rf at 18:26, 7 May 2017

2017-05-07T18:26:27Z

Rf at 18:23, 7 May 2017

2017-05-07T18:23:26Z

@@ Line 37: / Line 37: @@
   cd ~/i2rda_data/03_Quality_of_Mapping
-We should see the file SRR769316.bam in there, a binary alignment file. We want to have a quick look at it using <code>less</code>, but <code>bam</code> is a binary file. It will throw gibberish on the screen, and may even reset your character encoding, which woul dforce you to quite the session. So, we also need <code>samtools view</code> to view it properly. Also, it has very long lines so we will use a special option of <code>less</code>, the <code>-S</code> option.
+We should see the file <code>SRR769316.bam</code> in there, a binary alignment file. We want to have a quick look at it using <code>less</code>, but <code>bam</code> is a binary file. It will throw gibberish on the screen, and may even reset your character encoding, which would force you to quite the session. So, we also need <code>samtools view</code> to view it properly. Also, it has very long lines so we will use a special option of <code>less</code>, the <code>-S</code> option.
   samtools view SRR769316.bam | less -S

@@ Line 16: / Line 16: @@
 We will use the following software tools:
 * samtools v0.1.19: http://samtools.sourceforge.net/
-* BamQC git commit version 480c091
+* BamQC git commit version 480c091: https://github.com/s-andrews/BamQC
 * Picard v2.7.1  http://picard.sourceforge.net/
 * RSeQC v2.6.4: http://rseqc.sourceforge.net/
@@ Line 51: / Line 51: @@
 One concern is to quickly see how many of the reads have aligned. We've talked about the SAM format's flag value, well 4 is the value we want to check.
-  samtools view -F 4 SRR769316.bam
+  samtools view -F 4 -c SRR769316.bam
 Unmapped reads are given by
-  samtools view -f 4 SRR769316.bam
+  samtools view -f 4 -c SRR769316.bam
 Simon Andrews' lab gave us FastQC which we saw in the beginning, and BamQC is also from them and has a very similar way of working:
@@ Line 116: / Line 116: @@
   samtools index SRR769316.bam
-This will automatically create a <code>SRR769316.bai</code> file
+This will automatically create a <code>SRR769316.bam.bai</code> file
 Now we ca go ahead an launch the script.

@@ Line 88: / Line 88: @@
 RSeQC's <code>inner_distance.py</code> script calculates the inner distance (or insert size) between two paired RNA reads:
-  inner_distance.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
+  inner_distance.py -i SRR769316.bam -r Reference_files/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
 <ins>where</ins>:
@@ Line 120: / Line 120: @@
 Now we ca go ahead an launch the script.
-  geneBody_coverage.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
+  geneBody_coverage.py -i SRR769316.bam -r Reference_files/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
 When it's finished, open the plot:
@@ Line 133: / Line 133: @@
 RSeQC's <code>junction_annotation.py</code> script compares detected splice junctions to reference gene models:
-  junction_annotation.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
+  junction_annotation.py -i SRR769316.bam -r Reference_files/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
 Splicing annotation is performed at two levels: splice event level and splice junction level. All detected junctions can be grouped into 3 exclusive categories:
@@ Line 150: / Line 150: @@
 RSeQC's <code>junction_saturation.py</code> script checks for saturation by resampling 5%, 10%, 15%, ..., 95% of total alignments, and then detects splice junctions from each subset and compares them to reference gene model:
-  junction_saturation.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316 -v 5
+  junction_saturation.py -i SRR769316.bam -r Reference_files/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316 -v 5
 <ins>where</ins>:

@@ Line 16: / Line 16: @@
 We will use the following software tools:
 * samtools v0.1.19: http://samtools.sourceforge.net/
 * Picard v2.7.1  http://picard.sourceforge.net/
 * RSeQC v2.6.4: http://rseqc.sourceforge.net/
@@ Line 21: / Line 22: @@
 = Loading the software =
 * RSeQC is a python program and python is installed by default, so it needs no loading.
-* The other two are loaded as follows:
+* The other three are loaded as follows:
-  module load samtools picard-tools
+  module load samtools BamQC picard-tools
 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. Reads were aligned to the first 20 Mb of chromosome 19 of the mouse reference genome (GRCm38/mm10) using TopHat.
@@ Line 46: / Line 47: @@
 You can navigate through a line by using the <code>[Right arrow]</code> and <code>[Left arrow]</code> keys. You can navigate through the file line-by-line using the
 <code>Up-arrow</code> and <code>Down-arrow</code> keys, or advance page-by-page using the <code>Enter</code> key, <code>[Space</code>-bar and <code>B</code>-key. You can exit the file using the <code>Q</code>-key.
 = Marking duplicates =

@@ Line 9: / Line 9: @@
 In this part you will learn to:
 * mark duplicated read pairs
 * check the fragment size distribution
 * check the read coverage over the gene body
 * check if we have enough data
-You will use the following software tools:
+We will use the following software tools:
 * samtools v0.1.19: http://samtools.sourceforge.net/
 * Picard v2.7.1  http://picard.sourceforge.net/
@@ Line 26: / Line 27: @@
 You will use the following files:
-m10_chr19-1-20000000_Ensembl.bed: Ensembl mouse gene models
+* <code>m10_chr19-1-20000000_Ensembl.bed</code>: Ensembl mouse gene models
 = Data =
-The data is available in the directory <code>)03_Quality_of_Mapping</code>:
+The data is available in the directory <code>03_Quality_of_Mapping</code>:
   cd ~/i2rda_data/03_Quality_of_Mapping
@@ Line 103: / Line 105: @@
 When it's finished, open the plot:
-  xpdf SRR769316.geneBodyCoverage.pdf &
+  xpdf SRR769316.geneBodyCoverage.curves.pdf &
 <ins>Questions</ins>:
@@ Line 114: / Line 116: @@
   junction_annotation.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
-Splicing annotation is performed at two levels: splice event level and splice junction level.
+Splicing annotation is performed at two levels: splice event level and splice junction level. All detected junctions can be grouped into 3 exclusive categories:
 All detected junctions can be grouped into 3 exclusive categories:
 # annotated: Known junction. Both 5’ and 3' splice site are present in the reference gene model.
 # complete_novel: Complete new junction. Neither of the two splice sites is present in the reference gene model.

@@ Line 26: / Line 26: @@
 You will use the following files:
-* SRR769316.bam: aligned reads
+m10_chr19-1-20000000_Ensembl.bed: Ensembl mouse gene models
 = Data =
@@ Line 132: / Line 131: @@
   junction_saturation.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316 -v 5
-where:
+<ins>where</ins>:
 * <code>-v</code>: minimum number of supporting reads to call a junction (default=1)
 This shows if more sequencing will enable the discovery of new events.
@@ Line 140: / Line 140: @@
   xpdf SRR769316.junctionSaturation_plot.pdf &
-Do you think more sequencing is needed?
+<ins>Question</ins>:
-Some issues are only detectable in the context of the
+Some issues are only detectable in the context of the genome:
 * Duplicate reads

@@ Line 31: / Line 31: @@
 = Data =
-The data is available in the directory <code>Quality_of_Mapping</code>:
+The data is available in the directory <code>)03_Quality_of_Mapping</code>:
-  cd ~/i2rda_data/Quality_of_Mapping
+  cd ~/i2rda_data/03_Quality_of_Mapping
-We want to have a quick look at the file SRR769316.bam using less, but bam is a binary file. It will throw gibberish on the screen. So, we need <code>samtools view</code>. Also, it has very long lines so we will use a special option of <code>less</code>, the <code>-S</code> option.
+We should see the file SRR769316.bam in there, a binary alignment file. We want to have a quick look at it using <code>less</code>, but <code>bam</code> is a binary file. It will throw gibberish on the screen, and may even reset your character encoding, which woul dforce you to quite the session. So, we also need <code>samtools view</code> to view it properly. Also, it has very long lines so we will use a special option of <code>less</code>, the <code>-S</code> option.
   samtools view SRR769316.bam | less -S
@@ Line 44: / Line 44: @@
 You can navigate through a line by using the <code>[Right arrow]</code> and <code>[Left arrow]</code> keys. You can navigate through the file line-by-line using the
-<code>Up-arrow</code> and <code>Down-arrow</code> keys, or advance page-by-page using the <code>Enter</code> key, <code>[Space</code>-bar and <code>B</code>-key. You can exit the file using the <code>Q-</code>key.
+<code>Up-arrow</code> and <code>Down-arrow</code> keys, or advance page-by-page using the <code>Enter</code> key, <code>[Space</code>-bar and <code>B</code>-key. You can exit the file using the <code>Q</code>-key.
 = Marking duplicates =
@@ Line 57: / Line 57: @@
 * <code>M=</code>: file to write duplication metrics to
-MarkDuplicates creates two files, i.e. a BAM file with all duplicate records flagged, and a file containing the duplication metrics.
+This manner of presenting inputs and outputs to a program is quite particular to <code>picard-tools</code>, and we also see in its widely used sister program <code>GATK</code>. <code>MarkDuplicates</code> creates two files, i.e. a <code>bam</code> file with all duplicate records flagged, and a file containing the duplication metrics.
 * What is the duplicate rate of sample SRR769316?
-Try:
+To get the result quickly, try:
   awk 'BEGIN{FS="\t"}{if(NF>9) printf "%20s\t%20s\t%20s\n",$3,$7,$9}' SRR769316_duplicates.metrics.csv
@@ Line 87: / Line 87: @@
   xpdf SRR769316.inner_distance_plot.pdf &
 * Knowing that the read length is 99 bases what is the mean length of the fragments?
 = Checking the read coverage over the gene body =
-RSeQC's geneBody_coverage.py script reports the read coverage of all genes in the percentile of their length:
+RSeQC's geneBody_coverage.py script reports the read coverage of all genes in the percentile of their length. For its intensive search however, it would like us to generate an index first:
   geneBody_coverage.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
-Open the plot:
+When it's finished, open the plot:
   xpdf SRR769316.geneBodyCoverage.pdf &
 * Is there a 3' or 5' bias?
 * How would you explain the bell shaped curve?
-Checking if we have enough data
+= Checking if we have enough data =
-RSeQC's junction_annotation.py script compares detected splice junctions to reference gene models:
+RSeQC's <code>junction_annotation.py</code> script compares detected splice junctions to reference gene models:
-  junction_annotation.py -i SRR769316.bam \
+  junction_annotation.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
 Splicing annotation is performed at two levels: splice event level and splice junction level.
@@ Line 114: / Line 121: @@
 # partial_novel: Partially new junction. Only one of the splice sites is present in the reference gene model.
-Open the plots:
+Now open the plots, we use a heavier pdf viewer now to open two plots at the same time.
-  xpdf SRR769316.splice_*.pdf &
+  evince SRR769316.splice_*.pdf &
-Is the data mostly supporting known or novel junctions?
+<ins>Question</ins>:
-RSeQC's junction_saturation.py script checks for saturation by resampling 5%, 10%, 15%, ..., 95% of total alignments, and then detects splice junctions from each subset and compares them to reference gene model:
+RSeQC's <code>junction_saturation.py</code> script checks for saturation by resampling 5%, 10%, 15%, ..., 95% of total alignments, and then detects splice junctions from each subset and compares them to reference gene model:
-  junction_saturation.py -i SRR769316.bam \
+  junction_saturation.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316 -v 5
 where:

@@ Line 48: / Line 48: @@
 = Marking duplicates =
-Picard is a software suite programmed in Java. It started out a a "samtools for Java", so it is very similar, but it has evolved quite alot. Its <code>MarkDuplicates</code> subcommand is very widely used. It marks all duplicate read pairs in our dataset. We can find out how it may be used with
+Picard is a software suite programmed in Java. It started out as a "samtools for Java", so it is very similar, but it has evolved quite alot. Its <code>MarkDuplicates</code> subcommand is very widely used. It marks all duplicate read pairs in our dataset. We can find out how it may be used with
   java -jar $PICARDJARPATH/picard.jar MarkDuplicates I=SRR769316.bam O=SRR769316_duplicates_marked.bam M=SRR769316_duplicates.metrics.csv
@@ Line 58: / Line 58: @@
 MarkDuplicates creates two files, i.e. a BAM file with all duplicate records flagged, and a file containing the duplication metrics.
- gnumeric SRR769316_duplicates.metrics.csv &
+* What is the duplicate rate of sample SRR769316?
-* What is the duplicate rate of sample SRR769316?
+Try:
 = Checking the fragment size distribution =
-inner_distance.py -i SRR769316.bam \
+RSeQC's <code>inner_distance.py</code> script calculates the inner distance (or insert size) between two paired RNA reads:
 -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
 <ins>where</ins>:
@@ Line 75: / Line 75: @@
 * -o: prefix of output file(s)
-The file with reference gene models is used to correct for any introns that are present between the reads.
+The <code>bed</code>-file with reference gene models is used to correct for any introns that are present between the reads.
-The script first determines the genomic (DNA) size between two paired reads: D_size = read2_start - read1_end, then
+The script first determines the genomic DNA size (D_size) between two paired reads: D_size = read2_start - read1_end, then
 * if two paired reads map to the same exon: inner distance = D_size
 * if two paired reads map to different exons: inner distance = D_size - intron_size
@@ Line 89: / Line 89: @@
 * Knowing that the read length is 99 bases what is the mean length of the fragments?
-Checking the read coverage over the gene body
+= Checking the read coverage over the gene body =
 RSeQC's geneBody_coverage.py script reports the read coverage of all genes in the percentile of their length:
-geneBody_coverage.py -i SRR769316.bam \
+ geneBody_coverage.py -i SRR769316.bam -r Reference/mm10_chr19-1-20000000_Ensembl.bed -o SRR769316
 Open the plot: