Data_in_R

Summary table

First extract the Summary numbers and create a summary table

• First extract the right rows, these begin with (^) SN.
• Then turn it into a table using the function separate (View the help of separate)
  • with 3 columns (ID, Name, and Value)
  • separate by the tab character “\t”
  • and remove the first column ‘[,-1]’, the SN
• Print the table using kable from the knitr package (makes a pretty looking table)

?separate

sn <- grep("^SN",data, value=TRUE)
sn <- separate(data.frame(sn),col=1, into=c("ID", "Name","Value"), sep="\t")[,-1]
kable(sn, caption="Summary numbers")

?kable

Exercise 1: While the Value column is numeric, by default it is being read in as characters. Use kable align parameter to left justify name and right justify value.

Get the next section, read lengths

First extract the read length data and create a table

• First extract the right rows, these begin (^) with RL.
• Then turn it into a table using the function separate (View the help of separate)
  • with 6 columns (ID, insert size, pairs total, inward oriented pairs, outward oriented pairs, other pairs)
  • separate by the tab character “\t”
  • and remove the first column ‘[,-1]’, the RL

rl <- grep("^RL",data, value=TRUE)
rl <- separate(data.frame(rl),col=1, into=c("ID", "read_length", "count"), sep="\t", convert = TRUE)[,-1]

Get the insert size of mapped pairs

First extract the insert sizes data and create a table

• First extract the right rows, these begin (^) with IS.
• Then turn it into a table using the function separate (View the help of separate)
  • with 6 columns (ID, insert size, pairs total, inward oriented pairs, outward oriented pairs, other pairs)
  • separate by the tab character “\t”
  • and remove the first column ‘[,-1]’, the IS

is <- grep("^IS",data, value=TRUE)
is <- separate(data.frame(is),col=1, into=c("ID", "insert size","all pairs", "inward", "outward", "other"), sep="\t", convert=TRUE)[,-1]

Get the ACGT content per cycle

First extract the base composition of first and last pairs and create a table

• First extract the right rows, these begin (^) with GCC.
• Then turn it into a table using the function separate (View the help of separate)
  • with 6 columns (ID, insert size, pairs total, inward oriented pairs, outward oriented pairs, other pairs)
  • separate by the tab character “\t”
  • and remove the first column ‘[,-1]’, the GCC

actg <- grep("^GCC",data, value=TRUE)
actg <- separate(data.frame(actg),col=1, into=c("ID", "cycle", "A", "C", "G", "T", "N", "O"), sep="\t",  convert=TRUE)[,-1]

Get the fragment qualities of mapped pairs

First extract the fragment qualities of first and last pairs and create a table

• First extract the right rows, these begin (^) with FFQ or LFQ.
• Then turn it into a table using the function separate (View the help of separate)
  • with 3 columns (Pair, Cycle, 1 … 43)
  • separate by the tab character “\t”

fq <- grep("^FFQ|^LFQ",data, value=TRUE)
fq <- separate(data.frame(fq),col=1, into=c("Pair", "Cycle", seq(41)), sep="\t", convert=TRUE)

## Warning: Expected 43 pieces. Missing pieces filled with `NA` in 279
## rows [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
## 20, ...].

We get a message here, saying data is missing. This is because there are no 38,39,40,41 quality scores (the typical range for Illumina qualities).

Get the GC content of mapped pairs

First extract the GC content of first and last pairs and create a table

• First extract the right rows, these begin (^) with GCF or GCL.
• Then turn it into a table using the function separate (View the help of separate)
  • with 3 columns (Pair, GC, Count)
  • separate by the tab character “\t”

gc <- grep("^GCF|^GCL",data, value=TRUE)
gc <- separate(data.frame(gc),col=1, into=c("Pair", "GC", "Count"), sep="\t", convert=TRUE)

Get the Indel Distribution

First extract the indel distribution data and create a table

• First extract the right rows, these begin (^) with ID.
• Then turn it into a table using the function separate (View the help of separate)
  • with 4 columns (ID, length, insertion_count, deletion_count)
  • separate by the tab character “\t”
  • and remove the first column ‘[,-1]’, the ID

id <- grep("^ID",data, value=TRUE)
id <- separate(data.frame(id),col=1, into=c("ID", "length", "insertion_count", "deletion_count"), sep="\t", covert=TRUE)[,-1]

Get the Indel per cycle

First extract the indel by cycle data and create a table

• First extract the right rows, these begin (^) with IC.
• Then turn it into a table using the function separate (View the help of separate)
  • with 6 columns (ID, cycle, ins_fwd, ins_rev, del_fwd, del_rev)
  • separate by the tab character “\t”
  • and remove the first column ‘[,-1]’, the IC

ic <- grep("^IC",data, value=TRUE)
ic <- separate(data.frame(ic),col=1, into=c("ID", "cycle", "ins_fwd", "ins_rev", "del_fwd", "del_rev"), sep="\t", convert=TRUE)[,-1]

Get the coverage data and GC Coverage data

Exercise 2: Use what you learned above to extract these 2 sections from the file.

Coverage data

• First extract the right rows, these begin (^) with COV.
• Then turn it into a table using the function separate (View the help of separate)
  • with 4 columns (ID, coverage_range, coverage, bases)
  • separate by the tab character “\t”
  • and remove the first column ‘[,-1]’, the COV

GC Coverage data

• First extract the right rows, these begin (^) with GCD.
• Then turn it into a table using the function separate (View the help of separate)
  • with 8 columns (ID, GC, GC_percentile, gc.10, gc.25, gc.50, gc.75, gc.90)
  • separate by the tab character “\t”
  • and remove the first column ‘[,-1]’, the GCD

Some Summary stats with

summarize(is,low=min(`insert size`), max=max(`insert size`), average=mean(`all pairs`), noutward=sum(outward), ninward=sum(inward))

##   low max average noutward   ninward
## 1   0 576  272522 32273430 120770240

new_is <- mutate(is,poutward=outward/`all pairs`, pinward=inward/`all pairs`)

Exercise 3: Try using “distinct”, on is (or new_is) on the outward and inward columns.

Exercise 4:

1. View the head/tail of some (or even all) of the objects.

2. Use dim to get an idea of the table dimentions.

3. Use summary to summarize and produce summary statistics (min, max, means, 1st and 3rd quartile boundaries) of the columns.

4. Any other summaries?

So now we have new objects (data.frames) that hold the data that we are interested in plotting

• Summary Numbers -> SN -> sn
• First Fragment Qualitites -> FFQ -> fq
• Last Fragment Qualitites -> LFQ -> fq
• GC Content of first fragments -> GCF -> gc
• GC Content of last fragments -> GCL -> gc
• ACGT content per cycle -> GCC -> actg
• Insert sizes -> IS -> is
• Read lengths -> RL -> rl
• Indel distribution -> ID -> id
• Indels per cycle -> IC -> ic
• Coverage distribution -> COV -> cov
• Coverage distribution -> GCD -> gccov

Our first plot, plotting the insert size of mapped fragments

We use the ggplot function and define the data as ‘is’ and x, y as get(“insert size”), get(“all pairs”), respectively. We use “get” because they have spaces in the names. There is another way to do it. Hint: look at the documentation above for clues.

g <- ggplot(data = is)
g + geom_line( aes(x=get("insert size"), y=get("all pairs")))

Ok, now add some labels to the plot

g + geom_line( aes(x=get("insert size"), y=get("all pairs"))) + 
  labs( x = "insert size", y = "all pairs", title ="Mapped insert sizes", subtitle = "All Pairs", caption = "all pairs insert size")

Ok, what about plotting multiple data objects on the same plot (multiple lines)? We can specifically set the y axis in geom_line and color, then call geom_lines a second time (or more times).

g <- ggplot(data = is, aes(x=get("insert size")))
g + geom_line(aes(y=get("inward")),color="blue") +  
    geom_line(aes(y=get("outward")),color="orange") + 
    labs( x = "insert size", y = "all pairs", title ="Mapped insert sizes", subtitle = "All Pairs", caption = "all pairs insert size")

try adjusting the x/y limits to 0,600 and 0,20000 respectively.

g + geom_line(aes(y=get("inward")),color="blue") +
  geom_line(aes(y=get("outward")),color="orange") + 
  coord_cartesian(xlim=c(0,500), ylim=c(0,600000))

Ok so now put all these elements together into a single plot, save final plot as ‘g’

g <- ggplot(data = is, aes(x=get("insert size")))
g <- g + geom_line(aes(y=get("all pairs")), color="black") +  
    geom_line(aes(y=get("inward")),color="blue") +  
    geom_line(aes(y=get("outward")),color="orange") +
    geom_line(aes(y=get("other")), color="green")
g <- g + 
    labs( x = "insert size", y = "all pairs", title ="Mapped insert sizes", subtitle = "All Pairs", caption = "all pairs insert size")
g <- g + coord_cartesian(xlim=c(0,500), ylim=c(0,600000))
g <- g + theme_light()
plot(g)

Exercise 5: Put it all together, plot all four columns of the insert size data object, add in legends, reasonable coordinate limits.

Exercise 6: Play with ggplot2 themes (ex. theme_classic() )

Plotting GC content

In order to plot GC percentage we first need to convert the counts to proportions, to do so we can divide the counts by the sum of counts.

head(gc)

##   Pair   GC Count
## 1  GCF 0.25  9986
## 2  GCF 1.01  6442
## 3  GCF 1.76  6816
## 4  GCF 2.51  8029
## 5  GCF 3.27  9586
## 6  GCF 4.02 11557

h <- ggplot(gc, aes(GC, Count/sum(Count),color=Pair))
h <- h + geom_line()
h

Exercise 7: Finish the plot (add labels, etc.). Save the final graph object in h.

Plotting the base composition by cycle

Sometimes we may need to transform our data before plotting. The melt funciton from reshape2 takes data in wide format (data are in columns) and stacks a set of columns into a single column of data. In the ACTG object we can stack bases values by cycle.

actgm <- melt(actg,id="cycle")

now head the new actgm object. What did melt do?

ic <- ggplot(actgm, aes(cycle, value, by=variable, colour=variable))
i <- ic + geom_line() + coord_cartesian(ylim=c(0,100))
i

Exercise 8: Using what you have learned so far, finish the plot, save it as object i.

A boxplot of basepair

i2 <- ic + geom_boxplot()
i2

Exercise 9: Try some other geometries (Ex. bin2d, col, count, which generate an ‘interpretable’ plot).

Plotting a heatmap of qualities

First melt the quality scores

fqm <- melt(fq,id=c("Pair","Cycle"))

Take a look at the new object

j <- ggplot(fqm, aes(Cycle, variable)) 
j + geom_tile(aes(fill = as.numeric(value)))

Now try changing the gradient colors and modify the legend, add labels. The ggplot2 ‘theme’ function can be used to modify individual components of a theme.

?theme

j = j + geom_tile(aes(fill = as.numeric(value))) + 
  scale_fill_gradient(low = "red", high = "green") +
  ylab("Cycle") +
  xlab("Quality") +
  theme(legend.title = element_text(size = 10),
        legend.text = element_text(size = 12),
        plot.title = element_text(size=16),
        axis.title=element_text(size=14,face="bold"),
        axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(fill = "Quality value")
j

Exercise 10: Try modifying scale_fill_gradient to scale_fill_distiller.

Exercise 11: Play with parts of the plotting function, see how the change modifies the plot.

Plotting indel lengths

k <- ggplot(id, aes(x=as.numeric(length)))
k <- k + geom_line(aes(y=as.numeric(insertion_count)), color = "red", size=1.5)
k <- k + geom_line(aes(y=as.numeric(deletion_count)), color = "black", size=1.5)
k

Try changing the Y axis to log scale

k <- k + scale_y_log10()
k

Tweak the grid elments using theme

k <- k + theme(panel.grid.minor = element_blank(), 
  panel.grid.major = element_line(color = "gray50", size = 0.5),
  panel.grid.major.x = element_blank())
k

## Warning: Transformation introduced infinite values in continuous y-axis

Fancy Multiple plots in a grid

The gridExtra package is great for plotting multiple object in one plot.

include_graphics("https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2018-September-Bioinformatics-Prerequisites/master/thursday/Data_in_R/grid_plot.png")

full <- grid.arrange(
  g, h, i, i2, 
  widths = c(2, 1, 1),
  layout_matrix = rbind(c(1, 2, NA),
                        c(3, 3, 4))
)

Exercise 12: Play with th grid.arrange function, using the plots you’ve created to create you own final combined plot.

Saving plots as png or pdf

This might have to be done outside of the Notebook as the notebook may expect you to plot in the notebook only, so run on the Console if you have trouble running in the Notebook.

Saving plots to pdf ** do on the console if having trouble in Notebook **

ggsave("multi_plot.pdf",full,device="pdf",width=6,height=4, units="in", dpi=300)

Saving plots to png ** do on the console if having trouble in Notebook **

ggsave("multi_plot.png",full,device="png",width=6,height=4, units="in", dpi=300)

View the help documentation for ggsave, what other

With any remaining time (or homework), use the ggplot cheat sheet to further expand and modify the plots.