--- title: "RORa effect on exercise - RNAseq DE Analysis" author: "wandrille Duchemin" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ## setting up environnement ```{r , message=F} baseDir="." library(DESeq2) library(RColorBrewer) library(ggplot2) library(gridExtra) library(pheatmap) ## recipe for a heatmap makeHeatMap <- function( VSTdata , sampleNames , main='' ) { sampleDists <- dist(t(assay(VSTdata))) sampleDistMatrix <- as.matrix(sampleDists) rownames(sampleDistMatrix) <- sampleNames colnames(sampleDistMatrix) <- NULL colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255) X = pheatmap(sampleDistMatrix, clustering_distance_rows=sampleDists, clustering_distance_cols=sampleDists, col=colors , main=main) return(X) } ``` ## setting up data ```{r } ## count matrix. # - rows are genes # - columns are sample countMatrixFile = paste(baseDir,'05_expression_matrix.txt',sep='/') cts_raw <- read.csv(countMatrixFile,sep="\t", header=T) # samples MUST be in the same order in count matrix and info table # (MUST?) have the same name in count matrix and info table cts = as.matrix(cts_raw) #head(cts) dim(cts) ``` ```{r} colnames(cts) ``` ```{r } ## information table # - rows are sample # - columns are variable library(stringr) library(plyr) genotype = as.factor( str_detect( colnames(cts) , 'CTRL' ) ) genotype = revalue( genotype , c( 'TRUE'="WT" , 'FALSE' = 'KO' ) ) time_point = factor( substring(colnames(cts), regexpr("_Z",colnames(cts)) + 2) , levels = c("T0",'T4','T8','T12','T16','T20')) mouse = sapply( strsplit( colnames( cts ) , "_") , '[' , 5) replicate = substring(colnames(cts) , 1, nchar( "HL3YYBGX5" ) ) coldata <- data.frame( sample = colnames(cts) , genotype = genotype , time_point = time_point, mouse = mouse, replicate = replicate ) head( coldata ) ``` ### ensembl to gene name table ```{r} library(clusterProfiler) library(org.Mm.eg.db) entrezId = bitr( rownames( cts ) , fromType="ENSEMBL", toType="SYMBOL", OrgDb="org.Mm.eg.db") ## problem -> n:n mapping from ensembl to entrez ## -> for now remove all duplicates on both sides nbefore=nrow(entrezId) entrezId = entrezId[!duplicated(entrezId$ENSEMBL),] entrezId = entrezId[!duplicated(entrezId$SYMBOL),] nafter=nrow(entrezId) print(paste("removed",nbefore-nafter,"multiplicates")) print(paste( ' genes with a n Entrez id', length( entrezId$SYMBOL ) , '/' , length( rownames( cts ) ) , '(' ,length( entrezId$SYMBOL ) / length( rownames( cts ) ) , ')' ) ) rownames( entrezId ) = entrezId$ENSEMBL ``` ## 1. pre-processing : looking at batch effect and potential outliers ### count data filtering Let's apply a filtering of at least 10 reads counted in at least a number of samples equivalent to the number of replicates : 3. ```{r} totalNbGenes = nrow(cts) NbGenesNonEmpty = sum(apply(cts, 1, sum)>0) threshold = 10 nbSample = 3 filterOK <- function( row , t = threshold , n = nbSample ){ sum(row>=t)>=n } genesPassing = apply(cts, 1, filterOK) nbGenesPassing = sum(genesPassing) cts_filtered = cts[genesPassing,] print( paste( 'total number of genes' , totalNbGenes )) print( paste( 'number of genes with at least 1 read' , NbGenesNonEmpty )) print( paste( 'number of genes passing the filter' , nbGenesPassing )) dim(cts_filtered) #head(cts_filtered) ``` ### DESeq ```{r} # design should be a formula using the column of the info table. # In order to benefit from the default settings of the package, # you should put the variable of interest at the end of the formula # and make sure the control level is the first level. RNAdata <- DESeqDataSetFromMatrix(countData = cts_filtered, colData = coldata, design = ~ mouse + replicate) ``` ### PCA ```{r} library(plotly) # variance stabilisation transform RNAdata_VST <- vst(RNAdata, blind=FALSE) ``` ```{r} pairedPCAplot = function (object, intgroup = "condition", ntop = 500, returnData = FALSE , pairing = '') { rv <- rowVars(assay(object)) select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))] pca <- prcomp(t(assay(object)[select, ])) percentVar <- pca$sdev^2/sum(pca$sdev^2) if (!all(intgroup %in% names(colData(object)))) { stop("the argument 'intgroup' should specify columns of colData(dds)") } intgroup.df <- as.data.frame(colData(object)[, intgroup, drop = FALSE]) group <- if (length(intgroup) > 1) { factor(apply(intgroup.df, 1, paste, collapse = ":")) } else { colData(object)[[intgroup]] } d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], group = group, intgroup.df, name = colnames(object) , pairing = colData(object)[[pairing]] ) if (returnData) { attr(d, "percentVar") <- percentVar[1:2] return(d) } ggplot(data = d, aes_string(x = "PC1", y = "PC2", color = "group" )) + geom_point(size = 3, aes(label=d$group)) + xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) + ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance")) + coord_fixed() + geom_line(aes(group=pairing)) } ``` ```{r} g = pairedPCAplot(RNAdata_VST, intgroup=c("replicate") , pairing = 'mouse' ) + ggtitle("PCA - filtered 10reads 3samples - VST") ggplotly(g) ``` Zooming in. we can see that mice group very well by replicate. However, The relationship between these does not look random. Let's investigate if a simple model can find DE genes associated to replicate: ```{r} RNAdata <- DESeq(RNAdata) ``` ```{r} FDRthreshold = 0.01 resultsNames(RNAdata) res <- results(RNAdata, name= "replicate_HL3YYBGX5_vs_HL3YTBGX5" ,alpha = FDRthreshold ) print(summary(res)) res <- results(RNAdata, name= "replicate_HL7CHBGX5_vs_HL3YTBGX5" ,alpha = FDRthreshold ) print(summary(res)) ``` Indeed we can find a significant signal associated with technical replicates. As each mouse is present once per technical replication, we are not prevented from answering any questions we may have, but we will not be collapsing technical replicates and have to take them into account as a co-variable in our model. ```{r} g = pairedPCAplot(RNAdata_VST, intgroup=c("mouse",'genotype',"time_point") , pairing = 'time_point' ) + ggtitle("PCA - filtered 10reads 3samples - VST") ggplotly(g) ``` Another issue to address is outliers : mice 26, 15, 23, 12. They do not particularly pop out in the QC report. We will be removing them ### removing 26,15,23,12 ```{r} toKeep = !( str_detect( colnames(cts) , '_26_') | str_detect( colnames(cts) , '_15_') | str_detect( colnames(cts) , '_23_') | str_detect( colnames(cts) , '_12_') ) cts_no_outlier = cts[ , toKeep ] print( dim( cts_no_outlier ) ) coldata_no_outlier = coldata[ !(coldata$mouse %in% c(26,15,23,12)) ,] print( dim( coldata_no_outlier ) ) ``` ## 2. replicates non collapsed, removed "outliers" ### 2.1. count data filtering Let's apply a filtering of at least 10 reads counted in at least a number of samples equivalent to the number of replicates : 3. ```{r} totalNbGenes = nrow(cts_no_outlier) NbGenesNonEmpty = sum(apply(cts_no_outlier, 1, sum)>0) threshold = 10 nbSample = 3 filterOK <- function( row , t = threshold , n = nbSample ){ sum(row>=t)>=n } genesPassing = apply(cts_no_outlier, 1, filterOK) nbGenesPassing = sum(genesPassing) cts_filtered = cts_no_outlier[genesPassing,] print( paste( 'total number of genes' , totalNbGenes )) print( paste( 'number of genes with at least 1 read' , NbGenesNonEmpty )) print( paste( 'number of genes passing the filter' , nbGenesPassing )) dim(cts_filtered) #head(cts_filtered) ``` ### 2.2. DESeq ```{r} coldata_no_outlier$genotype_time = as.factor( paste( coldata_no_outlier$genotype , coldata_no_outlier$time_point, sep='.') ) # design should be a formula using the column of the info table. # In order to benefit from the default settings of the package, # you should put the variable of interest at the end of the formula # and make sure the control level is the first level. RNAdata <- DESeqDataSetFromMatrix(countData = cts_filtered, colData = coldata_no_outlier, design = ~ replicate + genotype_time) ``` ### PCA ```{r} library(plotly) # variance stabilisation transform RNAdata_VST <- vst(RNAdata, blind=TRUE) ``` ```{r} g = pairedPCAplot(RNAdata_VST, intgroup=c("genotype_time") , pairing = 'mouse' ) + ggtitle("PCA - filtered 10reads 3samples - VST") ggplotly(g) ``` The first axis (44%) seems related to genotype, while the second (19%) relates to time points, with the appearance of a trajectory (increase in axis 2 from T0 to T12, then decrease from T12 to T20 to come back "close" to T4/T0). ```{r} getNtopPCA = function (object, ntop = 500) { rv <- rowVars(assay(object)) select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))] pca <- prcomp(t(assay(object)[select, ])) return(pca) } pca = getNtopPCA( RNAdata_VST ) percentVar <- pca$sdev^2/sum(pca$sdev^2) plot(percentVar[1:10] , main='PCA' , xlab='PCA axis' , ylab="percent of variance") ``` ```{r} a1=2 a2=3 pcaAxis = data.frame( pca$x ) ggplot(data=pcaAxis, aes_string(x = paste0("PC",a1), y = paste0("PC",a2), color = coldata_no_outlier$time_point , pch= coldata_no_outlier$genotype)) + geom_point(size = 3) + xlab(paste0("PC",a1,": ", round(percentVar[a1] * 100), "% variance")) + ylab(paste0("PC",a2,": ", round(percentVar[a2] * 100), "% variance")) + coord_fixed() ``` PC3 seems to differentiate T4,T8 ("morning") from T16,T20 ("afternoon"). Together, PC2 and PC3 exhibit the so-called cycle. ```{r} a1=2 a2=4 pcaAxis = data.frame( pca$x ) ggplot(data=pcaAxis, aes_string(x = paste0("PC",a1), y = paste0("PC",a2), color = coldata_no_outlier$time_point , pch= coldata_no_outlier$genotype)) + geom_point(size = 3) + xlab(paste0("PC",a1,": ", round(percentVar[a1] * 100), "% variance")) + ylab(paste0("PC",a2,": ", round(percentVar[a2] * 100), "% variance")) + coord_fixed() ``` It is not immediately obvious what PC4 relates to. PC1, 2, and 3 shall make a very nice 3D plot : ```{r} fig <- plot_ly(pcaAxis, x = ~PC1, y = ~PC2, z = ~PC3, color = coldata_no_outlier$time_point , symbol= coldata_no_outlier$genotype , symbols=c('x','circle')) fig <- fig %>% add_markers() fig <- fig %>% layout(scene = list(xaxis = list(title = paste0("PC1: ", round(percentVar[1] * 100), "% variance")), yaxis = list(title = paste0("PC1: ", round(percentVar[2] * 100), "% variance")), zaxis = list(title = paste0("PC1: ", round(percentVar[3] * 100), "% variance")))) fig ``` ### 2.3. plotting the counts of different marker genes / genes of interest potential marker genes : * Bmal1 ENSMUSG00000055116 * RORĪ± ENSMUSG00000032238 * RORb ENSMUSG00000036192 * Cry1 ENSMUSG00000020038 * Cry2 ENSMUSG00000068742 ```{r , fig.height=12} intGenes = c( 'Bmal1' ='ENSMUSG00000055116', 'RORĪ±' = 'ENSMUSG00000032238', 'RORb' = 'ENSMUSG00000036192', 'Cry1' = 'ENSMUSG00000020038', 'Cry2' = 'ENSMUSG00000068742') ps = list() for(i in 1:length(intGenes)){ print(paste(i,intGenes[i], names(intGenes)[i])) if( !( intGenes[i] %in% rownames( cts_filtered ) ) ){ next } x = plotCounts(RNAdata, gene=intGenes[i], intgroup=c('time_point'), returnData=TRUE) y = ggplot(x, aes(x=time_point, y=count , color = coldata_no_outlier$genotype)) + geom_point(position=position_jitter(w=0.1,h=0)) + scale_y_log10() + ggtitle(names(intGenes)[i]) ps[[names(intGenes)[i]]] = y#ps = c(ps , y) } grid.arrange( grobs=ps, nrow=length(ps) ) ``` ### 3. DE analysis ```{r} RNAdata <- DESeq(RNAdata) ``` ```{r} # checking for outliers par(mar=c(8,5,2,2)) boxplot(log10(assays(RNAdata)[["cooks"]]), range=0, las=2, main="Boxplot of Cook's distance - filtered 10reads 3samples ") ``` No samples appear to be systematically above the others. ```{r} # check dispersion estimates plotDispEsts(RNAdata, main="dispersion estimates - filtered 10reads 3samples") ``` Seems fairly usual. ```{r} resultsNames(RNAdata) ``` ```{r} save(RNAdata , file='deseq2results.replicate_genotype_time.Rdata') load('deseq2results.replicate_genotype_time.Rdata') ``` Using a **FDRthreshold of 0.01**. ### contrast 1 : effect of KO at different time points ```{r} FDRthreshold = 0.01 for( T in c(0,4,8,12,16,20)){ res <- results(RNAdata, contrast=c("genotype_time", paste0("KO.T",T), paste0("WT.T",T)),alpha = FDRthreshold ) print(summary(res)) resNorm <- lfcShrink(RNAdata, res=res, type='ashr') DESeq2::plotMA(resNorm, alpha = FDRthreshold, main=paste( "effect of ROR-alpha at time ", T , "\n- filtered 10reads 4samples"), ylim=c(-0.5,0.5)) #replacing NA p-values by 1s resNorm$padj <- ifelse(is.na(resNorm$padj), 1, resNorm$padj) orderedResNorm <- resNorm[ order( resNorm$padj ) , ] orderedResNorm$SYMBOL = entrezId[ rownames(orderedResNorm) , 'SYMBOL' ] write.table( orderedResNorm , paste0(baseDir,"/",paste("RORalpha_KO_vs_WT",paste0("T",T),"filter_10reads_3samples","replicate.genotype_type","all","txt" , sep='.'))) } ``` ### cohesiveness of KO effect across time I will perform a simple exploration by counting the genes which are DE in different combinations of time points First upregulated genes ```{r} deSets = data.frame( "geneId" = rownames(cts_filtered) ) rownames( deSets ) = deSets$geneId for( T in c(0,4,8,12,16,20)){ res = read.table( paste0(baseDir,"/", paste("RORalpha_KO_vs_WT",paste0("T",T), "filter_10reads_3samples", "replicate.genotype_type","all","txt" , sep='.') ), header=TRUE ) deSets[paste0("T",T)] = 0 sig = rownames(res)[ ( res$padj<0.01 )&(res$log2FoldChange>0) ] deSets[sig, paste0("T",T)] = 1 } deSets = deSets[ , -1] print( colSums(deSets) ) library(UpSetR) upset(deSets , nsets=6 , sets = c('T0','T4','T8','T12','T16','T20'), keep.order=TRUE, mainbar.y.label = "KO.vs.WT - upregulated genes Intersections", sets.x.label = "time point", order.by='freq') ``` ```{r} common = apply(deSets , 1 , sum) == 6 write( rownames( deSets )[common] , paste0(baseDir,"/", paste("RORalpha_KO_vs_WT","COMMON_ALL_TIME","UP", "filter_10reads_3samples", "replicate.genotype_type","txt" , sep='.') ) ) onlyOne = apply(deSets , 1 , sum) == 1 for( T in c(0,4,8,12,16,20) ){ cName = paste0('T',T) m = ( deSets[cName] == 1 ) & onlyOne write( rownames( deSets )[m] , paste0(baseDir,"/", paste("RORalpha_KO_vs_WT",paste0("UNIQUE_",cName),"UP", "filter_10reads_3samples", "replicate.genotype_type","txt" , sep='.') ) ) } ``` Now downregulated genes ```{r} deSets = data.frame( "geneId" = rownames(cts_filtered) ) rownames( deSets ) = deSets$geneId for( T in c(0,4,8,12,16,20)){ res = read.table( paste0(baseDir,"/", paste("RORalpha_KO_vs_WT",paste0("T",T), "filter_10reads_3samples", "replicate.genotype_type","all","txt" , sep='.') ), header=TRUE ) deSets[paste0("T",T)] = 0 sig = rownames(res)[ ( res$padj<0.01 )&(res$log2FoldChange<0) ] deSets[sig, paste0("T",T)] = 1 } deSets = deSets[ , -1] print( colSums(deSets) ) library(UpSetR) upset(deSets , nsets=6 , sets = c('T0','T4','T8','T12','T16','T20'), keep.order=TRUE, mainbar.y.label = "KO.vs.WT - down genes Intersections", sets.x.label = "time point", order.by='freq') ``` ```{r} common = apply(deSets , 1 , sum) == 6 write( rownames( deSets )[common] , paste0(baseDir,"/", paste("RORalpha_KO_vs_WT","COMMON_ALL_TIME","DOWN", "filter_10reads_3samples", "replicate.genotype_type","txt" , sep='.') ) ) onlyOne = apply(deSets , 1 , sum) == 1 for( T in c(0,4,8,12,16,20) ){ cName = paste0('T',T) m = ( deSets[cName] == 1 ) & onlyOne write( rownames( deSets )[m] , paste0(baseDir,"/", paste("RORalpha_KO_vs_WT",paste0("UNIQUE_",cName),"DOWN", "filter_10reads_3samples", "replicate.genotype_type","txt" , sep='.') ) ) } ``` ## 4. downstream analysis reference : http://yulab-smu.top/clusterProfiler-book/ , https://yulab-smu.top/biomedical-knowledge-mining-book/reactomepa.html ```{r} library(clusterProfiler) library(ReactomePA) entrezId = bitr( rownames( cts_filtered ) , fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Mm.eg.db") ## problem -> n:n mapping from ensembl to entrez ## -> for now remove all duplicates on both sides nbefore=nrow(entrezId) entrezId = entrezId[!duplicated(entrezId$ENSEMBL),] entrezId = entrezId[!duplicated(entrezId$ENTREZID),] nafter=nrow(entrezId) print(paste("removed",nbefore-nafter,"multiplicates")) print(paste( ' genes with a n Entrez id', length( entrezId$ENTREZID ) , '/' , length( rownames( cts_filtered ) ) , '(' ,length( entrezId$ENTREZID ) / length( rownames( cts_filtered ) ) , ')' ) ) rownames( entrezId ) = entrezId$ENSEMBL ``` ```{r} ## the label_format argument does not work as intended... I had to grab the function code from https://rdrr.io/bioc/enrichplot/src/R/dotplot.R ## and correct it... myDotplot = function( object , x = "GeneRatio", color = "p.adjust", showCategory=10, split = NULL,font.size=12,title = "",orderBy="x",decreasing=TRUE,label_format = 30,size = "Count"){ colorBy <-"p.adjust" df <- fortify(object, showCategory = showCategory, split=split) if (orderBy != 'x' && !orderBy %in% colnames(df)) { message('wrong orderBy parameter; set to default `orderBy = "x"`') orderBy <- "x" } if (orderBy == "x") { df <- dplyr::mutate(df, x = eval(parse(text=x))) } subLabel = function(s){ if( str_length(s)