-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy path99-RNAseq_DE_analysis_with_R.html
More file actions
748 lines (734 loc) · 59.8 KB
/
99-RNAseq_DE_analysis_with_R.html
File metadata and controls
748 lines (734 loc) · 59.8 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<meta name="generator" content="pandoc">
<title>RNAseq Using R </title>
<link rel="shortcut icon" type="image/x-icon" href="/favicon.ico" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
<link rel="stylesheet" type="text/css" href="css/bootstrap/bootstrap.css" />
<link rel="stylesheet" type="text/css" href="css/bootstrap/bootstrap-theme.css" />
<link rel="stylesheet" type="text/css" href="css/swc.css" />
<meta charset="UTF-8" />
<!-- HTML5 shim, for IE6-8 support of HTML5 elements -->
<!--[if lt IE 9]>
<script src="http://html5shim.googlecode.com/svn/trunk/html5.js"></script>
<![endif]-->
</head>
<body class="lesson">
<div class="container card">
<div class="banner" align="right" style="margin-top: 5px">
<a href="http://monash.edu/bioinformatics" title="Monash Bioinformatics Platform">
<div style="margin-right: 5px;" class="label swc-blue-bg">Monash Bioinformatics Platform</div>
</a>
</div>
<article>
<div class="row">
<div class="col-md-10 col-md-offset-1">
<h1 class="title"></h1>
<h1 id="rna-seq-differential-gene-expression-analysis">RNA-seq: differential gene expression analysis</h1>
<section class="objectives panel panel-warning">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-certificate"></span>Learning Objectives</h2>
</div>
<div class="panel-body">
<p>This course is an introduction to differential expression analysis from RNAseq data. It will take you from the raw fastq files all the way to the list of differentially expressed genes, via the mapping of the reads to a reference genome and statistical analysis using the limma package.</p>
<p>The Monash Bioinformatics Platform thanks and acknowledges QFAB (<a href="http://www.qfab.org">http://www.qfab.org</a>) for original materials in this section.</p>
</div>
</section>
<h3 id="install-and-load-packages">Install and load packages</h3>
<p>Most generic R packages are hosted on the Comprehensive R Archive Network (CRAN, <a href="http://cran.us.r-project.org/">http://cran.us.r-project.org/</a>). To install one of these packages, you would use <code>install.packages("packagename")</code>. You only need to install a package once, then load it each time using <code>library(packagename)</code>.</p>
<p>Bioconductor packages work a bit differently, and are not hosted on CRAN. Go to <a href="http://bioconductor.org/">http://bioconductor.org/</a> to learn more about the Bioconductor project. To use any Bioconductor package, you’ll need a few “core” Bioconductor packages. Run the following commands to (1) download the installer script, and (2) install some core Bioconductor packages. You’ll need internet connectivity to do this, and it’ll take a few minutes, but it only needs to be done once.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Download the installer script</span>
<span class="kw">source</span>(<span class="st">"http://bioconductor.org/biocLite.R"</span>)
<span class="co"># biocLite() is the bioconductor installer function. Run it without any</span>
<span class="co"># arguments to install the core packages or update any installed packages. This</span>
<span class="co"># requires internet connectivity and will take some time!</span>
<span class="kw">biocLite</span>()</code></pre>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pushpin"></span>NOTE</h2>
</div>
<div class="panel-body">
<p>All packages required for this course have already been installed on our training server so you won’t need to run this part today.</p>
</div>
</aside>
<h3 id="data-pre-processing">Data pre-processing</h3>
<p>To start with, delete all previously saved R objects and define your working directory for the RNAseq data analysis.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Delete all previously saved R objects</span>
<span class="kw">rm</span>(<span class="dt">list=</span><span class="kw">ls</span>())</code></pre>
<p>The data considered for the RNAseq part of the course have been downloaded from ArrayExpress (<a href="http://www.ebi.ac.uk/arrayexpress">http://www.ebi.ac.uk/arrayexpress</a>) and correspond to 8 RNA sequencing libraries from Human brain and liver.</p>
<p>Raw sequencing data are usually available in FASTQ format which is a well defined text-based format for storing both biological sequences (usually nucleotide sequences) and their corresponding quality scores. The raw data from this study have been downloaded (8Gb / fastq file) into the shared directory “~/data/RNAseq/raw_data”.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># define shared directory for RNAseq data</span>
RNAseqDATADIR <-<span class="st"> "data/raw_data"</span>
<span class="co">#list the fastq files in the raw data directory</span>
<span class="kw">dir</span>(RNAseqDATADIR)</code></pre>
<pre><code>## [1] "ERR420388_mini_1.fastq.gz" "ERR420388_mini_2.fastq.gz"
## [3] "ERR420388_subsamp_1.fastq.gz" "ERR420388_subsamp_2.fastq.gz"
## [5] "experiment_design.txt"</code></pre>
<p>The first step in a RNAseq analysis is to run a quick quality check on your data, this will give you an idea of the quality of your raw data in terms of number of reads per library, read length, average quality score along the reads, GC content, sequence duplication level, adaptors that might have not been removed correctly from the data etc.</p>
<p>The fastQC tool is quick and easy to run and can be downloaded from here: <a href="http://www.bioinformatics.babraham.ac.uk/projects/fastqc/">http://www.bioinformatics.babraham.ac.uk/projects/fastqc/</a>.</p>
<p>To ensure highest quality of the sequences for subsequent mapping and differential expression analysis steps, the reads can also be trimmed using the Trimmomatic tool (Lohse et al. 2012, <a href="http://www.usadellab.org/cms/?page=trimmomatic">http://www.usadellab.org/cms/?page=trimmomatic</a>).</p>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pushpin"></span>NOTE</h2>
</div>
<div class="panel-body">
<p>For the scope of this course we will focus on the R-based steps and will assume that the data are fit for purpose.</p>
<p>Also the next two sections (Mapping and Read counts) will be performed on small subsets of those files due to time and I/O limitations.</p>
</div>
</aside>
<h3 id="mapping-reads-to-a-reference-genome">Mapping reads to a reference genome</h3>
<p>Once the reads have been quality checked and trimmed, the next step is to map the reads to the reference genome (in our case the human genome “hg19”). This can be done with the Bioconductor package <em>Rsubread</em> (Y et al. 2013).</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(Rsubread)</code></pre>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pushpin"></span>NOTE</h2>
</div>
<div class="panel-body">
<p>Please remember that you can also perform those steps using command line tools (such as BWA and featureCounts) as described in the Unix shell course.</p>
</div>
</aside>
<p>If the package is not already installed on your system, you can install it by typing:</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">source</span>(<span class="st">"http://bioconductor.org/biocLite.R"</span>)
<span class="kw">biocLite</span>(<span class="st">"Rsubread"</span>)</code></pre>
<p>Rsubread provides reference genome indices for the most common organisms: human and mouse. If you are working with a different organism you can build your own index using the <em>buildindex</em> command.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># define the reference genome fasta file</span>
REF_GENOME <-<span class="st"> "hg19.fa"</span>
<span class="co"># define the output directory for the Rsubread index</span>
<span class="co"># (admin note: requires data/ref_data/download_hg19.sh to be run first)</span>
RSUBREAD_INDEX_PATH <-<span class="st"> "data/ref_data"</span>
<span class="co"># define the basename for the index</span>
RSUBREAD_INDEX_BASE <-<span class="st"> "hg19"</span>
<span class="co"># check what is in the reference directory</span>
<span class="kw">dir</span>(RSUBREAD_INDEX_PATH)</code></pre>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pushpin"></span>NOTE</h2>
</div>
<div class="panel-body">
<p>For the purpose of this course the index has already been built and the mapping will be done on small subset of reads since this step can take up to a couple of hours per library. <strong>Please do not run the command below</strong>.</p>
</div>
</aside>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># build the index</span>
<span class="kw">buildindex</span>(<span class="dt">basename=</span><span class="kw">file.path</span>(RSUBREAD_INDEX_PATH,RSUBREAD_INDEX_BASE), <span class="dt">reference=</span>REF_GENOME)</code></pre>
<p>Once the Rsubread index has been created you can map your reads to the genome by running the <em>align</em> command.</p>
<p>The code below will be used to map the reads for a specific library against the genome for which the index has been built.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># list files in the raw data directory</span>
<span class="kw">dir</span>(RNAseqDATADIR)
<span class="co"># define the fastq file with forward reads</span>
inputfilefwd <-<span class="st"> </span><span class="kw">file.path</span>(RNAseqDATADIR,<span class="st">"ERR420388_subsamp_1.fastq.gz"</span>)
<span class="co"># define the fastq file with reverse reads</span>
inputfilervs <-<span class="st"> </span><span class="kw">file.path</span>(RNAseqDATADIR,<span class="st">"ERR420388_subsamp_2.fastq.gz"</span>)
<span class="co"># run the align command to map the reads</span>
<span class="kw">align</span>(<span class="dt">index=</span><span class="kw">file.path</span>(RSUBREAD_INDEX_PATH,RSUBREAD_INDEX_BASE), <span class="dt">readfile1=</span>inputfilefwd, <span class="dt">readfile2=</span>inputfilervs, <span class="dt">output_file=</span><span class="st">"ERR420388.sam"</span>, <span class="dt">output_format=</span><span class="st">"SAM"</span>)</code></pre>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pushpin"></span>NOTE:</h2>
</div>
<div class="panel-body">
<p>The nthreads parameter can be used in the <em>align</em> command to speed up the process and run the alignment using several CPUs in parallel.</p>
</div>
</aside>
<p>The function <em>propmapped</em> returns the proportion of mapped reads in the output SAM file: total number of input reads, number of mapped reads and proportion of mapped reads. Let’s have a look at a small SAM file as an example:</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># define the path to SAM file </span>
outputsamfile <-<span class="st"> "data/mapping/ERR420388.sam"</span>
<span class="kw">propmapped</span>(outputsamfile)</code></pre>
<h3 id="count-reads-for-each-feature">Count reads for each feature</h3>
<p>Rsubread provides a read summarization function <em>featureCounts</em>, which takes as input the SAM or BAM files and assigns them to genomic features. This gives the number of reads mapped per gene, which can then be transformed into RPKM values (Read Per Killobase per Million), normalised and tested for differential expression.</p>
<p>For the purpose of this course we will use the index previously built and its corresponding GTF file to identify and count the reads mapped to each feature (gene).</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Getting read counts using the index previously built</span>
mycounts<-<span class="kw">featureCounts</span>(outputsamfile, <span class="dt">annot.ext=</span><span class="kw">file.path</span>(RSUBREAD_INDEX_PATH,<span class="st">"hg19.genes.gtf"</span>), <span class="dt">isGTFAnnotationFile=</span><span class="ot">TRUE</span>, <span class="dt">isPairedEnd=</span><span class="ot">TRUE</span>)
<span class="co"># Checking your output object</span>
<span class="kw">summary</span>(mycounts)
<span class="kw">dim</span>(mycounts$counts)
<span class="kw">head</span>(mycounts$annotation)
mycounts$targets
mycounts$stat</code></pre>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pushpin"></span>NOTE:</h2>
</div>
<div class="panel-body">
<p>In-built annotations for mm9, mm10 and hg19 are also provided for users’ convenience, but we won’t be using it for this course. <strong>Please do not run the command below</strong></p>
<p>Getting read counts using the inbuilt hg19 annotation:</p>
<p>mycounts.inbuilt <- featureCounts(myoutputfile, annot.inbuilt=“hg19”, isGTFAnnotationFile=FALSE, isPairedEnd=TRUE)</p>
</div>
</aside>
<p>For the purpose of this course the read summarisation step has already been performed for all libraries. You will need to load the corresponding RData file to get these read counts.</p>
<pre class="sourceCode r"><code class="sourceCode r">MAPPINGDIR <-<span class="st"> "data/mapping"</span>
<span class="co"># load the counts previously calculated</span>
<span class="kw">load</span>(<span class="kw">file.path</span>(MAPPINGDIR,<span class="st">"RawCounts.RData"</span>))
<span class="co"># check the presence of read counts for the 8 libraries</span>
<span class="kw">summary</span>(counts)</code></pre>
<pre><code>## Length Class Mode
## counts 205616 -none- numeric
## annotation 2 data.frame list
## targets 8 -none- character</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">counts$targets</code></pre>
<pre><code>## [1] "/data/Intro_to_R/RNAseq/mapping/ERR420386.sam"
## [2] "/data/Intro_to_R/RNAseq/mapping/ERR420387.sam"
## [3] "/data/Intro_to_R/RNAseq/mapping/ERR420388.sam"
## [4] "/data/Intro_to_R/RNAseq/mapping/ERR420389.sam"
## [5] "/data/Intro_to_R/RNAseq/mapping/ERR420390.sam"
## [6] "/data/Intro_to_R/RNAseq/mapping/ERR420391.sam"
## [7] "/data/Intro_to_R/RNAseq/mapping/ERR420392.sam"
## [8] "/data/Intro_to_R/RNAseq/mapping/ERR420393.sam"</code></pre>
<p>You can then print out these counts in a text file for future reference.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># print out counts table for every sample</span>
<span class="kw">write.table</span>(counts$counts,<span class="dt">file=</span><span class="st">"~/raw_read_counts.txt"</span>,<span class="dt">sep=</span><span class="st">"</span><span class="ch">\t</span><span class="st">"</span>, <span class="dt">quote=</span>F,<span class="dt">append=</span>F)</code></pre>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pushpin"></span>NOTE:</h2>
</div>
<div class="panel-body">
<p>To save disk space on your machine, remember to convert your SAM files to BAM format using samtools and to compress the fastq files. This can be done out of Rstudio, in a shell window:</p>
<p><code>samtools view -Shb myfile.sam -o myfile.bam</code></p>
<p><code>gzip myfile.fq</code></p>
<p>Visualisation of the BAM files you created can also be done via a genome browser such as IGV (<a href="http://www.broadinstitute.org/igv/">http://www.broadinstitute.org/igv/</a>) after sorting and indexing of those files.</p>
</div>
</aside>
<h3 id="qc-and-stats">QC and stats</h3>
<p>Rsubread provides the number of reads mapped to each gene which can then be used for ploting quality control figures and for differential expression analysis.</p>
<p>QC figures of the mapped read counts can be plotted and investigated for potential outlier libraries and to confirm grouping of samples.</p>
<p>Before plotting QC figures it is useful to get the experiment design. This will allow labeling the data with the sample groups they belong to, or any other parameter of interest.</p>
<p>The experiment design file corresponding to this study has been downloaded from the ArrayExpress webpage and formatted as a tab separated file for this analysis purposes. You can find it in the shared data directory.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># define the experiment design file (tab separated text file is best)</span>
EXPMT_DESIGN_FILE <-<span class="st"> </span><span class="kw">file.path</span>(RNAseqDATADIR,<span class="st">'experiment_design.txt'</span>)
<span class="co"># read the experiment design file and save it into memory</span>
experiment_design<-<span class="kw">read.table</span>(EXPMT_DESIGN_FILE,<span class="dt">header=</span>T,<span class="dt">sep=</span><span class="st">"</span><span class="ch">\t</span><span class="st">"</span>)
<span class="co">#</span>
<span class="co"># set the rownames to the sampleID to allow for ordering</span>
<span class="kw">rownames</span>(experiment_design) <-<span class="st"> </span>experiment_design$SampleID
<span class="co"># order the design following the counts sample order</span>
experiment_design.ord <-<span class="st"> </span>experiment_design[<span class="kw">colnames</span>(counts$counts),]
<span class="co"># look at the design</span>
experiment_design.ord</code></pre>
<pre><code>## SampleID Source_Name organism sex age tissue
## ERR420386 ERR420386 brain_sample_1 Homo_sapiens male 26 brain
## ERR420387 ERR420387 brain_sample_1 Homo_sapiens male 26 brain
## ERR420388 ERR420388 liver_sample_1 Homo_sapiens male 30 liver
## ERR420389 ERR420389 liver_sample_1 Homo_sapiens male 30 liver
## ERR420390 ERR420390 liver_sample_1 Homo_sapiens male 30 liver
## ERR420391 ERR420391 brain_sample_1 Homo_sapiens male 26 brain
## ERR420392 ERR420392 brain_sample_1 Homo_sapiens male 26 brain
## ERR420393 ERR420393 liver_sample_1 Homo_sapiens male 30 liver
## Extract_Name Material_Type Assay_Name technical_replicate_group
## ERR420386 GCCAAT RNA Assay4 group_2
## ERR420387 ACAGTG RNA Assay2 group_1
## ERR420388 GTGAAA RNA Assay7 group_4
## ERR420389 GTGAAA RNA Assay8 group_4
## ERR420390 CTTGTA RNA Assay6 group_3
## ERR420391 ACAGTG RNA Assay1 group_1
## ERR420392 GCCAAT RNA Assay3 group_2
## ERR420393 CTTGTA RNA Assay5 group_3</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># list the ordered samples for future use</span>
samples <-<span class="st"> </span><span class="kw">as.character</span>(experiment_design.ord$SampleID)
<span class="co"># create factors for future plotting</span>
group<-<span class="kw">factor</span>(experiment_design.ord$tissue)
group</code></pre>
<pre><code>## [1] brain brain liver liver liver brain brain liver
## Levels: brain liver</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">age<-<span class="kw">factor</span>(experiment_design.ord$age)
age</code></pre>
<pre><code>## [1] 26 26 30 30 30 26 26 30
## Levels: 26 30</code></pre>
<h4 id="basic-qc-plots">Basic QC plots</h4>
<p>Density plots of log-intensity distribution of each library can be superposed on a single graph for a better comparison between libraries and for identification of libraries with weird distribution. On the boxplots the density distributions of raw log-intensities are not expected to be identical but still not totally different.</p>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pushpin"></span>NOTE:</h2>
</div>
<div class="panel-body">
<p>The <em>png</em> function will create a png file where to save the plots created straight after, and will close this file when <em>dev.off()</em> is called.</p>
<p>To see your plots interactively, simply ommit those two lines.</p>
</div>
</aside>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># density plot of raw read counts (log10)</span>
<span class="kw">png</span>(<span class="dt">file=</span><span class="st">"~/Raw_read_counts_per_gene.density.png"</span>)
logcounts <-<span class="st"> </span><span class="kw">log</span>(counts$counts[,<span class="dv">1</span>],<span class="dv">10</span>)
d <-<span class="st"> </span><span class="kw">density</span>(logcounts)
<span class="kw">plot</span>(d,<span class="dt">xlim=</span><span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">8</span>),<span class="dt">main=</span><span class="st">""</span>,<span class="dt">ylim=</span><span class="kw">c</span>(<span class="dv">0</span>,.<span class="dv">45</span>),<span class="dt">xlab=</span><span class="st">"Raw read counts per gene (log10)"</span>, <span class="dt">ylab=</span><span class="st">"Density"</span>)
for (s in <span class="dv">2</span>:<span class="kw">length</span>(samples)){
logcounts <-<span class="st"> </span><span class="kw">log</span>(counts$counts[,s],<span class="dv">10</span>)
d <-<span class="st"> </span><span class="kw">density</span>(logcounts)
<span class="kw">lines</span>(d)
}
<span class="kw">dev.off</span>()</code></pre>
<p><img src="fig/RNAseq/Raw_read_counts_per_gene.density.png" alt="boxplotlog10" /></p>
<p>Boxplots of the raw read counts after log10 transformation.</p>
<pre class="sourceCode r"><code class="sourceCode r">## box plots of raw read counts (log10)
<span class="kw">png</span>(<span class="dt">file=</span><span class="st">"~/Raw_read_counts_per_gene.boxplot.png"</span>)
logcounts <-<span class="st"> </span><span class="kw">log</span>(counts$counts,<span class="dv">10</span>)
<span class="kw">boxplot</span>(logcounts, <span class="dt">main=</span><span class="st">""</span>, <span class="dt">xlab=</span><span class="st">""</span>, <span class="dt">ylab=</span><span class="st">"Raw read counts per gene (log10)"</span>,<span class="dt">axes=</span><span class="ot">FALSE</span>)
<span class="kw">axis</span>(<span class="dv">2</span>)
<span class="kw">axis</span>(<span class="dv">1</span>,<span class="dt">at=</span><span class="kw">c</span>(<span class="dv">1</span>:<span class="kw">length</span>(samples)),<span class="dt">labels=</span><span class="kw">colnames</span>(logcounts),<span class="dt">las=</span><span class="dv">2</span>,<span class="dt">cex.axis=</span><span class="fl">0.8</span>)
<span class="kw">dev.off</span>()</code></pre>
<p><img src="fig/RNAseq/Raw_read_counts_per_gene.boxplot.png" alt="boxplotlog10" /></p>
<p>In order to investigate the relationship between samples, hierarchical clustering can be performed using the <em>heatmap</em> function from the <em>stats</em> package. In this example <em>heatmap</em> calculates a matrix of euclidean distances from the mapped read counts for the 100 most highly expressed genes.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># select data for the 100 most highly expressed genes</span>
select =<span class="st"> </span><span class="kw">order</span>(<span class="kw">rowMeans</span>(counts$counts), <span class="dt">decreasing=</span><span class="ot">TRUE</span>)[<span class="dv">1</span>:<span class="dv">100</span>]
highexprgenes_counts <-<span class="st"> </span>counts$counts[select,]
<span class="co"># heatmap with sample name on X-axis</span>
<span class="kw">png</span>(<span class="dt">file=</span><span class="st">"~/High_expr_genes.heatmap.png"</span>)
<span class="kw">heatmap</span>(highexprgenes_counts, <span class="dt">col=</span><span class="kw">topo.colors</span>(<span class="dv">50</span>), <span class="dt">margin=</span><span class="kw">c</span>(<span class="dv">10</span>,<span class="dv">6</span>))
<span class="kw">dev.off</span>()</code></pre>
<p>You will notice that the samples clustering does not follow the original order in the data matrix (alphabetical order “ERR420386” to “ERR420393”). They have been re-ordered according to the similarity of the 100 genes expression profiles. To understand what biological effect lies under this clustering, one can use the samples annotation for labeling (samples group, age, sex etc).</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># heatmap with condition group as labels</span>
<span class="kw">colnames</span>(highexprgenes_counts)<-<span class="st"> </span>group
<span class="co"># plot</span>
<span class="kw">png</span>(<span class="dt">file=</span><span class="st">"~/High_exprs_genes.heatmap.group.png"</span>)
<span class="kw">heatmap</span>(highexprgenes_counts, <span class="dt">col =</span> <span class="kw">topo.colors</span>(<span class="dv">50</span>), <span class="dt">margin=</span><span class="kw">c</span>(<span class="dv">10</span>,<span class="dv">6</span>))
<span class="kw">dev.off</span>()</code></pre>
<p><img src="fig/RNAseq/High_expr_genes.heatmap.group.png" alt="heatmap" /></p>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pencil"></span>Exercise: Heatmap</h2>
</div>
<div class="panel-body">
<p>Produce a heatmap for the 50 most highly expressed genes and annotate the samples with with their age.</p>
<ul>
<li>Subset the read counts object for the 50 most highly expressed genes</li>
<li>Annotate the samples in the subset with their age (check order with design!)</li>
<li>Plot a heatmap with this subset of data, scaling genes and ordering both genes and samples</li>
</ul>
</div>
</section>
<h4 id="principal-component-analysis">Principal Component Analysis</h4>
<p>A Principal Component Analysis (PCA) can also be performed with these data using the <em>cmdscale</em> function (from the <em>stats</em> package) which performs a classical multidimensional scaling of a data matrix.</p>
<p>Reads counts need to be transposed before being analysed with the <em>cmdscale</em> functions, i.e. genes should be in columns and samples should be in rows. This is the code for transposing and checking the data before further steps:</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># select data for the 1000 most highly expressed genes</span>
select =<span class="st"> </span><span class="kw">order</span>(<span class="kw">rowMeans</span>(counts$counts), <span class="dt">decreasing=</span><span class="ot">TRUE</span>)[<span class="dv">1</span>:<span class="dv">100</span>]
highexprgenes_counts <-<span class="st"> </span>counts$counts[select,]
<span class="co"># annotate the data with condition group as labels</span>
<span class="kw">colnames</span>(highexprgenes_counts)<-<span class="st"> </span>group
<span class="co"># transpose the data to have variables (genes) as columns</span>
data_for_PCA <-<span class="st"> </span><span class="kw">t</span>(highexprgenes_counts)
<span class="kw">dim</span>(data_for_PCA)</code></pre>
<pre><code>## [1] 8 100</code></pre>
<p>The <em>cmdscale</em> function will calculate a matrix of dissimilarities from your transposed data and will also provide information about the proportion of explained variance by calculating Eigen values.</p>
<pre class="sourceCode r"><code class="sourceCode r">## calculate MDS (matrix of dissimilarities)
mds <-<span class="st"> </span><span class="kw">cmdscale</span>(<span class="kw">dist</span>(data_for_PCA), <span class="dt">k=</span><span class="dv">3</span>, <span class="dt">eig=</span><span class="ot">TRUE</span>)
<span class="co"># k = the maximum dimension of the space which the data are to be represented in</span>
<span class="co"># eig = indicates whether eigenvalues should be returned</span></code></pre>
<p>The variable mds$eig provides the Eigen values for the first 8 principal components:</p>
<pre class="sourceCode r"><code class="sourceCode r">mds$eig</code></pre>
<pre><code>## [1] 9.490938e+13 1.099639e+13 1.125271e+11 1.026586e+10 1.500500e+07
## [6] 6.240239e+06 3.206875e+06 2.285361e-03</code></pre>
<p>Plotting this variable as a percentage will help you determine how many components can explain the variability in your dataset and thus how many dimensions you should be looking at.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># transform the Eigen values into percentage</span>
eig_pc <-<span class="st"> </span>mds$eig *<span class="st"> </span><span class="dv">100</span> /<span class="st"> </span><span class="kw">max</span>(mds$eig)
<span class="co"># plot the PCA</span>
<span class="kw">png</span>(<span class="dt">file=</span><span class="st">"~/PCA_PropExplainedVariance.png"</span>)
<span class="kw">plot</span>(eig_pc,
<span class="dt">type=</span><span class="st">"h"</span>, <span class="dt">lwd=</span><span class="dv">15</span>, <span class="dt">las=</span><span class="dv">1</span>,
<span class="dt">xlab=</span><span class="st">"Dimensions"</span>,
<span class="dt">ylab=</span><span class="st">"Proportion of explained variance"</span>, <span class="dt">y.axis=</span><span class="ot">NULL</span>,
<span class="dt">col=</span><span class="st">"darkgrey"</span>)
<span class="kw">dev.off</span>()</code></pre>
<p><img src="fig/RNAseq/PCA_PropExplainedVariance.png" alt="prop" /></p>
<p>In most cases, the first 2 components explain more than half the variability in the dataset and can be used for plotting. The <em>cmdscale</em> function run with default parameters will perform a principal components analysis on the given data matrix and the <em>plot</em> function will provide scatter plots for individuals representation.</p>
<pre class="sourceCode r"><code class="sourceCode r">## calculate MDS
mds <-<span class="st"> </span><span class="kw">cmdscale</span>(<span class="kw">dist</span>(data_for_PCA)) <span class="co"># Performs MDS analysis </span></code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co">#Samples representation</span>
<span class="kw">png</span>(<span class="dt">file=</span><span class="st">"~/PCA_Dim1vsDim2.png"</span>)
<span class="kw">plot</span>(mds[,<span class="dv">1</span>], -mds[,<span class="dv">2</span>], <span class="dt">type=</span><span class="st">"n"</span>, <span class="dt">xlab=</span><span class="st">"Dimension 1"</span>, <span class="dt">ylab=</span><span class="st">"Dimension 2"</span>, <span class="dt">main=</span><span class="st">""</span>)
<span class="kw">text</span>(mds[,<span class="dv">1</span>], -mds[,<span class="dv">2</span>], <span class="kw">rownames</span>(mds), <span class="dt">cex=</span><span class="fl">0.8</span>)
<span class="kw">dev.off</span>()</code></pre>
<p><img src="fig/RNAseq/PCA_Dim1vsDim2.group.png" alt="PCA1" /></p>
<p>The PCA plot of the first two components show a clear separation of the Brain and Liver samples across the 1st dimension. Within each sample group we can also notice a split between the 4 samples of each group, which seem to cluster in pair. This observation can be explained by another factor of variability in the data, commonly batch effect or another biological biais such as age or sex.</p>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pencil"></span>Exercise: PCA</h2>
</div>
<div class="panel-body">
<p>Produce a PCA plot from the read counts of the 500 most highly expressed genes and change the labels until you can identify the reason for the split between samples from the same tissue.</p>
<ul>
<li>Get the read counts for the 500 most highly expressed genes</li>
<li>Transpose this matrix of read counts</li>
<li>Check the number of dimensions explaining the variability in the dataset</li>
<li>Run the PCA with an appropriate number of components</li>
<li>Annotate the samples with their age & re-run the PCA & plot the main components</li>
<li>Annotate the samples with other clinical data & re-run the PCA & plot the main components until you can separate the samples within each tissue group</li>
</ul>
</div>
</section>
<h3 id="differential-expression">Differential Expression</h3>
<p>Before proceeding with differential expression analysis, it is useful to filter out very lowly expressed genes. This will help increasing the statistical power of the analysi while keeping genes of interest. A common way to do this is by filtering out genes having less than 1 count-per-million reads (cpm) in half the samples. The “edgeR” library provides the <em>cpm</em> function which can be used here.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># load required libraries</span>
<span class="kw">library</span>(edgeR)</code></pre>
<pre><code>## Loading required package: limma</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># get the expression counts from previous alignment step</span>
mycounts <-<span class="st"> </span>counts$counts
<span class="kw">dim</span>(mycounts)</code></pre>
<pre><code>## [1] 25702 8</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">mycounts[<span class="dv">1</span>:<span class="dv">5</span>,<span class="dv">1</span>:<span class="dv">3</span>]</code></pre>
<pre><code>## ERR420386 ERR420387 ERR420388
## 1 2 56 5461
## 2 3723 2270 27665
## 9 14 3 148
## 10 1 1 373
## 12 42 34 20969</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># filtering</span>
<span class="co">#Keep genes with least 1 count-per-million reads (cpm) in at least 4 samples</span>
isexpr <-<span class="st"> </span><span class="kw">rowSums</span>(<span class="kw">cpm</span>(mycounts)><span class="dv">1</span>) >=<span class="st"> </span><span class="dv">4</span>
<span class="kw">table</span>(isexpr)</code></pre>
<pre><code>## isexpr
## FALSE TRUE
## 10686 15016</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">mycounts <-<span class="st"> </span>mycounts[isexpr,]
genes <-<span class="st"> </span><span class="kw">rownames</span>(mycounts)
<span class="kw">dim</span>(mycounts)</code></pre>
<pre><code>## [1] 15016 8</code></pre>
<p>The <em>limma</em> package (since version 3.16.0) offers the <em>voom</em> function that will normalise read counts and apply a linear model to the normalised data before computing moderated t-statistics of differential expression.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># load required libraries</span>
<span class="kw">library</span>(limma)
<span class="co"># check your samples grouping</span>
experiment_design.ord[<span class="kw">colnames</span>(mycounts),]$tissue ==<span class="st"> </span>group</code></pre>
<pre><code>## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># create design matrix for limma</span>
design <-<span class="st"> </span><span class="kw">model.matrix</span>(~<span class="dv">0</span>+group)
<span class="co"># substitute "group" from the design column names</span>
<span class="kw">colnames</span>(design)<-<span class="st"> </span><span class="kw">gsub</span>(<span class="st">"group"</span>,<span class="st">""</span>,<span class="kw">colnames</span>(design))
<span class="co"># check your design matrix</span>
design</code></pre>
<pre><code>## brain liver
## 1 1 0
## 2 1 0
## 3 0 1
## 4 0 1
## 5 0 1
## 6 1 0
## 7 1 0
## 8 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># calculate normalization factors between libraries</span>
nf <-<span class="st"> </span><span class="kw">calcNormFactors</span>(mycounts)
<span class="co">#</span>
<span class="co"># normalise the read counts with 'voom' function</span>
y <-<span class="st"> </span><span class="kw">voom</span>(mycounts,design,<span class="dt">lib.size=</span><span class="kw">colSums</span>(mycounts)*nf)
<span class="co"># extract the normalised read counts</span>
counts.voom <-<span class="st"> </span>y$E
<span class="co"># save normalised expression data into output dir</span>
<span class="kw">write.table</span>(counts.voom,<span class="dt">file=</span><span class="st">"~/counts.voom.txt"</span>,<span class="dt">row.names=</span>T,<span class="dt">quote=</span>F,<span class="dt">sep=</span><span class="st">"</span><span class="ch">\t</span><span class="st">"</span>)
<span class="co"># fit linear model for each gene given a series of libraries</span>
fit <-<span class="st"> </span><span class="kw">lmFit</span>(y,design)
<span class="co"># construct the contrast matrix corresponding to specified contrasts of a set of parameters</span>
cont.matrix <-<span class="st"> </span><span class="kw">makeContrasts</span>(liver-brain,<span class="dt">levels=</span>design)
cont.matrix </code></pre>
<pre><code>## Contrasts
## Levels liver - brain
## brain -1
## liver 1</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># compute estimated coefficients and standard errors for a given set of contrasts</span>
fit <-<span class="st"> </span><span class="kw">contrasts.fit</span>(fit, cont.matrix)
<span class="co"># compute moderated t-statistics of differential expression by empirical Bayes moderation of the standard errors</span>
fit <-<span class="st"> </span><span class="kw">eBayes</span>(fit)
<span class="kw">options</span>(<span class="dt">digits=</span><span class="dv">3</span>)
<span class="co"># check the output fit</span>
<span class="kw">dim</span>(fit)</code></pre>
<pre><code>## [1] 15016 1</code></pre>
<p>The <em>topTable</em> function summarises the output from limma in a table format. Significant DE genes for a particular comparison can be identified by selecting genes with a p-value smaller than a chosen cut-off value and/or a fold change greater than a chosen value in this table. By default the table will be sorted by increasing adjusted p-value, showing the most significant DE genes at the top.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># set adjusted pvalue threshold and log fold change threshold</span>
mypval=<span class="fl">0.01</span>
myfc=<span class="dv">3</span>
<span class="co"># get the coefficient name for the comparison of interest</span>
<span class="kw">colnames</span>(fit$coefficients)</code></pre>
<pre><code>## [1] "liver - brain"</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">mycoef=<span class="st">"liver - brain"</span>
<span class="co"># get the output table for the 10 most significant DE genes for this comparison</span>
<span class="kw">topTable</span>(fit,<span class="dt">coef=</span>mycoef)</code></pre>
<pre><code>## logFC AveExpr t P.Value adj.P.Val B
## 5265 10.45 7.33 86.7 1.89e-13 1.37e-09 18.6
## 3240 11.28 7.62 89.4 1.46e-13 1.37e-09 18.5
## 2335 6.82 8.56 67.2 1.53e-12 2.93e-09 18.4
## 213 11.77 10.69 69.3 1.19e-12 2.93e-09 18.3
## 338 11.51 8.31 70.4 1.04e-12 2.93e-09 18.0
## 2243 11.56 7.17 80.0 3.64e-13 1.37e-09 17.7
## 716 7.49 6.75 61.0 3.38e-12 4.22e-09 17.7
## 229 10.90 6.45 82.0 2.98e-13 1.37e-09 17.4
## 1571 9.64 6.62 63.8 2.33e-12 3.18e-09 17.3
## 125 11.54 7.08 64.3 2.19e-12 3.18e-09 16.9</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># get the full table ("n = number of genes in the fit")</span>
limma.res <-<span class="st"> </span><span class="kw">topTable</span>(fit,<span class="dt">coef=</span>mycoef,<span class="dt">n=</span><span class="kw">dim</span>(fit)[<span class="dv">1</span>])
<span class="co"># get significant DE genes only (adjusted p-value < mypval)</span>
limma.res.pval <-<span class="st"> </span><span class="kw">topTable</span>(fit,<span class="dt">coef=</span>mycoef,<span class="dt">n=</span><span class="kw">dim</span>(fit)[<span class="dv">1</span>],<span class="dt">p.val=</span>mypval)
<span class="kw">dim</span>(limma.res.pval)</code></pre>
<pre><code>## [1] 8183 6</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># get significant DE genes with low adjusted p-value high fold change</span>
limma.res.pval.FC <-<span class="st"> </span>limma.res.pval[<span class="kw">which</span>(<span class="kw">abs</span>(limma.res.pval$logFC)>myfc),]
<span class="kw">dim</span>(limma.res.pval.FC)</code></pre>
<pre><code>## [1] 3044 6</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># write limma output table for significant genes into a tab delimited file</span>
filename =<span class="st"> </span><span class="kw">paste</span>(<span class="st">"~/DEgenes_"</span>,mycoef,<span class="st">"_pval"</span>,mypval,<span class="st">"_logFC"</span>,myfc,<span class="st">".txt"</span>,<span class="dt">sep=</span><span class="st">""</span>)
<span class="kw">write.table</span>(limma.res.pval.FC,<span class="dt">file=</span>filename,<span class="dt">row.names=</span>T,<span class="dt">quote=</span>F,<span class="dt">sep=</span><span class="st">"</span><span class="ch">\t</span><span class="st">"</span>)</code></pre>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pencil"></span>Exercise: Limma</h2>
</div>
<div class="panel-body">
<p>Get the number of DE genes between technical group 1 and technical group 2 (all Brain samples) with adj pvalue<0.01.</p>
<ul>
<li>Create a new design matrix for limma with the technical replicate groups</li>
<li>Re-normalise the read counts with ‘voom’ function with new design matrix</li>
<li>Fit a linear model on these normalised data</li>
<li>Make the contrast matrix corresponding to the new set of parameters</li>
<li>Fit the contrast matrix to the linear model</li>
<li>Compute moderated t-statistics of differential expression</li>
<li>Get the output table for the 10 most significant DE genes for this comparison</li>
</ul>
</div>
</section>
<h3 id="gene-annotation">Gene Annotation</h3>
<p>The annotation of EntrezGene IDs from RNAseq data can be done using the BioMart database which contains many species including Human, Mouse, Zebrafish, Chicken and Rat.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># get the Ensembl annotation for human genome</span>
<span class="kw">library</span>(biomaRt)</code></pre>
<pre><code>## Loading required package: methods</code></pre>
<pre class="sourceCode r"><code class="sourceCode r">mart<-<span class="st"> </span><span class="kw">useDataset</span>(<span class="st">"hsapiens_gene_ensembl"</span>, <span class="kw">useMart</span>(<span class="st">"ENSEMBL_MART_ENSEMBL"</span>,<span class="dt">host=</span><span class="st">"www.ensembl.org"</span>))
<span class="co"># get entrez gene IDs from limma output table</span>
entrez_genes <-<span class="st"> </span><span class="kw">as.character</span>(<span class="kw">rownames</span>(limma.res.pval.FC))
<span class="kw">length</span>(entrez_genes)</code></pre>
<pre><code>## [1] 3044</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># interrogate the BioMart database to get gene symbol and description for these genes </span>
detags.IDs <-<span class="st"> </span><span class="kw">getBM</span>(
<span class="dt">filters=</span> <span class="st">"entrezgene"</span>,
<span class="dt">attributes=</span> <span class="kw">c</span>(<span class="st">"entrezgene"</span>,<span class="st">"hgnc_symbol"</span>,<span class="st">"description"</span>),
<span class="dt">values=</span> entrez_genes,
<span class="dt">mart=</span> mart
)
<span class="kw">dim</span>(detags.IDs)</code></pre>
<pre><code>## [1] 2880 3</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">head</span>(detags.IDs)</code></pre>
<pre><code>## entrezgene hgnc_symbol
## 1 10 NAT2
## 2 10000 AKT3
## 3 100033415 SNORD116-3
## 4 100033421 SNORD116-9
## 5 100033427 SNORD116-15
## 6 100033442 SNORD115-5
## description
## 1 N-acetyltransferase 2 (arylamine N-acetyltransferase) [Source:HGNC Symbol;Acc:HGNC:7646]
## 2 AKT serine/threonine kinase 3 [Source:HGNC Symbol;Acc:HGNC:393]
## 3 small nucleolar RNA, C/D box 116-3 [Source:HGNC Symbol;Acc:HGNC:33069]
## 4 small nucleolar RNA, C/D box 116-9 [Source:HGNC Symbol;Acc:HGNC:33075]
## 5 small nucleolar RNA, C/D box 116-15 [Source:HGNC Symbol;Acc:HGNC:33081]
## 6 small nucleolar RNA, C/D box 115-5 [Source:HGNC Symbol;Acc:HGNC:33024]</code></pre>
<p>In many cases, several annotations ar available per entrez gene ID. This results in duplicate entries in the output table from <em>getBM</em>. The simplest way to deal with this issue is to remove duplicates, although they can also be concatenated in some ways.</p>
<p>Once the annotation has been obtained for all DE genes, this table can be merged with the output table from limma for a complete result and an easier interpretation.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># remove duplicates</span>
detags.IDs.matrix<-detags.IDs[-<span class="kw">which</span>(<span class="kw">duplicated</span>(detags.IDs$entrezgene)),]
<span class="co"># select genes of interest only</span>
<span class="kw">rownames</span>(detags.IDs.matrix)<-detags.IDs.matrix$entrezgene
entrez_genes.annot <-<span class="st"> </span>detags.IDs.matrix[<span class="kw">as.character</span>(entrez_genes),]
<span class="co"># join the two tables</span>
<span class="kw">rownames</span>(limma.res.pval.FC) <-<span class="st"> </span>limma.res.pval.FC$ID
limma.res.pval.FC.annot <-<span class="st"> </span><span class="kw">cbind</span>(entrez_genes.annot,limma.res.pval.FC)
<span class="co"># check the annotated table</span>
<span class="kw">head</span>(limma.res.pval.FC.annot)</code></pre>
<pre><code>## entrezgene hgnc_symbol
## 5265 5265 SERPINA1
## 3240 3240 HP
## 2335 2335 FN1
## 213 213 ALB
## 338 338 APOB
## 2243 2243 FGA
## description logFC
## 5265 serpin family A member 1 [Source:HGNC Symbol;Acc:HGNC:8941] 10.45
## 3240 haptoglobin [Source:HGNC Symbol;Acc:HGNC:5141] 11.28
## 2335 fibronectin 1 [Source:HGNC Symbol;Acc:HGNC:3778] 6.82
## 213 albumin [Source:HGNC Symbol;Acc:HGNC:399] 11.77
## 338 apolipoprotein B [Source:HGNC Symbol;Acc:HGNC:603] 11.51
## 2243 fibrinogen alpha chain [Source:HGNC Symbol;Acc:HGNC:3661] 11.56
## AveExpr t P.Value adj.P.Val B
## 5265 7.33 86.7 1.89e-13 1.37e-09 18.6
## 3240 7.62 89.4 1.46e-13 1.37e-09 18.5
## 2335 8.56 67.2 1.53e-12 2.93e-09 18.4
## 213 10.69 69.3 1.19e-12 2.93e-09 18.3
## 338 8.31 70.4 1.04e-12 2.93e-09 18.0
## 2243 7.17 80.0 3.64e-13 1.37e-09 17.7</code></pre>
<h3 id="gene-set-enrichment">Gene Set Enrichment</h3>
<p>Gene Ontology (GO) enrichment is a method for investigating sets of genes using the Gene Ontology system of classification, in which genes are assigned to a particular set of terms for three major domains: cellular component, biological process and molecular function.</p>
<p>The <em>GOstats</em> package can test for both over and under representation of GO terms using the Hypergeometric test. The output of the analysis is typically a ranked list of GO terms, each associated with a p-value.</p>
<p>The Hypergeometric test will require both a list of selected genes (i.e. your DE genes) and a “universe” list (e.g. all genes annotated in the genome you are working with), all represented by their “EntrezGene” ID.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># load the library</span>
<span class="kw">library</span>(GOstats)</code></pre>
<pre><code>## Loading required package: Biobase</code></pre>
<pre><code>## Loading required package: BiocGenerics</code></pre>
<pre><code>## Loading required package: parallel</code></pre>
<pre><code>##
## Attaching package: 'BiocGenerics'</code></pre>
<pre><code>## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB</code></pre>
<pre><code>## The following object is masked from 'package:limma':
##
## plotMA</code></pre>
<pre><code>## The following objects are masked from 'package:stats':
##
## IQR, mad, xtabs</code></pre>
<pre><code>## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, cbind, colnames,
## do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, lengths, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff,
## sort, table, tapply, union, unique, unsplit</code></pre>
<pre><code>## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.</code></pre>
<pre><code>## Loading required package: Category</code></pre>
<pre><code>## Loading required package: stats4</code></pre>
<pre><code>## Loading required package: AnnotationDbi</code></pre>
<pre><code>## Loading required package: IRanges</code></pre>
<pre><code>## Loading required package: S4Vectors</code></pre>
<pre><code>##
## Attaching package: 'S4Vectors'</code></pre>
<pre><code>## The following objects are masked from 'package:base':
##
## colMeans, colSums, expand.grid, rowMeans, rowSums</code></pre>
<pre><code>## Loading required package: Matrix</code></pre>
<pre><code>##
## Attaching package: 'Matrix'</code></pre>
<pre><code>## The following object is masked from 'package:S4Vectors':
##
## expand</code></pre>
<pre><code>## Loading required package: graph</code></pre>
<pre><code>## </code></pre>
<pre><code>##
## Attaching package: 'GOstats'</code></pre>
<pre><code>## The following object is masked from 'package:AnnotationDbi':
##
## makeGOGraph</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Define list of genes of interest (DE genes - EntrezGene IDs)</span>
entrezgeneids <-<span class="st"> </span><span class="kw">as.character</span>(<span class="kw">rownames</span>(limma.res.pval.FC))
<span class="kw">length</span>(entrezgeneids)</code></pre>
<pre><code>## [1] 3044</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Define the universe</span>
universeids <-<span class="st"> </span><span class="kw">rownames</span>(mycounts)
<span class="kw">length</span>(universeids)</code></pre>
<pre><code>## [1] 15016</code></pre>
<p>Before running the hypergeometric test with the <em>hyperGTest</em> function, you need to define the parameters for the test (gene lists, ontology, test direction) as well as the annotation database to be used. The ontology to be tested can be any of the three GO domains: biological process (“BP”), cellular component (“CC”) or molecular function (“MF”).</p>
<p>In the example below we will test for over-represented biological processes in our list of differentially expressed genes.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># define the p-value cut off for the hypergeometric test</span>
hgCutoff <-<span class="st"> </span><span class="fl">0.05</span>
params <-<span class="st"> </span><span class="kw">new</span>(<span class="st">"GOHyperGParams"</span>,<span class="dt">annotation=</span><span class="st">"org.Hs.eg"</span>,<span class="dt">geneIds=</span>entrezgeneids,<span class="dt">universeGeneIds=</span>universeids,<span class="dt">ontology=</span><span class="st">"BP"</span>,<span class="dt">pvalueCutoff=</span>hgCutoff,<span class="dt">testDirection=</span><span class="st">"over"</span>)</code></pre>
<pre><code>## Loading required package: org.Hs.eg.db</code></pre>
<pre><code>## </code></pre>
<pre><code>## Warning in makeValidParams(.Object): removing geneIds not in
## universeGeneIds</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Run the test</span>
hg <-<span class="st"> </span><span class="kw">hyperGTest</span>(params)
<span class="co"># Check results</span>
hg</code></pre>
<pre><code>## Gene to GO BP test for over-representation
## 9093 GO BP ids tested (3464 have p < 0.05)
## Selected gene set size: 1528
## Gene universe size: 12123
## Annotation package: org.Hs.eg</code></pre>
<p>You can get the output table from the test for significant GO terms only by adjusting the pvalues with the <em>p.adjust</em> function.</p>
<pre class="sourceCode r"><code class="sourceCode r">## Get the p-values of the test
hg.pv <-<span class="st"> </span><span class="kw">pvalues</span>(hg)
## Adjust p-values for multiple test (FDR)
hg.pv.fdr <-<span class="st"> </span><span class="kw">p.adjust</span>(hg.pv,<span class="st">'fdr'</span>)
## select the GO terms with adjusted p-value less than the cut off
sigGO.ID <-<span class="st"> </span><span class="kw">names</span>(hg.pv.fdr[hg.pv.fdr <<span class="st"> </span>hgCutoff])
<span class="kw">length</span>(sigGO.ID)</code></pre>
<pre><code>## [1] 2527</code></pre>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># get table from HyperG test result</span>
df <-<span class="st"> </span><span class="kw">summary</span>(hg)
<span class="co"># keep only significant GO terms in the table</span>
GOannot.table <-<span class="st"> </span>df[df[,<span class="dv">1</span>] %in%<span class="st"> </span>sigGO.ID,]
<span class="kw">head</span>(GOannot.table)</code></pre>
<pre><code>## GOBPID Pvalue OddsRatio ExpCount Count Size
## 1 GO:0042221 1.56e-104 3.45 395 762 3131
## 2 GO:0050896 1.79e-86 3.10 753 1108 5977
## 3 GO:0044707 5.56e-69 2.64 562 877 4458
## 4 GO:0032501 2.60e-68 2.62 599 915 4756
## 5 GO:0010033 1.11e-65 2.85 298 563 2362
## 6 GO:0044699 2.31e-63 5.07 1242 1453 9851
## Term
## 1 response to chemical
## 2 response to stimulus
## 3 single-multicellular organism process
## 4 multicellular organismal process
## 5 response to organic substance
## 6 single-organism process</code></pre>
<p>The Gene Ontology enrichment result can be saved in a text file or an html file for future reference.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Create text report of the significantly over-represented GO terms</span>
<span class="kw">write.table</span>(GOannot.table,<span class="dt">file=</span><span class="st">"~/GOterms_OverRep_BP.txt"</span>,<span class="dt">sep=</span><span class="st">"</span><span class="ch">\t</span><span class="st">"</span>,<span class="dt">row.names=</span>F)
<span class="co"># Create html report of all over-represented GO terms</span>
<span class="kw">htmlReport</span>(hg, <span class="dt">file=</span><span class="st">"~/GOterms_OverRep_BP.html"</span>)</code></pre>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2><span class="glyphicon glyphicon-pencil"></span>Exercise: GOstats</h2>
</div>
<div class="panel-body">
<p>Identify the GO terms in the Molecular Function domain that are over-represented (pvalue<0.01) in your list of DE genes.</p>
<ul>
<li>Get your list of DE genes (Entrez Gene IDs)</li>
<li>Set the new parameters for the hypergeometric test</li>
<li>Run the test and adjust the pvalues in the output object</li>
<li>Identify the significant GO terms at pvalue 0.01</li>
</ul>
</div>
</section>
<p>Other softwares can be used to investigate over-represented pathways, such as GeneGO (<a href="https://portal.genego.com/">https://portal.genego.com/</a>) and Ingenuity (<a href="http://www.ingenuity.com/products/ipa">http://www.ingenuity.com/products/ipa</a>). The advantage of these softwares is that they maintain curated and up-to-date extensive databases. They also provide intuitive visualisation and network modelling tools.</p>
<p>You can save an image of your RNAseq analysis before moving on to the next part of this course.</p>
<pre class="sourceCode r"><code class="sourceCode r">RDataFile <-<span class="st"> "~/RNAseq_DE_analysis_with_R.RData"</span>
<span class="kw">save.image</span>(RDataFile)</code></pre>
<h3 id="record-package-and-version-info-with-sessioninfo">Record package and version info with <code>sessionInfo()</code></h3>
<p>The <em>sessionInfo</em> function prints version information about R and any attached packages. It’s a good practice to always run this command at the end of your R session and record it for the sake of reproducibility in the future.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">sessionInfo</span>()</code></pre>
<pre><code>## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.4 LTS
##
## locale:
## [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
## [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
## [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel methods stats graphics grDevices utils
## [8] datasets base
##
## other attached packages:
## [1] GO.db_3.3.0 org.Hs.eg.db_3.3.0 GOstats_2.38.1
## [4] graph_1.50.0 Category_2.38.0 Matrix_1.2-6
## [7] AnnotationDbi_1.34.4 IRanges_2.6.1 S4Vectors_0.10.2
## [10] Biobase_2.32.0 BiocGenerics_0.18.0 biomaRt_2.28.0
## [13] edgeR_3.14.0 limma_3.28.17 Rsubread_1.22.2
##
## loaded via a namespace (and not attached):
## [1] knitr_1.13 magrittr_1.5 splines_3.3.1
## [4] xtable_1.8-2 lattice_0.20-33 stringr_1.0.0
## [7] tools_3.3.1 grid_3.3.1 AnnotationForge_1.14.2
## [10] DBI_0.4-1 genefilter_1.54.2 survival_2.39-5
## [13] RBGL_1.48.1 GSEABase_1.34.0 formatR_1.4
## [16] bitops_1.0-6 RCurl_1.95-4.8 RSQLite_1.0.0
## [19] evaluate_0.9 stringi_1.1.1 XML_3.98-1.4
## [22] annotate_1.50.0</code></pre>
</div>
</div>
</article>
<div class="footer">
<a class="label swc-blue-bg" href="https://github.com/MonashBioinformaticsPlatform/RNAseq-DE-analysis-with-R">Source</a>
</div>
</div>
<!-- Javascript placed at the end of the document so the pages load faster -->
<script src="http://code.jquery.com/jquery-1.9.1.js"></script>
<script src="css/bootstrap/bootstrap-js/bootstrap.js"></script>
</body>
</html>