RNA seq: Differential Gene Expression analysis - AnimaTardeb/Meiogenix-UGDR GitHub Wiki
This wiki is just notes of my work. However, it could be useful for some beginners in differential gene analysis (DGE). Several tutorials are present to help in DGE analysis. Following those tutorials, I was able to create my own pipeline to analyse RNASeq data. The most relevant for me are:
Many tools/methods are available to analyse differential gene expression and each is builds on specific model. Various studies described comparisons between those tools/methods see this article.
In this wiki, I will be using DESeq2 to analyse my RNASeq data for DGE.
- I like loading all the packages, inputs and output files (specify the path) and folders before starting my analysis
source("http://www.Bioconductor.org/biocLite.R")
#biocLite( c("ShortRead","DESeq", "edgeR") )
library( "DESeq2" )
library("RColorBrewer")
library("ggplot2") #Best plots
library("ggrepel") #Avoid overlapping labels
library("ggpubr")
#Read input file
inputcondition = "smooth_muscle_contractility3.csv"
inputMcount = "MergedCounts.csv"
inputlengthgene= "genelengthMeng.txt"
gf = "Homo_sapiens.GRCh37.75.gtf"
Length=read.table("R-gene_name_length.txt", header=T,sep="\t",row.names = 1)
outputcpm = "2019_ABgorupsCPM2.txt"
outputrpkm = "2019_ABgorupsrpkm2.txt"
#read metadata
sri = read.csv(inputcondition, stringsAsFactors=FALSE, sep=";")
sri = sri[order(sri$donors),]
#read counts matrix
MergedCount=read.table(inputMcount, header=T,sep=";",row.names = 1)
gene=read.table(inputlengthgene, header=T,sep="\t",row.names = 1)- Create DESeq2Dataset object from matrix:
dds = DESeqDataSetFromMatrix(countData = MergedCount , colData = sri , design =~groups+donors))
Important: informations on how to design the formula are here
- Normalize the counts, estimate dispersion and fit curve to gene-wide dispersion. All these steps are carried out by Deseq() function. for more information about these steps See .
NOTE 1: DESeq2 will model the raw counts, using normalization factors (size factors) to account for differences in library depth. Then, it will estimate the gene-wise dispersions and shrink these estimates to generate more accurate estimates of dispersion to model the counts. Finally, DESeq2 will fit the negative binomial model and perform hypothesis testing using the Wild test or Likelihood Ratio Test. Note Credit: hbctraining.github.io.
However, command line to :
dds = estimateSizeFactors(dds)#send back the results to dds (DESeqDataSet object).
dds = normalizationFactors(dds)
sizeFactors(dds)
plotDispEsts(dds)
NOTE 2: DESeq2 doesn’t actually use normalized counts, rather it uses the raw counts and models the normalization inside the Generalized Linear Model (GLM). These normalized counts will be useful for downstream visualization of results, but cannot be used as input to DESeq2 or any other tools that peform differential expression analysis which use the negative binomial model. Note Credit: hbctraining.github.io.
in other words we can't use deseq2(normalized_counts).
- PCA Analysis
The PCA plot is useful to see the similarity between samples, the overall effect of the experiments and batch effect. Moreover, PCA plots help in identifying samples outliers.
Deseq2 used regulated log transform (rlog)
library(ggplot2)
rld <- rlogTransformation(dds, blind=TRUE)
# vst() would run faster for more than 20 samples
#blind=TRUE => transformation unbiased by prior information on samples .
png("PCA.png", width = 7*300, height = 5*300, res = 400, pointsize = 12)
plotPCA(rld, intgroup=c('groups', 'samples"),ntop = 6000, returnData = F)
p <- plotPCA(rld,intgroup=c('groups'),ntop = 6000, returnData = F)+theme_light()
p + geom_text(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black") + theme_light()
dev.off()- RUN Deseq2()
know, You can run the Deseq2 for differential gene expression analysis.
dds <- DESeq(dds)
- Results
-
Plot the dispersion
plotDispEsts(dds)that should be identical to the one plotted earlier. -
Extract the results.
contrast=c("groups","B","A")
res <- results(dds, contrast=contrast)
#inspecting before showing the results.
mcols(res, use.names=TRUE)
summary(res)
#select the genes with FDR>0.05 and FDR<-0.05
up=sum(res$padj < 0.05 & res$log2FoldChange >= 0.5 , na.rm=TRUE)
down=sum(res$padj < 0.05 & res$log2FoldChange <= -0.5, na.rm=TRUE)
up
down- It would be interesting to see the distribution of the results P-value.
#plot the P_value
png("mRNA_pval.png", width = 7*300, height = 5*300, res = 400, pointsize = 12)
hist(res$pvalue,breaks=100,col="grey", xlab="p-value",main="p-value distribution")
dev.off() #can do better plot- Retrieve the results into files
#Write tables:
respadjfc <- res[order(res$padj, decreasing = FALSE), ]
CountsN<-counts(dds, normalized=TRUE)# or data2=fpm(dds, normalized=TRUE)
write.table(respadjfc,file="mRNA-results.txt",sep="\t",quote=F)
write.table(CountsN, file="mRNA-normalizedCounts.txt", sep="\t")
- At this step, I will have two files. The first file
mRNA-results.txtcontains: baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, the second one contains normalized counts (it can be FPM or FPKM depends on what you did earlier). The good thing is to merge these tow files into one file (See figure 1).
- To Merge the two first files see bellow:
#first define the output file (you can write all the files for each step, or just the last one
#need some information from the gene length file.
length=data.frame(rownames(Length),Length$name, Length$length)
#Merge table of results to gene name table:
idx=match(length$rownames.Length,rownames(respadjfc))
Lengthmatch = as.data.frame(length[!is.na(idx),])
tDESEQmatch = respadjfc[idx[!is.na(idx)],]
Table1=cbind(tDESEQmatch, Lengthmatch)
#write.table(Table1[order(Table1$padj),],file=outputfile2,sep=";",quote=F)
#Merge table of counts to gene name table
idx=match(length$rownames.Length,rownames(data2))
Lengthmatch = as.data.frame(length[!is.na(idx),])
tDESEQmatch = data2[idx[!is.na(idx)],]
Table2=cbind(tDESEQmatch, Lengthmatch)
#write.table(Table2[order(Table2$padj),],file=outputfile1,sep=";",quote=F)
#Merge Table1 and Table2 by matching the genenames
idx=match(rownames(Table1),rownames(Table2))
Lengthmatch = as.data.frame(Table1[!is.na(idx),])
tDESEQmatch = Table2[idx[!is.na(idx)],]
Table3=cbind(tDESEQmatch, Lengthmatch)
#The last file that merges the normalized counts and the results matched by gene name.
write.table(Table3[order(Table3$padj),],file=outputfile3,sep="\t",quote=F)- Plot Significant Results: Volcano plot and/or MAplot.
#Use the outputfile2 (table1).
name=read.table("outputfile2", header=T,sep=";",row.names = 1)
#another type of maplot
png("mRNA_MAggplot.png", width = 7*300, height = 5*300, res = 400, pointsize = 12)
ggmaplot(name, main = expression("B. vs A."),
fdr = 0.05, fc = 0.75, size = 1,
palette = c("#B31B21", "#1465AC", "darkgray"),
genenames = as.vector(name[,"Length.name"]),#for mRNA
legend = "top", top = 30,
font.label = c("bold", 5),
font.legend = "bold",
font.main = "bold",
ggtheme = ggplot2::theme_light())
dev.off()
- Or Volcano plot
#Script modified from the original one. However I don't remember where I get the original script. SORRY.
png("mRNA_volcanoggplot.png", width = 7*400, height = 5*300, res = 400, pointsize = 12)
mydata<-name%>%mutate(threshold = ifelse(log2FoldChange >= 1 & padj<0.05 ,"UP", ifelse(log2FoldChange<=-1 & padj<0.05, "Down", "NA")))
mydata2=subset(mydata, mydata$threshold=="UP" | mydata$threshold=="Down")
ggplot(mydata, aes(x=log2FoldChange, y=-log10(padj)),fdr = 0.05, fc = 1, size = 100,) + xlim(-2, 2)+#ylim(0,10)+
geom_point(aes(colour = threshold), size=1.5) +theme_light()+
theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)), text = element_text(size=14))+
scale_colour_manual(values = c("UP"= "#B31B21", "Down"="#1465AC", "NA"= "darkgray"))+
geom_text_repel(data=head(mydata2, 40), aes(label=Length.name) ,size = 2)+ #adding text for the top 20 genes
geom_vline(xintercept=0, linetype="dashed",color = "black", size=.5)
dev.off()- TO BE CONTINUED ...