Plotting Exercise 5

These exercises cover the sections on Plotting in R PlottingInR.

Section 1 - barplot

The paired-end RNAseq data for our case study was aligned to mouse genome (mm9) using Tophat2. The alignment summary can be found here "data/tophat_alignment_overallstat.txt" with the following columns.

• sample: sample name
• num.input.pairs: number of original pairs
• num.aligned.pairs: number of aligned pairs
• num.uniquely_mapped.pairs: number of properly aligned paired-reads, uniquely mapped
• num.multiple_mapped.pairs: number of properly aligned paired-reads, multiple mapped
• num.discordant_mapped.pairs: number of discordant aligned paired-reads

num.aligned.pairs = num.uniquely_mapped.pairs + num.multiple_mapped.pairs + num.discordant_mapped.pairs

• hint: use read.table() function to load tophat_alignment_overallstat.txt

1. load tophat_alignment_overallstat.txt to R as a data.frame named aligned_res

head(aligned_res)

##     sample num.input.pairs num.aligned.pairs num.uniquely_mapped.pairs
## 1 DOX24H_1        41113500          30527725                  29078260
## 2 DOX24H_2        38166396          29575699                  28188066
## 3 DOX24H_3        42506262          31376922                  29572848
## 4  DOX7D_1        46156954          34147421                  32392402
## 5  DOX7D_2        40372425          32876374                  31409139
## 6  DOX7D_3        46160430          34722360                  32896496
##   num.multiple_mapped.pairs num.discordant_mapped.pairs
## 1                   1249793                      199672
## 2                   1173844                      213789
## 3                   1553296                      250778
## 4                   1551960                      203059
## 5                   1201468                      265767
## 6                   1615952                      209912

2. add a new column 'num.unaligned.pairs' to aligned_res

3. create the bar plot below using geom_bar(stat = "identity",position=position_dodge())

• hints:

  • you need to convert the ‘wide’ data.frame to a ‘long’ data.frame using melt() function from reshape2 package

  • check ?geom_bar and see description for position, and see what happens after you put position=position_dodge() in geom_bar()

  • The x-axis labels can be rotated 45 degrees by adding layer theme(axis.text.x = element_text(angle=45, hjust=1))

## Warning in register(): Can't find generic `scale_type` in package ggplot2 to
## register S3 method.

4. create the bar plot below: using geom_bar(position = "fill")

• hint: you need to understand how factor data type works

5. create the bar plot below: using facet_grid()

• hint: check the Examples from ?facet_grid

6. create the bar plot below: using facet_wrap()

• hint: check the Examples from ?facet_wrap

Section 2 - bubble plot [bonus exercise]

Gene Set Enrichment Analysis (GSEA) is one of the popular functional analysis tools for bulkRNAseq data. More details please see [https://www.gsea-msigdb.org/gsea/index.jsp]. File "data/GSEA_hallmark_Dox7D_Minus_Dox1D.csv" summarizes the GSEA results for bulkRNAseq data based on Dox_7D vs Dox_24h using the hallmark gene sets.

GSEA_hallmark_Dox7D_Minus_Dox1D.csv contains the following columns:

• NAME: full Gene set name
• SIZE: Number of genes in the gene set after filtering out those genes not in the expression datase.
• ES: Enrichment score for the gene set; that is, the degree to which this gene set is overrepresented at the top or bottom of the ranked list of genes in the expression dataset.
• NES: Normalized enrichment score; that is, the enrichment score for the gene set after it has been normalized across analyzed gene sets.
• NOM.p.val: Nominal p value; that is, the statistical significance of the enrichment score. The nominal p value is not adjusted for gene set size or multiple hypothesis testing; therefore, it is of limited use in comparing gene sets.
• FDR.q.val: False discovery rate; that is, the estimated probability that the normalized enrichment score represents a false positive finding.
• GS: short Gene set name
• no.genes.in.Core: number of genes contributes to the leading-edge subset within the gene set and contributes most to the enrichment result.

More details please see [http://www.gsea-msigdb.org/gsea/doc/GSEAUserGuideFrame.html?_Interpreting_GSEA_Results]

Please create the bubble plot below with the following criteria

• use 50 hallmark gene sets as y-axis and NES as x-axis.

  • if the NES > 0, order the Gene Sets based on NES values largest to smallest

  • if the NES < 0, order the Gene Sets based on the NES values smallest to largest

• FDR < 0.25 as threshold to show significant gene set with solid circle

  • hint: use scale_shape_manual(values=c(1,19)). 1: hollow circle; 19: solid circle

• use continuous colours to show FDR.q.val

  • hint: use scale_colour_gradient(name = "FDR.q.val",low = "#D55E00", high = "#0072B2")

• use the circle size to show no.genes.in.Core

Then save it as a pdf with width=8 and height=10