notebook curso

mine.v2.Rmd

+---
+title: "R Notebook"
+output: html_notebook
+---
+
+
+Vamos a comenzar cargando los paquetes necesarios
+```{r echo=FALSE, message=FALSE, warning=FALSE}
+require(edgeR)
+require(WGCNA)
+require(flashClust)
+require(org.At.tair.db)
+
+#para que el calculo wgcna sea multicore
+enableWGCNAThreads() 
+
+```
+
+Esta notebook esta pensada para ilustrar el tipo de analisis que puede realizarse utilizando los conceptos introducido por Horvath y colaboradores en relacion al uso de Redes pesadas de correlacion (Weighted Gene Correlation Networks) para la caracterizacion y analisis de estados transcripcionales.
+
+Para ilustrar los conceptos introducidos haremos uso de los datos RNAseq (GSE88798) presentados en el trabajo de Mine y colaboradores [Plant Cell 2018](http://www.plantcell.org/content/30/6/1199). El paper presenta un estudio transcripcional sobre los mecanismos de respuesta inmune en plantas (Arabidospis thaliana). Para ello analiza la dinamica transcripcional de plantas salvajes y knowouts de diferentes componentes moleculares involucrados en dos tipos de respuesta: pattern-triggered immunity (PTY) y effector-triggered immunity (ETI). En ambos casos se realizan ensayos utilizando alternativamente placebo, cepas virulentas y no-virulentas sobre distintos background geneticos. Nosotros nos concentraremos en los experimentos realizados sobre plantas wild-type
+
+
+Aqui esta el abstract del trabajo:
+
+###### *The phytohormone network consisting of jasmonate, ethylene, PHYTOALEXIN-DEFICIENT4, and salicylic acid signaling is
+required for the two modes of plant immunity, pattern-triggered immunity (PTI), and effector-triggered immunity (ETI). A
+previous study showed that during PTI, the transcriptional responses of over 5000 genes qualitatively depend on complex
+interactions between the network components. However, the role of the network in transcriptional reprogramming during
+ETI and whether it differs between PTI and ETI remain elusive. Here, we generated time-series RNA-sequencing data of
+Arabidopsis thaliana wild-type and combinatorial mutant plants deficient in components of the network upon challenge with
+virulent or ETI-triggering avirulent strains of the foliar bacterial pathogen Pseudomonas syringae. Resistant plants such as
+the wild type achieved high-amplitude transcriptional reprogramming 4 h after challenge with avirulent strains and sustained
+this transcriptome response. Strikingly, susceptible plants including the quadruple network mutant showed almost identical
+transcriptome responses to resistant plants but with several hours delay. Furthermore, gene coexpression network structure
+was highly conserved between the wild type and quadruple mutant. Thus, in contrast to PTI, the phytohormone network is
+required only for achieving high-amplitude transcriptional reprogramming within the early time window of ETI against this
+bacterial pathogen.*
+
+
+
+La matriz de expresion que utilizaremos esta disponible en el archivo *GSE88798_ReadCountTable_M001_348.txt.gz* . 
+
+## Carga de datos
+```{r}
+data <- read.table(file="data/GSE88798_ReadCountTable_M001_348.txt.gz",header=TRUE,stringsAsFactors=FALSE)
+
+#decodificamos los nombres de las samples para entender los fenotipos que disponemos
+samples           <- strsplit2(colnames(data),"_")
+colnames(samples) <- c("genotype_id","genotype","treatment_id","treatment","time","rep")
+
+```
+
+Entendemos que cargamos...
+```{r}
+#cuantas muestras tengo?
+dim(samples)
+
+#a ver ...
+head(samples)
+
+#cuantas muestras tengo de cada genotipo?
+table(samples[,"genotype"])
+
+#cuantas muestras tengo de cada tratamiento?
+table(samples[,"treatment"])
+
+#tabla de doble entrada: cuantos tratamientos sobre cada genotipo
+table(samples[,"genotype"],samples[,"treatment"])
+
+#armo una variable genotype_treatment, por ejemplo: Col_mock, Col_EV ,...
+gt <- apply(samples,1,function(x){paste(x[2],x[4],sep="_")})
+table(gt)
+
+#Cuantas replicas por punto temporal tengo para cada combinacion genotype_treatment 
+table(gt,samples[,"time"])
+```
+
+## Estadistica de las cuentas. 
+
+* cuantas reads hay por muestra (profundidad de la muestra
+* cuantas reads recibe cada gen a lo largo de todas las muestras
+```{r}
+#cuento cuantas reads tiene cada muestra (cada columna de dataCol)
+countsPerSample <- apply(data,2,sum)
+hist(countsPerSample,breaks=20)
+
+#Entre la muestra con mas profundidad y menos profundidad hay un factor *~3x*
+max(countsPerSample)/min(countsPerSample)
+
+#cuantos reads recibe cada gen a lo largo de todas las muestras
+# (linea punteada: 10*nro de samples)
+countsPerGene <- apply(data,1,sum)
+plot(density(log10(countsPerGene)),xlab="log10(#reads)")
+abline(v=log10(10*ncol(data)),col="gray",lty=2)
+```
+
+
+## Normalizacion y procesamiento de counts
+En esta parte nos focalizamos en el analisis de plantas salvajes sometidas a los diferentes tratamientos: *mock* (placebo), *AvrRpt2*, *AvrRpm1* (cepas avirulentas Pto DC3000, asociadas a respuesta inmune tipo *ETI* (effector-triggered immunity)) , y *EV* (virulent Pto DC3000) 
+
+```{r}
+#Nos centramos en el analisis de Col + mock/AvrRpt2 
+#Seteo filtros para elegir las samples de interes
+#y genero indice que las selecciona
+bGenotype  <- samples[,"genotype"]=="Col"
+bTreatment <- samples[,"treatment"]%in%c("mock","AvrRpt2","AvrRpm1","EV")
+bTime      <- samples[,"time"]%in%c("01h","03h","04h","06h","09h","12h","16h","24h")
+iSamples   <- which( bGenotype & bTreatment & bTime  )
+
+#me quedp con los datos de interes
+data1    <- data[,iSamples]
+samples1 <- samples[iSamples,]
+
+#antes de empezar los reordeno de una forma razonable: tratment:time:rep
+imock <- which(samples1[,"treatment"]=="mock")
+imock <- imock[order(samples1[imock,"time"],samples1[imock,"rep"])]
+iEV <- which(samples1[,"treatment"]=="EV")
+iEV <- iEV[order(samples1[iEV,"time"],samples1[imock,"rep"])]
+iavr1  <- which(samples1[,"treatment"]=="AvrRpm1")
+iavr1  <- iavr1[order(samples1[iavr1,"time"],samples1[iavr1,"rep"])]
+iavr2  <- which(samples1[,"treatment"]=="AvrRpt2")
+iavr2  <- iavr2[order(samples1[iavr2,"time"],samples1[iavr2,"rep"])]
+iordered <- c(imock,iavr1,iavr2,iEV)
+
+#aver como quedo...
+# colnames(data1)[iordered]
+# head(samples1[iordered,])
+
+data1   <- data1[,iordered]
+samples1<- samples1[iordered,]
+
+
+#Para normalizar los datos usamos edgeR + voom (diferente tamanio de libreria)
+#Primero metemos los datos (coutns e informacion sobre samples) en un objeto DGEList
+dge  <- DGEList(counts=data1,samples = samples1[,c("genotype","treatment","time","rep")])
+
+#voy a remover genes que no tengan un nro suficientes de reads
+# (10 reads en promedio por sample)
+minReads <- log2(10*ncol(data1))
+
+keep <- filterByExpr(dge, min.count = 10, min.total.count = minReads)
+table(keep)
+dge  <- dge[keep,,keep.lib.sizes=FALSE]
+
+#ahora voy a calcular los factores norm.factors
+#son necesarios para 'equalizar' las librerias (i.e. muestras) y que sean comparables a pesar de presnetar eventualente
+#diferencias en profundidad (nro de reads generados)
+dge <- calcNormFactors(dge)
+# a ver...
+hist(dge$samples$norm.factors,breaks=20,xlab="samples' norm factors")
+
+
+```
+
+
+
+
+Estimaciones de voom para las ~log2(cpm)
+```{r}
+
+  #uso voom (primero defino grupos experimentales para armar una matriz de disenio que describa la estructura de replicas)
+  group <- apply(dge$samples[,c("genotype","treatment","time")],1,function(x){
+    paste(x,collapse="_")
+  })
+  design <- model.matrix(~0+group)
+  v      <- voom(dge,design)
+  
+   
+  #calculo directo de cpm
+  logCPM <- cpm(dge,log=TRUE,prior.count = 3)
+  #Comparacion logCPM vs Voom: son estimaciones muy similares...diferencias para bajo nro de cuentas
+  plot(logCPM[,11],v$E[,11])
+  abline(0,1)
+  
+   #promediando triplicados para visualizar perfiles
+  xprofiles<-t(apply(v$E,1,function(x){
+                   y <- apply(matrix(x,byrow=FALSE,nrow=3),2,mean) #promedio: pongo como matriz
+                                                                   # (cada fila es un replicado) 
+                                                                   #y luego tomo la media por columnas
+                   return(y)
+                  }))
+  #fijo los nombres de las columnas:
+  treatment <- unique(samples1[,"treatment"])
+  time      <- unique(samples1[,"time"])
+  #convierto 'time' a numeric, sacandole la  'h'
+  time <- as.numeric(substr(time,1,2))
+  #ahora si...pongo los nombres correctos a las columnas
+  colnames(xprofiles)<-paste(rep(treatment,each=length(time)),time,sep="_")
+ } 
+
+```
+
+## Analisis WGCNA - 1
+El workflow de WGCNA involucra:
+0. Data cleaning
+1. prefiltrado de genes informativos
+2. calcular correlaciones entre todos los pares de genes
+3. Elegir el valor del exponente,$p$, que se usara para estimar pesos a partir de correlaciones: $w_{ij}=cor_{ij}^p$ 
+4.  Estimar similaridad entre genes a partir del calculo de TOM  (Topological overlap measure)
+Esta es una medidad mas robusta que la correlacion original
+5. Detectar clusters
+
+Hasta ahora tenemos los perfiles temporales de los diferentes tratamientos en una unica matriz de expresion (~log2(cpm)).
+En el enfoque WGCNA se estiman redes/modulos para cada tratamiento 
+
+```{r}
+colnames(xprofiles)
+
+#para este tutorial vamos a continuar con un subset de genes informativos a lo largo de los tratamientos considerados
+geneVar <- apply(xprofiles,1,var)
+if(FALSE){ # voy a remover genes que no  varien demasiado (poco informativos)
+  varOK   <- geneVar > 0.2
+  table(varOK)
+}else{     # me quedo con los 3000 que mas varian
+  varOK   <- order(geneVar,decreasing=TRUE)[1:5000]
+}
+```
+
+```{r}
+#vamos a trabajar por cada tratamiento por separado
+# asi que identifico las columnas del tratamiento de interes en la matriz de expresion
+troi <- treatment[1]
+itr <- grep(troi,colnames(xprofiles))
+
+# Matriz de expresion para analisis WGCNA
+# (wgcna espera que los samples esten como filas y los genes como columnas (!))
+datExpr <- t(xprofiles[varOK,itr])
+
+# Busqueda de soft-threshold (nos nos estressamos tanto por esto...)
+powers <- c(c(1:10), seq(from = 12, to=20, by=2));
+sft    <- pickSoftThreshold(datExpr,dataIsExpr = TRUE,powerVector = powers,corFnc = cor,corOptions = list(use = 'p'),networkType = "signed")
+
+#
+layout(matrix(1:2,1,2))
+plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],typ="n")
+text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=0.9,col="red")
+abline(h=0.80,col="red")
+
+plot(sft$fitIndices[,1], sft$fitIndices[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main = paste("Mean connectivity"))
+text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=0.9,col="red")
+
+
+softPower = 18;
+type <- "signed"
+
+#calclute the adjacency matrix
+adj= adjacency(datExpr,type = type, power = softPower);
+
+#turn adjacency matrix into topological overlap to minimize the effects of noise and spurious associations
+dissTOM <- TOMdist(adj,TOMType=type)
+colnames(dissTOM)<-rownames(dissTOM)<-colnames(datExpr)
+
+
+#comparacion de similaridad directa vs TOM
+set.seed(123457)
+isample <- sample(1:10000)
+plot(adj[upper.tri(adj)][isample],
+     1-dissTOM[upper.tri(dissTOM)][isample],pch=20,col=rgb(.4,.4,.4,.2),
+     xlab="adj",ylab="1-disTOM",xlim=c(0,1),ylim=c(0,1))
+abline(0,1,col="gray")
+
+
+# Set the minimum module size
+minModuleSize = 20;
+
+# Module identification using dynamic tree cut
+geneTree    = flashClust(as.dist(dissTOM),method="average")
+
+hyb         = cutreeHybrid(dendro = geneTree,  dissTOM,  minClusterSize = minModuleSize, deepSplit = 1);
+ccut<-hyb$labels
+table(table(ccut))
+
+#ploteemos el dendrograma
+dynamicColors = labels2colors(ccut)
+plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")
+
+
+# Vamos a visualizar los clusters detectados
+# Primero estandarizamos los perfiles para plotear
+zdatExpr <- apply(datExpr,2,function(x){(x-mean(x))/sd(x)})
+
+#elegimos algunos clusters (las etiquetas de menor valor en general corresponden a los clusters mas gdes)
+layout(matrix(1:9,3,3,byrow=TRUE))
+uccut <- unique(ccut)
+uccut <- uccut[uccut!=0]
+for(i in 1:length(uccut)){
+  ig <- which(ccut==i)
+  matplot(time,zdatExpr[,ig],col=rgb(.8,0.8,.8,.1),typ="l",main=paste(troi,"Cluster",i))
+}
+
+# Hay algunos clusters que pueden presentar perfiles parecidos: los mergeo
+# Calculate eigengenes
+MEList = moduleEigengenes(datExpr, colors = dynamicColors)
+MEs    = MEList$eigengenes
+# Calculate dissimilarity of module eigengenes
+MEDiss = 1-cor(MEs);
+# Cluster module eigengenes
+METree = flashClust(as.dist(MEDiss), method = "average");
+
+plot(METree)
+MEDissThres = 0.15
+# Plot the cut line into the dendrogram
+abline(h=MEDissThres, col = "red")
+
+#junto a los que esten a una distancia menor a MEDissThres 
+# Call an automatic merging function
+merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
+# The merged module colors
+mergedColors = merge$colors;
+# Eigengenes of the new merged modules:
+mergedMEs = merge$newMEs;
+
+
+plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
+                    c("Dynamic Tree Cut", "Merged dynamic"),
+                    dendroLabels = FALSE, hang = 0.03,
+                    addGuide = TRUE, guideHang = 0.05)
+
+
+# Vamos a seguir el analisis con los modulos mergeados
+
+# Rename to moduleColors
+moduleColors = mergedColors
+# Construct numerical labels corresponding to the colors
+colorOrder = c("grey", standardColors(50));
+moduleLabels = match(moduleColors, colorOrder)-1
+
+
+
+# GO analysis
+require(GOstats)
+require(annotate)
+require(org.At.tair.db)
+
+
+myGOSummary<-function(hg){
+  if(!exists("myGOTerms"))myGOTerms <- Term(GOTERM)
+
+ 
+  
+  a<-lapply(nodeData(hg@goDag),function(x){
+    c(size     = length(x$geneIds),
+      expCount = x$expCount,
+      oddsRatio = x$oddsRatio,
+      pvalue = x$pvalue)
+  })
+  mysum<-matrix(unlist(a),byrow=TRUE,ncol=4)
+  rownames(mysum)<-names(a)
+  colnames(mysum)<-names(a[[1]])
+  mysum <- data.frame(term=myGOTerms[names(a)],hits=geneCounts(hg)[names(a)],mysum,
+                 fdr=p.adjust(pvalues(hg),method="fdr")[names(a)],stringsAsFactors = FALSE)
+  mysum <- mysum[order(mysum[,"fdr"],decreasing=FALSE),]
+  return(mysum)
+}
+
+
+#recordemos la distribucion de tamanios de los modulos:
+table(moduleColors)
+
+universe <- keys(org.At.tairGO)
+params <- new("GOHyperGParams",
+              geneIds         = universe[1:11],
+              universeGeneIds = universe,
+              annotation      = "org.At.tair.db",
+              ontology        = "BP",
+              pvalueCutoff    = 0.1,
+              conditional     = FALSE,
+              testDirection   ="over")
+
+genes <-colnames(datExpr)[moduleColors %in% modules[i]]
+a     <- rep(1,length(genes))
+names(a)<-genes
+u <- rep(1,length(universe))
+names(u)<-universe
+
+modules <- unique(moduleColors)
+modules <- modules[modules!="grey"]
+moduleGOterm<-c()
+for(i in seq_along(modules)){ 
+  geneIds(params) <- colnames(datExpr)[moduleColors %in% modules[i]]
+  hg              <- hyperGTest(params)
+  
+  a <- myGOSummary(hg)
+
+  #filtro los conceptos GO
+  a <- a[a[,"fdr"] < 0.05 & (a[,"size"]>10 & a[,"size"]<500),]
+  if(nrow(a)>0){
+   plot(-log10(a[,"fdr"]),a[,"hits"])
+   
+   #asigno al cluster el concepto mas enriquecido
+   moduleGOterm <- c(moduleGOterm,a[which.min(a[,"fdr"]),"term"])
+  }else{
+    moduleGOterm <- c(moduleGOterm,NA)
+  }  
+}
+names(moduleGOTerm)<-modules
+
+
+
+#Visualizacion de eigengenes
+MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
+
+simME <- (1+cor(MEs))/2
+simME[simME<0.8]<-0
+
+require(igraph)
+g <- graph_from_adjacency_matrix(simME,mode="undirected",weighted=TRUE,diag=FALSE)
+
+
+plot(g)
+
+
+plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
+= 90)
+
+plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
+plotDendrograms = FALSE, xLabelsAngle = 90)
+
+
+
+
+
+
+```
+
+
+
+
+
+```{r}
+#Relating modules to external traits 
+#(en este caso vamos a definir un trait artificial que en realidad no es independiente: respuesta temprana - media -tardia)
+datTraits <- data.frame(temprano=time<4, media=time>3 & time<9, tardia=time>12)
+
+# Define numbers of genes and samples
+nGenes = ncol(datExpr);
+nSamples = nrow(datExpr);
+# Recalculate MEs with color labels
+MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
+MEs = orderMEs(MEs0)
+moduleTraitCor = cor(MEs, datTraits, use = "p");
+moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
+
+textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
+                   signif(moduleTraitPvalue, 1), ")", sep = "");
+dim(textMatrix) = dim(moduleTraitCor)
+par(mar = c(6, 8.5, 3, 3));
+# Display the correlation values within a heatmap plot
+labeledHeatmap(Matrix = moduleTraitCor,
+               xLabels = names(datTraits),
+               yLabels = names(MEs),
+               ySymbols = names(MEs),
+               colorLabels = FALSE,
+               colors = greenWhiteRed(50),
+               textMatrix = textMatrix,
+               setStdMargins = FALSE,
+               cex.text = 0.5,
+               zlim = c(-1,1),
+               main = paste("Module-trait relationships"))
+
+
+```
+