Standard edgeR Workflow: edgeR Video Tutorial 1

[qhawenid]

Great stuff as always!! Please also do pseudobulk RNA with DESeq2. Thanks much!

[Collection_of_online_tutorials]

### Standard edgeR Workflow for Bulk/Pseudo RNA-seq Data Analysis
library(tidyverse)
library(limma)
library(edgeR)

### Step 1: load RNA-seq data (oral carcinomas vs matched normal tissue)
Oral_RNA <- read.csv("../edgeR/Oral_RNA.csv", header = T, sep = ", ")


### Step 2: create a DGEList
Tissue <- factor(c("N", "T", "N", "T", "N", "T"))
y <- DGEList(counts = Oral_RNA[, 3:8], genes= Oral_RNA[, 1:2], group = Tissue)


### Step 3: Filtering
keep <- filterByExpr(y)
table(keep)
y <- y[keep, , keep.lib.sizes = FALSE]


### Step 4: Normalization
y <- normLibSizes(y)


### Step 5: Data exploration
plotMDS(y)


### Step 6: Design matrix (according to your research question)
Patient <- factor(c(8, 8, 33, 33, 51, 51))
design <- model.matrix(~Patient + Tissue)
rownames(design) <- colnames(y)


### Step 7: Dispersion estimation
y <- estimateDisp(y, design, robust = TRUE)
plotBCV(y)


### Step 8: Differential expression
# To perform quasi-likelihood F-tests:
fit <- glmQLFit(y, design) #Genewise negative binomial generalized linear models
qlf <- glmQLFTest(fit) # perform quasi-likelihood F-tests (qlf)
topTags(qlf, n=20)

A <- qlf[["table"]]
A <- filter(A, PValue<0.05)
write.csv(A, "../edgeR/A.csv")

# To perform likelihood ratio tests:
fit <- glmFit(y, design) # Genewise negative binomial generalized linear models
lrt <- glmLRT(fit) # perform likelihood ratio tests (lrt)
topTags(lrt, n=10)

B <- lrt[["table"]]
B <- filter(B, PValue<0.05)
write.csv(B, "../edgeR/B.csv")

# Common genes
A_genes <- row.names(A)
B_genes <- row.names(B)
c <- intersect(A_genes, B_genes)

# look at the counts-per-million in individual samples for the top genes
o <- order(lrt$table$PValue)
cpm(y)[o[1:10], ]
summary(decideTests(lrt))

# Plot log-fold change against log-counts per million
plotMD(lrt)
abline(h=c(-1, 1), col = "blue")


### Step 9: Pathway analysis
# Gene ontology analysis
go <- goana(lrt)
topGO(go, ont = "BP", sort = "Up", n = 30, truncate = 30)

# KEGG analysis
keg <- kegga(lrt)
topKEGG(keg, sort = "up")

[HopeOverDebt]

share link from where you have downloaded data please