--- title: "12-Apul-miRNA-mRNA-WGCNA" author: "Kathleen Durkin" date: "2025-01-01" output: bookdown::html_document2: theme: cosmo toc: true toc_float: true number_sections: true code_folding: show code_download: true html_document: theme: cosmo toc: true toc_float: true number_sections: true code_folding: show code_download: true github_document: toc: true number_sections: true bibliography: references.bib --- To-Do/Questions: - why am I getting no significant module correlation with timepoint, but clear PCA clustering by timepoint? - Correct formatting of WGCNA input data. This code currently performs count + pOverA filtering separately on raw gene counts and raw miRNA counts (filtered from all sRNA counts)), then combines for DESeq2 variance stabilization and WGCNA. I think it'd be better to separately variance stabilize counts since the two data sets were generated using different sequencing and alignment methods. Also think, if doing separate variance stabilization, may need tovariance stabilize full sRNA dataset before filtering down to miRNAs, since the miRNAs are a tiny subset of the full sRNA data. - How to make useful visualizations in Cytoscvape? Was able to import file, but had trouble working with the network Running Weighted Gene Correlation Network Analysis (WGCNA) to assess patterns of miRNA-mRNA coexpression in A.pulchra. ```{r setup, include=FALSE} library(knitr) knitr::opts_chunk$set( echo = TRUE, # Display code chunks eval = TRUE, # Evaluate code chunks warning = FALSE, # Hide warnings message = FALSE, # Hide messages comment = "" # Prevents appending '##' to beginning of lines in code output ) ``` # Install and load packages ```{r load_libraries, inlcude = TRUE} library(tidyverse) library(ggplot2) library(WGCNA) library(magrittr) library(genefilter) library(DESeq2) library(ggfortify) library(RColorBrewer) library(pheatmap) library(factoextra) library(vegan) library(dendsort) library(ComplexHeatmap) ``` # Load data Load in count matrices for RNAseq. ```{r load-data} # raw gene counts data (will filter and variance stabilize) Apul_genes <- read_csv("../output/02.20-D-Apul-RNAseq-alignment-HiSat2/apul-gene_count_matrix.csv") Apul_genes <- as.data.frame(Apul_genes) # format gene IDs as rownames (instead of a column) rownames(Apul_genes) <- Apul_genes$gene_id Apul_genes <- Apul_genes%>%select(!gene_id) # load and format metadata metadata <- read_csv("../../M-multi-species/data/rna_metadata.csv")%>%select(AzentaSampleName, ColonyID, Timepoint)%>% filter(grepl("ACR", ColonyID)) metadata$Sample <- paste(metadata$AzentaSampleName, metadata$ColonyID, metadata$Timepoint, sep = "_") colonies <- unique(metadata$ColonyID) # Load physiological data phys<-read_csv("https://github.com/urol-e5/timeseries/raw/refs/heads/master/time_series_analysis/Output/master_timeseries.csv")%>%filter(colony_id_corr %in% colonies)%>% select(colony_id_corr, species, timepoint, site, Host_AFDW.mg.cm2, Sym_AFDW.mg.cm2, Am, AQY, Rd, Ik, Ic, calc.umol.cm2.hr, cells.mgAFDW, prot_mg.mgafdw, Ratio_AFDW.mg.cm2, Total_Chl, Total_Chl_cell, cre.umol.mgafdw) # format timepoint phys$timepoint <- gsub("timepoint", "TP", phys$timepoint) #add column with full sample info phys <- merge(phys, metadata, by.x = c("colony_id_corr", "timepoint"), by.y = c("ColonyID", "Timepoint")) %>% select(-AzentaSampleName) #add site information into metadata metadata$Site<-phys$site[match(metadata$ColonyID, phys$colony_id_corr)] # Rename gene column names to include full sample info (as in miRNA table) colnames(Apul_genes) <- metadata$Sample[match(colnames(Apul_genes), metadata$AzentaSampleName)] # raw miRNA counts (will filter and variance stabilize) Apul_miRNA <- read.table(file = "../output/03.10-D-Apul-sRNAseq-expression-DESeq2/Apul_miRNA_ShortStack_counts_formatted.txt", header = TRUE, sep = "\t", check.names = FALSE) ``` Ensure there are no genes or miRNAs with 0 counts across all samples. ```{r} nrow(Apul_genes) Apul_genes<-Apul_genes %>% mutate(Total = rowSums(.[, 1:40]))%>% filter(!Total==0)%>% dplyr::select(!Total) nrow(Apul_genes) # miRNAs nrow(Apul_miRNA) Apul_miRNA<-Apul_miRNA %>% mutate(Total = rowSums(.[, 1:40]))%>% filter(!Total==0)%>% dplyr::select(!Total) nrow(Apul_miRNA) ``` Removing genes with only 0 counts reduced number from 44371 to 35869. Retained all 51 miRNAs. # Physiology filtering Run PCA on physiology data to see if there are phys outliers Export data for PERMANOVA test. ```{r} test<-as.data.frame(phys) test<-test[complete.cases(test), ] ``` Build PERMANOVA model. ```{r} scaled_test <-prcomp(test%>%select(where(is.numeric)), scale=TRUE, center=TRUE) fviz_eig(scaled_test) # scale data vegan <- scale(test%>%select(where(is.numeric))) # PerMANOVA permanova<-adonis2(vegan ~ timepoint*site, data = test, method='eu') permanova ``` ```{r} pca1<-ggplot2::autoplot(scaled_test, data=test, frame.colour="timepoint", loadings=FALSE, colour="timepoint", shape="site", loadings.label.colour="black", loadings.colour="black", loadings.label=FALSE, frame=FALSE, loadings.label.size=5, loadings.label.vjust=-1, size=5) + geom_text(aes(x = PC1, y = PC2, label = paste(colony_id_corr, timepoint)), vjust = -0.5)+ theme_classic()+ theme(legend.text = element_text(size=18), legend.position="right", plot.background = element_blank(), legend.title = element_text(size=18, face="bold"), axis.text = element_text(size=18), axis.title = element_text(size=18, face="bold"));pca1 ``` Remove ACR-173, timepoint 3 sample from analysis. This is Azenta sample 1B2. ```{r} Apul_genes <- Apul_genes%>% select(!`1B2_ACR-173_TP3`) Apul_miRNA <- Apul_miRNA%>% select(!`1B2_ACR-173_TP3`) metadata <- metadata %>% filter(Sample != "1B2_ACR-173_TP3") ``` We also do not have phys data for colony 1B9 ACR-265 at TP4, so I'll remove that here as well. ```{r} Apul_genes <- Apul_genes%>% select(!`1B9_ACR-265_TP4`) Apul_miRNA <- Apul_miRNA%>% select(!`1B9_ACR-265_TP4`) metadata <- metadata %>% filter(Sample != "1B9_ACR-265_TP4") ``` # pOverA filtering *pOverA*: Specifying the minimum count for a proportion of samples for each gene. Here, we are using a pOverA of 0.1. This is because we have 40 samples with a minimum of n=4 samples per timepoint per site. Therefore, we will accept genes that are present in 4/40 = 0.1 of the samples because we expect different expression by life stage. We are further setting the minimum count of genes and miRNA to 10, such that 12.5% of the samples must have a gene count of >10 in order for the gene to remain in the data set. Filter in the package "genefilter". Pre-filtering our dataset to reduce the memory size dataframe, increase the speed of the transformation and testing functions, and improve quality of statistical analysis by removing low-coverage counts. Removed counts could represent outliers in the data and removing these improves sensitivity of statistical tests. genes: ```{r} filt <- filterfun(pOverA(0.1,10)) #create filter for the counts data gfilt <- genefilter(Apul_genes, filt) #identify genes to keep by count filter gkeep <- Apul_genes[gfilt,] #identify genes to keep by count filter gkeep <- Apul_genes[gfilt,] #identify gene lists gn.keep <- rownames(gkeep) #gene count data filtered in PoverA, P percent of the samples have counts over A Apul_genes_filt <- as.data.frame(Apul_genes[which(rownames(Apul_genes) %in% gn.keep),]) #How many rows do we have before and after filtering? nrow(Apul_genes) #Before nrow(Apul_genes_filt) #After ``` We had 35869 genes before, and 23459 genes after filtering. miRNA: ```{r} mifilt <- filterfun(pOverA(0.1,10)) #create filter for the counts data mifilt <- genefilter(Apul_miRNA, mifilt) #identify genes to keep by count filter mikeep <- Apul_miRNA[mifilt,] #identify genes to keep by count filter mikeep <- Apul_miRNA[mifilt,] #identify gene lists mi.keep <- rownames(mikeep) #gene count data filtered in PoverA, P percent of the samples have counts over A Apul_miRNA_filt <- as.data.frame(Apul_miRNA[which(rownames(Apul_miRNA) %in% mi.keep),]) #How many rows do we have before and after filtering? nrow(Apul_miRNA) #Before nrow(Apul_miRNA_filt) #After ``` Of our original 51 miRNAs, 43 are retained. # Combine counts data ```{r} # combine gene counts and miRNA counts # make sure column names in both are in the same order Apul_miRNA_filt <- Apul_miRNA_filt[, colnames(Apul_genes_filt)] # Combine by rows Apul_g_mi_filt <- rbind(Apul_genes_filt, Apul_miRNA_filt) ``` # Assign metadata and arrange order of columns Display current order of metadata and gene count matrix. ```{r} metadata$Sample colnames(Apul_g_mi_filt) ``` Order metadata the same as the column order in the gene matrix. ```{r} list<-colnames(Apul_g_mi_filt) list<-as.factor(list) metadata$Sample<-as.factor(metadata$Sample) # Re-order the levels metadata$Sample <- factor(as.character(metadata$Sample), levels=list) # Re-order the data.frame metadata_ordered <- metadata[order(metadata$Sample),] metadata_ordered$Sample ``` Metadata and gene count matrix are now ordered the same. # Conduct variance stabilized transformation ```{r} #Set DESeq2 design dds <- DESeqDataSetFromMatrix(countData = Apul_g_mi_filt, colData = metadata_ordered, design = ~Timepoint+ColonyID) ``` Check size factors. ```{r} SF.dds <- estimateSizeFactors(dds) #estimate size factors to determine if we can use vst to transform our data. Size factors should be less than 4 for us to use vst print(sizeFactors(SF.dds)) #View size factors all(sizeFactors(SF.dds)) < 4 ``` All size factors are less than 4, so we can use VST transformation. ```{r} vsd <- vst(dds, blind=FALSE) #apply a variance stabilizing transformation to minimize effects of small counts and normalize with respect to library size head(assay(vsd), 3) #view transformed gene count data for the first three genes in the dataset. ``` # PCA Keep in mind this is of both the genes and miRNAs. since there are vastly more genes than miRNAs, the PCA will be dictated by genes. for miRNas alone, see `03.10`. ```{r} plotPCA(vsd, intgroup = c("ColonyID")) plotPCA(vsd, intgroup = c("Timepoint")) ``` # Sample clustering ```{r} sampleDists <- dist(t(assay(vsd))) #calculate distance matrix sampleDistMatrix <- as.matrix(sampleDists) #distance matrix rownames(sampleDistMatrix) <- colnames(vsd) #assign row names colnames(sampleDistMatrix) <- NULL #assign col names colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255) #assign colors pht <- pheatmap(sampleDistMatrix, #plot matrix clustering_distance_rows=sampleDists, #cluster rows clustering_distance_cols=sampleDists, #cluster columns col=colors) #set colors print(pht) ``` # Outlier checks Look for outliers by examining tree of samples ```{r} datExpr <- as.data.frame(t(assay(vsd))) sampleTree = hclust(dist(datExpr), method = "average"); plot(sampleTree, main = "Sample clustering to detect outliers: genes", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2) ``` Soft threshold selection. The soft thresholding power (β) is the number to which the co-expression similarity is raised to calculate adjacency. The function pickSoftThreshold performs a network topology analysis. The user chooses a set of candidate powers, however the default parameters are suitable values. ```{r, message=FALSE, warning=FALSE} allowWGCNAThreads() # # Choose a set of soft-thresholding powers powers <- c(seq(from = 1, to=19, by=2), c(21:30)) #Create a string of numbers from 1 through 10, and even numbers from 10 through 20 # # # Call the network topology analysis function sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 5) ``` Plot the results. ```{r} sizeGrWindow(9, 5) par(mfrow = c(1,2)); cex1 = 0.9; # # # Scale-free topology fit index as a function of the soft-thresholding power plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence")); text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=cex1,col="red"); # # # this line corresponds to using an R^2 cut-off abline(h=0.9,col="red") # # # Mean connectivity as a function of the soft-thresholding power 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=cex1,col="red") ``` I used a scale-free topology fit index **R^2 of 0.9**. The lowest recommended R^2 by Langfelder and Horvath is 0.8. I chose 0.9 because we want to use the smallest soft thresholding power that maximizes with model fit. It appears that our **soft thresholding power is 5** because it is the lowest power above the R^2=0.9 threshold that maximizes with model fit. I will use a **signed network**. # WGCNA Now we're ready to run WGCNA! ```{r, cache=TRUE} picked_power = 5 temp_cor <- cor cor <- WGCNA::cor # Force it to use WGCNA cor function (fix a namespace conflict issue) netwk_Apul <- blockwiseModules(datExpr, # <= input here # == Adjacency Function == power = picked_power, # <= power here networkType = "signed", # == Tree and Block Options == deepSplit = 2, pamRespectsDendro = F, # detectCutHeight = 0.75, minModuleSize = 30, maxBlockSize = 5000, # == Module Adjustments == reassignThreshold = 1e-6, mergeCutHeight = 0.15, # == TOM == Archive the run results in TOM file (saves time) saveTOMs = F, saveTOMFileBase = "ER", # == Output Options numericLabels = T, verbose = 3) cor <- temp_cor # Return cor function to original namespace ``` Take a look at dendrogram. ```{r} # Convert labels to colors for plotting mergedColors = labels2colors(netwk_Apul$colors) labels = table(netwk_Apul$colors) labels # Plot the dendrogram and the module colors underneath plotDendroAndColors( netwk_Apul$dendrograms[[1]], mergedColors[netwk_Apul$blockGenes[[1]]], "Module colors", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05 ) ``` Show module number before and after merging. ```{r} length(table(netwk_Apul$unmergedColors)) length(table(netwk_Apul$colors)) ``` There are 54 modules before merging to 85% similarity and 52 after. ```{r} MEs<-netwk_Apul$MEs moduleLabels<-netwk_Apul$colors moduleColors <- labels2colors(netwk_Apul$colors) ``` ```{r} # Get Module Eigengenes per cluster MEs0_Apul <- moduleEigengenes(datExpr, mergedColors)$eigengenes # Add treatment names MEs0_Apul$sample = row.names(MEs0_Apul) # Join metadata to add timepoint to mME_Apul mME_Apul <- mME_Apul %>% left_join(metadata, by = c("sample" = "Sample")) # Order samples by timepoint timepoint_order <- c("TP1", "TP2", "TP3", "TP4") # Specify the desired order of timepoints mME_Apul <- mME_Apul %>% mutate(sample = factor(sample, levels = metadata$Sample[order(match(metadata$Timepoint, timepoint_order))])) # Plot heatmap mME_Apul %>% ggplot(., aes(x = sample, y = name, fill = value)) + geom_tile() + theme_bw() + scale_fill_gradient2( low = "blue", high = "red", mid = "white", midpoint = 0, limit = c(-1,1) ) + theme(axis.text.x = element_text(angle = 90)) + labs(title = "Module-trait Relationships", y = "Modules", fill = "corr") ``` # Correlate to traits Ensure we've removed sample from expression data that is not present in physiological data. 1B9 ```{r} phys<-phys%>% filter(!Sample=="1B2_ACR-173_TP3")%>% filter(!Sample=="1B9_ACR-265_TP4") setdiff(metadata_ordered$Sample, phys$Sample) #all have 38 observations ``` Create a physiological/factor dataset with sample in rows and phys data in columns (row name = Azenta sample name) ```{r} # phys$AzentaSampleName<-metadata_ordered$AzentaSampleName[match(phys$Sample, metadata_ordered$Sample)] traits<-phys%>% select(Sample, timepoint, site, Host_AFDW.mg.cm2, Sym_AFDW.mg.cm2, Am, Rd, Ik, Ic, AQY, cells.mgAFDW, Total_Chl, Ratio_AFDW.mg.cm2, Total_Chl_cell) #time point #site traits<-traits%>% mutate(site=gsub("Mahana", "1", site)) %>% mutate(site=gsub("Manava", "2", site))%>% mutate(timepoint=gsub("TP", "", timepoint)) traits$timepoint<-as.numeric(traits$timepoint) traits$site<-as.numeric(traits$site) traits<-as.data.frame(traits) rownames(traits)<-traits$Sample datTraits<-traits%>%select(!Sample) #datTraits<-datTraits%>%select(!Timepoint) ``` From this tutorial: https://alexslemonade.github.io/refinebio-examples/04-advanced-topics/network-analysis_rnaseq_01_wgcna.html#410_Which_modules_have_biggest_differences_across_treatment_groups ```{r} # Create the design matrix from the `time_point` variable des_mat <- model.matrix(~ datTraits$site) # lmFit() needs a transposed version of the matrix fit <- limma::lmFit(t(MEs), design = des_mat) # Apply empirical Bayes to smooth standard errors fit <- limma::eBayes(fit) stats_df <- limma::topTable(fit, number = ncol(MEs)) %>% tibble::rownames_to_column("module") head(stats_df) ``` Add in temperature and light information for each time point from physiology manuscript. ```{r} # env<-read_csv("https://github.com/urol-e5/timeseries/raw/refs/heads/master/time_series_analysis/Output/environment_characteristics_RDA.csv")%>%select(mean_Temp_mean, mean_solar_rad_kwpm2_mean, timepoint)%>% # mutate(timepoint=if_else(timepoint=="timepoint1", "1", # if_else(timepoint=="timepoint2", "2", # if_else(timepoint=="timepont3", "3", "4"))))%>% # mutate(timepoint=as.numeric(timepoint))%>% # select(!timepoint) # # #merge into datTraits # # datTraits<-left_join(datTraits, env) # rownames(datTraits)<-traits$Sample # # str(datTraits) # # datTraits<-datTraits%>%select(!timepoint) ``` ```{r} nGenes = ncol(datExpr) nSamples = nrow(datExpr) nGenes nSamples ``` Generate labels for module eigengenes as numbers. ```{r} MEs0 = moduleEigengenes(datExpr, moduleLabels, softPower=5)$eigengenes MEs = orderMEs(MEs0) names(MEs) ``` Correlations of traits with eigengenes ```{r} moduleTraitCor = cor(MEs, datTraits, use = "p"); moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples); Colors=sub("ME","", names(MEs)) moduleTraitTree = hclust(dist(t(moduleTraitCor)), method = "average") # pdf(file="D-Apul/output/11.00-Apul-WGCNA/ModuleTraitClusterTree.pdf") plot(moduleTraitTree) # dev.off() ``` Correlations of genes with eigengenes. Calculate correlations between ME's and lifestages. ```{r} moduleGeneCor=cor(MEs,datExpr) moduleGenePvalue = corPvalueStudent(moduleGeneCor, nSamples); ``` Calculate kME values (module membership). ```{r} datKME = signedKME(datExpr, MEs, outputColumnName = "kME") head(datKME) ``` Generate a complex heatmap of module-trait relationships. ```{r} #bold sig p-values #dendrogram with WGCNA MEtree cut-off #colored y-axis #Create list of pvalues for eigengene correlation with specific life stages heatmappval <- signif(moduleTraitPvalue, 1) #Make list of heatmap row colors htmap.colors <- names(MEs) htmap.colors <- gsub("ME", "", htmap.colors) row_dend = dendsort(hclust(dist(moduleTraitCor))) col_dend = dendsort(hclust(dist(t(moduleTraitCor)))) pdf(file = "../output/12-Apul-miRNA-mRNA-WGCNA/Module-trait-relationship-heatmap.pdf", height = 14, width = 8) Heatmap(moduleTraitCor, name = "Eigengene", row_title = "Gene Module", column_title = "Module-Trait Eigengene Correlation", col = blueWhiteRed(50), row_names_side = "left", #row_dend_side = "left", width = unit(5, "in"), height = unit(8.5, "in"), #column_dend_reorder = TRUE, #cluster_columns = col_dend, row_dend_reorder = TRUE, #column_split = 6, row_split=3, #column_dend_height = unit(.5, "in"), #column_order = lifestage_order, cluster_rows = row_dend, row_gap = unit(2.5, "mm"), border = TRUE, cell_fun = function(j, i, x, y, w, h, col) { if(heatmappval[i, j] < 0.05) { grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 10, fontface = "bold")) } else { grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 10, fontface = "plain")) }}, column_names_gp = gpar(fontsize = 12, border=FALSE), column_names_rot = 35, row_names_gp = gpar(fontsize = 12, alpha = 0.75, border = FALSE)) #draw(ht) #dev.off() ``` This plot is cool, but I'm a little skeptical of the results. No modules show significant relationship with timepoint, despite the PCAs showing a pretty clear clustering by timepoint... # Plot eigengene values View module eigengene data and make dataframe for Strader plots. ```{r} head(MEs) names(MEs) Eigen_MEs <- MEs Eigen_MEs$AzentaSampleName <- rownames(Eigen_MEs) head(Eigen_MEs) plotTraits<-datTraits plotTraits$AzentaSampleName <- rownames(plotTraits) Eigen_MEs<-Eigen_MEs%>% droplevels() #drop unused level dim(Eigen_MEs) head(Eigen_MEs) ``` Plot mean module eigengene for each module. ```{r} #convert wide format to long format for plotting str(Eigen_MEs) plot_MEs <- Eigen_MEs %>% pivot_longer( cols = where(is.numeric), # Select only numeric columns names_to = "Module", # Name for the new column containing the column names values_to = "Mean" # Name for the new column containing the values ) str(plot_MEs) #join with phys data plot_MEs<-left_join(plot_MEs, plotTraits) expression_plots<-plot_MEs%>% group_by(Module) %>% ggplot(aes(x=mean_Temp_mean, y=Mean)) + facet_wrap(~ Module)+ geom_point()+ geom_smooth(method="lm")+ #ylim(-0.5,1) + geom_hline(yintercept = 0, linetype="dashed", color = "grey")+ theme_bw() + theme(axis.text.x=element_text(angle = 45, hjust=1, size = 12), #set x-axis label size axis.title.x=element_text(size = 14), #set x-axis title size axis.ticks.x=element_blank(), #No x-label ticks #axis.title.y=element_blank(), #No y-axis title axis.text.y=element_text(size = 14), #set y-axis label size, panel.border = element_rect(color = "black", fill = NA, size = 1), #set border panel.grid.major = element_blank(), #Set major gridlines panel.grid.minor = element_blank(), #Set minor gridlines axis.line = element_line(colour = "black"), #Set axes color plot.background=element_blank(), plot.title = element_text(size=22)); expression_plots ``` # Export to Cytoscape Export modules of interest for network visualization ```{r} mods_interest <- c("26", "9", "2", "50", "27", "10", "5", "38", "16", "14", "11", "6", "1", "49", "8", "31", "10", "39", "41", "51", "29", "11", "45", "21") ``` Note that, due to the size of many modules, it's probably most practical to export modules individually (unless interested in in connections between modules) ```{r} # Recalculate topological overlap if needed TOM = TOMsimilarityFromExpr(datExpr, power = 5); # Read in the annotation file annot = datTraits; # Select modules modules = c("50"); # Select module probes probes = names(datExpr) inModule = is.finite(match(moduleLabels, modules)); modProbes = probes[inModule]; # modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)]; # Select the corresponding Topological Overlap modTOM = TOM[inModule, inModule]; dimnames(modTOM) = list(modProbes, modProbes) # Export the network into edge and node list files Cytoscape can read cyt = exportNetworkToCytoscape(modTOM, edgeFile = paste("CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""), nodeFile = paste("CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""), weighted = TRUE, threshold = 0.02, nodeNames = modProbes, # altNodeNames = modGenes, nodeAttr = moduleLabels[inModule]); ``` ```{r, engine='bash'} # Move cytoscape files to output directory mv ./CytoscapeInput* ../output/12-Apul-miRNA-mRNA-WGCNA ``` To load network to Cytoscape: 1. File -> Import -> Network from File... 2. Select Edges file (created above) 3. Assign fromNode to "Source Node," toNode to "Target Node," weight to "Edge Attribute," and direction to "Interaction Type" 4. Once netwrok is loaded, can manipulate/style as desired. ```{r} # Print session info sessioninfo::session_info() ```