R_programming2016

Introduction

R-studio helps greatly in doing reproducable research via the R markdown language. In a Rmd document you can r code, r results and plain text in one document that then can be exported in many different formats.

Code blocks can be created both inline with text 4 + 2 = 6 or in blocks of code:

4 + 2

[1] 6

RNA seq NBIS report

In this document I will create two plots that one often creates in analysing RNA sequence data.

The data is from a recent paper on YAP and TAZ control of peripheral myelination and the expression of laminin receptors in Schwann cells (Poitelon et al. Nature Neurosci. 2016). In the study the genes YAP and TAZ were knocked-down in Schwann cells.

The material for RNA-seq was collected from three biological replicates of both knockdowns and controls.

Both read mapping and read counting has been done outside R and the code below reads in the read counts per gene and create an object for analysis of differential gene expression.

library(edgeR)
## Preparing samples group labels
sample.groups <- rep(c('KO','Wt'), each = 3)
## Reading in data
data.counts <- read.table("tableCounts.txt", header = TRUE)
#colnames(data.counts) <- sample.groups 
## edgeR
## building edgeR object
data.cds <- DGEList(data.counts, group = sample.groups)
## filtering low count reads (keeping rads with at least 1 read per milion in at least 2 samples)
data.cds <- data.cds[rowSums(cpm(data.cds)> 1) > 2,]
## normalizing data
data.cds <- calcNormFactors(data.cds)

MDS plot

Multidimensional scaling plots can be used to visualize distances between samples in a gene expression experiment. Note that there are a seperation of sample along the major axis of differentiation indicating that there like are genes showing differential gene expression between knockdown and wildtype.

plotMDS(data.cds, labels = sample.groups)

Detection of differentially expressed genes

my.design <- model.matrix(~0+sample.groups)
colnames(my.design) <- levels(data.cds$samples$group)
data.cds <- estimateDisp(data.cds, my.design, trend.method='movingave')
my.contrasts <- makeContrasts(KOvsWt = KO-Wt, levels = my.design)
fit <- glmFit(data.cds, my.design)
lrt <- glmLRT(fit, contrast=my.contrasts[,"KOvsWt"])
lrt.top <- topTags(lrt, n=Inf, adjust.method='BH')
lrt.table <- lrt.top$table
de <- decideTestsDGE(lrt, lfc=0)
detags <- rownames(data.cds)[as.logical(de)]
summary(de)

   [,1] 
-1   920
0  11939
1    672

MA-plot

plotSmear(lrt, de.tags = detags)

R session

A useful command for reproducable research in R is the sessioninfo. It prints a information on all loaded packages in your the current R session and will hence be helpful if you share the code or want to rerun the code at a later time and have forgotten the packages to load.

sessionInfo()

R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] edgeR_3.14.0  limma_3.28.18

loaded via a namespace (and not attached):
[1] formatR_1.4   tools_3.3.1   splines_3.3.1 knitr_1.14