Skip to content

Latest commit

 

History

History
427 lines (270 loc) · 10.3 KB

pathway_analysis.md

File metadata and controls

427 lines (270 loc) · 10.3 KB

Pathway Analysis

based on a tutorial by Asela Wijeratne

###During this session you will learn about:

  1. Use of GAGE package to do pathway analysis

    a. "Native workflow"

    b. Combined workflow with RNAseq

  2. Interpret the GAGE output

  3. Use of Pathview to visualize the perturbed KEGG pathways

######For more details, please go to relevent reference manauls: GAGE package and tutorial:

http://bioconductor.org/packages/devel/bioc/html/gage.html

pathview and gageData:

http://bioconductor.org/packages/devel/bioc/html/pathview.html

http://bioconductor.org/packages/devel/data/experiment/html/gageData.html EdgeR user guide:

http://bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf

###Installing necessary packages

Before you start the tutorial following packages needs to be installed. Install required packages:

source("http://bioconductor.org/biocLite.R")
biocLite(c("pathview", "edgeR", "gage"))

Load required packages

library(edgeR)
library(gage)
library(pathview)

###Getting and preprocessing data

First, let's get some RNAseq data. This data comes from a study described in Zhang et al., 2011, Genome-wide mapping of the HY5-mediated gene-networks in Arabidopsis that involve both transcriptional and post-transcriptional regulation (http://onlinelibrary.wiley.com/doi/10.1111/j.1365-313X.2010.04426.x/full). This a comparison between a Arabidopsis mutant called hy5 and wild-type. For each genotype, there are two replicates. Vince Buffalo has re-analyzed this data as described here: (https://github.com/vsbuffalo/rna-seq-example). We are going to grab the count data file from this analysis.

rnaseq_URL <- "https://github.com/vsbuffalo/rna-seq-example/raw/master/results/raw-counts.txt"
download.file(rnaseq_URL, destfile = "rnaseq_data.txt")
rnaseq_data <- read.table("rnaseq_data.txt",
                          header=TRUE)

#####If above code work, ignore the rest of the next part of the code.

You can also get the data from the following link and save it in the Desktop/pathway_analysis folder.

http://bit.ly/150826_rnaseqdata

Read the data into R using the following code:

rnaseq_data <- read.table("rnaseq_data.txt",
                          header=TRUE)

Let's look at the head of the data set.


head(rnaseq_data)

This is how they are labelled:

GEO:GSM613465: wild type samples
        SRR070570: replicate 1
        SRR070571: replicate 2
GEO:GSM613466: hy5 mutant samples
        SRR070572: replicate 1
        SRR070573: replicate 2

Rename the columns so we can identify them easily.

colnames(rnaseq_data) <- c("WT_1", "WT_2", "hy_1", "hy_2")

###Getting KEGG data

#####1. Finding your species in KEGG First, we need to make sure the species we work with is found in KEGG database. You can go the following webpage and see the organism is list:

http://rest.kegg.jp/list/organism

In addition, pathview also carry a data matrix with KEGG supported species. You can explore this matrix as well.

data(korg)
head(korg)

You can get the short names for Arabidopsis by searching the korg data.

org <- "arabidopsis thaliana"
species <- unlist(sapply(1:ncol(korg), function(i) {
    agrep(org, korg[, i])
}))
korg[species, 1, drop = F]

#####2. Creating the KEGG dataset for GAGE analysis kegg.gsets can be used to get KEGG data for any species present in the KEGG database. GAGE manual recommends that you save this data as a .Rdata file. This way you don't need to download this each time you need to use and also increase the reproducibility.

kegg_arab <- kegg.gsets("ath")
kegg.gs <- kegg_arab$kg.sets[kegg_arab$sigmet.idx]

###GAGE "native" workflow

#####1. Remove zero counts in all samples.

rnaseq_counts <- rnaseq_data # rename the dataset 

dim(rnaseq_counts)
non_zero <- rowSums(rnaseq_counts) != 0
rnaseq_counts <- rnaseq_counts[non_zero,]
# number of genes left after removing zero counts 

dim(rnaseq_counts)

#####2. Normalize library sizes

libsizes <- colSums(rnaseq_counts)
size_factor <- libsizes / exp(mean(log(libsizes)))
norm_counts <- t(t(rnaseq_counts) / size_factor)
range(norm_counts)

#####3. Variance stabilizing transformation We need to add small (positive) value to avoid getting -inf during log transformation. We will add 8 as per GAGE manual.

norm_counts <- log2(norm_counts + 8)
range(norm_counts)

You can do MA and PCA plots to see the processed data variances and overall variances and similarity between samples, respectively.

#####4. Define reference and experiment samples

ref_idx <- 1:2 #wt
samp_idx <- 3:4 #hy mutant

#####5. Enrichment analysis Here we are going to perform the enrichment analysis using fold change (default) as the gene set statistics. In addition, control (wt) and experiment (hy) are not paired. Therefore, we going to use compare ="unpaired".

native_kegg_fc <- gage(norm_counts, 
                       gsets = kegg.gs, 
                       ref = ref_idx, 
                       samp = samp_idx, 
                       compare ="unpaired")

#####6. Some outputs

First four pathways that are up-regulated.

head(native_kegg_fc$greater[,1:5], 4)

######Meaning of the output:

p.geomean geometric mean of the individual p-values from multiple single array based gene set tests

stat.mean mean of the individual statistics from multiple single array based gene set tests

p.val gloal p-value or summary of the individual p-values from multiple single array based gene set tests. This is the default p-value being used.

q.val FDR q-value adjustment of the global p-value using the Benjamini & Hochberg procedure implemented in multtest package. This is the default q-value being used.

set.size the effective gene set size, i.e. the number of genes included in the gene set test

First four pathways that are down-regulated.

head(native_kegg_fc$less[,1:5], 4)

We can also do significance test to find the pathways that are up-regulated and down-regulated using a given cutoff (default is 0.1). sigGeneSet function will also create a heat map, showing the average fold change of expression in pathways that are perturbed.

wt_hy_sig_kegg <-sigGeneSet(native_kegg_fc, outname="wt_hy.kegg")

Write the significant pathway sets to a text file.

write.table(rbind(wt_hy_sig_kegg$greater,
                  wt_hy_sig_kegg$less), 
            file = "wt_hy_sig_kegg.txt", sep = "\t")

###Combine GAGE and EdgeR analysis

First, we will do a differential gene expression analysis using the edgeR package.

#####1. Make the study design as follows:

targets <-
        data.frame(sample_name=colnames(rnaseq_data),Group=rep(c("WT","hy"),each=2))

targets

#####2. Create DGElist object.

This is a R list object that can be easily manipulated.



d <- DGEList(counts=rnaseq_data, 
             group=targets$Group) # Constructs DGEList object

#####3. Remove lowly expressed genes Keep genes with more than 10 counts per million (calculated with cpm()) at least in two samples.

#before filtering:
dim(d)
keep <- rowSums(cpm(d)> 10) >= 2
d <- d[keep,]
#after filtering:
dim(d)

#####4. Normalize the data

d$samples$lib.size <- colSums(d$counts)
d <- calcNormFactors(d)
d

#####5. multidimentioanl scaling plot Quick sanity check to see how samples are related to each other using multidimensional scaling plot (MSD).


plotMDS(d, labels = targets$Group, 
        col = c("darkgreen","blue")[factor(targets$Group)])

#####6. Estimate the dispersion

d1 <- estimateCommonDisp(d, verbose=T)

d1 <- estimateTagwiseDisp(d1)

plotBCV(d1)

#####7. Differentail gene expression

de.com <- exactTest(d1, pair = c("WT", "hy"))
topTags(de.com, n = 5)


FDR <- p.adjust(de.com$table$PValue, method="BH")
de1 <- decideTestsDGE(de.com, adjust.method="BH", p.value=0.05)
summary(de1)

#####8. Making the data suitable for pathway analysis

isDE <- as.logical(de1) # covert to DE set to true/false set

DEnames <- rownames(d)[isDE] # get the DE gene names
edger.fc <- de.com$table$logFC # get the log fold change
names(edger.fc) <- rownames(de.com$table) # assign row names to fold change 
exp.fc <- edger.fc
head(exp.fc)

Getting DE gene set and fold change

DE_foldchange <- edger.fc[DEnames]

length(DE_foldchange)

#####9. GAGE analysis

fc.kegg.p <- gage(exp.fc, gsets = kegg.gs, ref = NULL, samp = NULL)


head(fc.kegg.p$greater[,1:5], 4)

Significance test

wt_hy_sig_kegg <-sigGeneSet(fc.kegg.p, outname="wt_hy.kegg")

###Visualize the data using Pathview

log_fc= norm_counts[, samp_idx]-rowMeans(norm_counts[, ref_idx])

#####1. Selecting the path ids for the upregulated set.

greater_set <- native_kegg_fc$greater[, "q.val"] < 0.1 &
        !is.na(native_kegg_fc$greater[, "q.val"])
greater_ids <- rownames(native_kegg_fc$greater)[greater_set]
head(greater_ids)

#####2. Selecting path ids for down-regulated set

less_set <- native_kegg_fc$less[, "q.val"] < 0.1 &
        !is.na(native_kegg_fc$less[,"q.val"])
less_ids <- rownames(native_kegg_fc$less)[less_set]

#####3. Combine up and down-regulated path ids.

combine_ids <- substr(c(greater_ids, less_ids), 1, 8)
head(combine_ids)

#####4. Visualization

Here we are going to get the first three pathways and visualize using pathview.

pv.out.list <- sapply(combine_ids[1:3], function(pid) pathview(
        gene.data =  exp.fc, pathway.id = pid,
        species = "ath", out.suffix="Wt_hy",
        gene.idtype="KEGG"))

#####5. Mapping DE genes to individual pathways

We can map genes that show statistically significant changes to individual pathways. As an example, we can take DNA replication pathway and see which of the DE genes from EdgeR is present.

pv_replication <- pathview(gene.data = DE_foldchange, 
                           gene.idtype = "KEGG", 
                           pathway.id = combine_ids[3], 
                           species = "ath", 
                           out.suffix = "DNA_replication", 
                           keys.align = "y", 
                           kegg.native = T, 
                           match.data = T, 
                           key.pos = "topright")

See how the output look like.

head(pv_replication)