## ----label= 'Load a sample dataset', eval=TRUE-------------------------------- library(RTN) data(tniData) ## ----label='Create a new TNI object', eval=TRUE, results='hide'--------------- # Input 1: 'expData', a named gene expression matrix (genes on rows, samples on cols); # Input 2: 'regulatoryElements', a vector listing genes regarded as TFs # Input 3: 'rowAnnotation', an optional data frame with gene annotation # Input 4: 'colAnnotation', an optional data frame with sample annotation tfs <- c("FOXM1","E2F2","E2F3","RUNX2","PTTG1") rtni <- tni.constructor(expData = tniData$expData, regulatoryElements = tfs, rowAnnotation = tniData$rowAnnotation, colAnnotation = tniData$colAnnotation) # p.s. alternatively, 'expData' can be a 'SummarizedExperiment' object ## ----label='Permutation', eval=TRUE, results='hide'--------------------------- # Please set nPermutations >= 1000 rtni <- tni.permutation(rtni, nPermutations = 100) ## ----label='Bootstrap', eval=TRUE, results='hide'----------------------------- rtni <- tni.bootstrap(rtni) ## ----label='Run DPI filter', eval=TRUE, results='hide'------------------------ rtni <- tni.dpi.filter(rtni) ## ----label='Check overall summary', eval=TRUE, collapse=TRUE------------------ tni.regulon.summary(rtni) ## ----label='Check summary of a regulon', eval=TRUE, collapse=TRUE------------- tni.regulon.summary(rtni, regulatoryElements = "FOXM1") ## ----label='Check results', eval=TRUE, collapse=TRUE-------------------------- regulons <- tni.get(rtni, what = "regulons.and.mode", idkey = "SYMBOL") head(regulons$FOXM1) ## ----label='Create a new TNA object (preprocess TNI-to-TNA)', eval=TRUE, results='hide'---- # Input 1: 'object', a TNI object with regulons # Input 2: 'phenotype', a named numeric vector, usually log2 differential expression levels # Input 3: 'hits', a character vector, usually a set of differentially expressed genes # Input 4: 'phenoIDs', an optional data frame with gene anottation mapped to the phenotype data(tnaData) rtna <- tni2tna.preprocess(object = rtni, phenotype = tnaData$phenotype, hits = tnaData$hits, phenoIDs = tnaData$phenoIDs) ## ----label='Run the MRA method', eval=TRUE, results='hide'-------------------- # Run the MRA method rtna <- tna.mra(rtna) ## ----label='Get MRA results', eval=TRUE, collapse=TRUE------------------------ # Get MRA results; #..setting 'ntop = -1' will return all results, regardless of a threshold mra <- tna.get(rtna, what="mra", ntop = -1) head(mra) ## ----label='Run GSEA method', eval=TRUE, results='hide'----------------------- # Run the GSEA method # Please set nPermutations >= 1000 rtna <- tna.gsea1(rtna, nPermutations=100) ## ----label='Get GSEA results', eval=TRUE, collapse=TRUE----------------------- # Get GSEA results gsea1 <- tna.get(rtna, what="gsea1", ntop = -1) head(gsea1) ## ----label='Plot GSEA results', eval=FALSE------------------------------------ # # Plot GSEA results # tna.plot.gsea1(rtna, labPheno="abs(log2 fold changes)", ntop = -1) ## ----label='Run the GSEA-2T method', eval=TRUE, results='hide'---------------- # Run the GSEA-2T method # Please set nPermutations >= 1000 rtna <- tna.gsea2(rtna, nPermutations = 100) ## ----label='Get GSEA-2T results', eval=TRUE, collapse=TRUE-------------------- # Get GSEA-2T results gsea2 <- tna.get(rtna, what = "gsea2", ntop = -1) head(gsea2$differential) ## ----label='Plot GSEA-2T results', eval=FALSE--------------------------------- # # Plot GSEA-2T results # tna.plot.gsea2(rtna, labPheno="log2 fold changes", tfs="PTTG1") ## ----eval=FALSE--------------------------------------------------------------- # library(RTN) # library(TCGAbiolinks) # library(SummarizedExperiment) # library(TxDb.Hsapiens.UCSC.hg38.knownGene) # library(snow) ## ----eval=FALSE--------------------------------------------------------------- # # Set GDCquery for the TCGA-BRCA cohort # # Gene expression data will be aligned against hg38 # query <- GDCquery(project = "TCGA-BRCA", # data.category = "Transcriptome Profiling", # data.type = "Gene Expression Quantification", # workflow.type = "HTSeq - FPKM-UQ", # sample.type = c("Primary solid Tumor")) # # # Get a subset for demonstration (n = 500 cases) # cases <- getResults(query, cols = "cases") # cases <- sample(cases, size = 500) # query <- GDCquery(project = "TCGA-BRCA", # data.category = "Transcriptome Profiling", # data.type = "Gene Expression Quantification", # workflow.type = "HTSeq - FPKM-UQ", # sample.type = c("Primary solid Tumor"), # barcode = cases) # GDCdownload(query) # tcgaBRCA_mRNA_data <- GDCprepare(query) ## ----eval=FALSE--------------------------------------------------------------- # # Subset by known gene locations # geneRanges <- genes(TxDb.Hsapiens.UCSC.hg38.knownGene) # tcgaBRCA_mRNA_data <- subsetByOverlaps(tcgaBRCA_mRNA_data, geneRanges) ## ----eval=FALSE--------------------------------------------------------------- # nrow(rowData(tcgaBRCA_mRNA_data)) # ## [1] ~30,000 ## ----eval=FALSE--------------------------------------------------------------- # # Change column names in gene annotation for better summarizations # colnames(rowData(tcgaBRCA_mRNA_data)) <- c("ENSEMBL", "SYMBOL", "OG_ENSEMBL") ## ----eval=FALSE--------------------------------------------------------------- # # Save the preprocessed data for subsequent analyses # save(tcgaBRCA_mRNA_data, file = "tcgaBRCA_mRNA_data_preprocessed.RData") ## ----eval=FALSE--------------------------------------------------------------- # # Load TF annotation # data("tfsData") ## ----eval=FALSE--------------------------------------------------------------- # # Check TF annotation: # # Intersect TFs from Lambert et al. (2018) with gene annotation # # from the TCGA-BRCA cohort # geneannot <- rowData(tcgaBRCA_mRNA_data) # regulatoryElements <- intersect(tfsData$Lambert2018$SYMBOL, geneannot$SYMBOL) ## ----eval=FALSE--------------------------------------------------------------- # # Run the TNI constructor # rtni_tcgaBRCA <- tni.constructor(expData = tcgaBRCA_mRNA_data, # regulatoryElements = regulatoryElements) ## ----eval=FALSE--------------------------------------------------------------- # # Compute the reference regulatory network by permutation and bootstrap analyses. # # Please set 'spec' according to your available hardware # options(cluster=snow::makeCluster(spec=4, "SOCK")) # rtni_tcgaBRCA <- tni.permutation(rtni_tcgaBRCA, pValueCutoff = 1e-7) # rtni_tcgaBRCA <- tni.bootstrap(rtni_tcgaBRCA) # stopCluster(getOption("cluster")) ## ----eval=FALSE--------------------------------------------------------------- # # Compute the DPI-filtered regulatory network # rtni_tcgaBRCA <- tni.dpi.filter(rtni_tcgaBRCA, eps = 0) ## ----eval=FALSE--------------------------------------------------------------- # # Save the TNI object for subsequent analyses # save(rtni_tcgaBRCA, file="rtni_tcgaBRCA.RData") ## ----eval=FALSE--------------------------------------------------------------- # # For example, to estimate 'alphaB' for 'nB', given 'nA', 'alphaA', and 'betaA' # alphaB <- tni.alpha.adjust(nB = 100, nA = 300, alphaA = 1e-5, betaA = 0.2) # alphaB # # [1] 0.029 ## ----eval=FALSE--------------------------------------------------------------- # library(RTN) # library(Fletcher2013b) # library(pheatmap) ## ----eval=FALSE--------------------------------------------------------------- # # Load 'rtni1st' data object, which includes regulons and expression profiles # data("rtni1st") ## ----eval=FALSE--------------------------------------------------------------- # # A list of transcription factors of interest (here 36 risk-associated TFs) # risk.tfs <- c("AFF3", "AR", "ARNT2", "BRD8", "CBFB", "CEBPB", "E2F2", "E2F3", "ENO1", "ESR1", "FOSL1", "FOXA1", "GATA3", "GATAD2A", "LZTFL1", "MTA2", "MYB", "MZF1", "NFIB", "PPARD", "RARA", "RB1", "RUNX3", "SNAPC2", "SOX10", "SPDEF", "TBX19", "TCEAL1", "TRIM29", "XBP1", "YBX1", "YPEL3", "ZNF24", "ZNF434", "ZNF552", "ZNF587") ## ----eval=FALSE--------------------------------------------------------------- # # Compute regulon activity for individual samples # rtni1st <- tni.gsea2(rtni1st, regulatoryElements = risk.tfs) # metabric_regact <- tni.get(rtni1st, what = "regulonActivity") ## ----eval=FALSE--------------------------------------------------------------- # # Get sample attributes from the 'rtni1st' dataset # metabric_annot <- tni.get(rtni1st, "colAnnotation") # # Get ER+/- and PAM50 attributes for pheatmap # attribs <- c("LumA","LumB","Basal","Her2","Normal","ER+","ER-") # metabric_annot <- metabric_annot[,attribs] ## ----eval=FALSE--------------------------------------------------------------- # # Plot regulon activity profiles # pheatmap(t(metabric_regact$dif), # main="METABRIC cohort 1 (n=977 samples)", # annotation_col = metabric_annot, # show_colnames = FALSE, annotation_legend = FALSE, # clustering_method = "ward.D2", fontsize_row = 6, # clustering_distance_rows = "correlation", # clustering_distance_cols = "correlation") ## ----eval=FALSE, include=FALSE------------------------------------------------ # pheatmap(t(metabric_regact$dif), # main="METABRIC cohort 1 (n=977 samples)", # treeheight_row = 40, treeheight_col = 20, # annotation_col = metabric_annot, # show_colnames = FALSE, annotation_legend = FALSE, # clustering_method = "ward.D2", # clustering_distance_rows = "correlation", # clustering_distance_cols = "correlation", # fontsize_row = 6, height = 5, width = 9, # filename = "fig4.pdf") ## ----eval=FALSE--------------------------------------------------------------- # # Replace samples # rtni1st_tcgasamples <- tni.replace.samples(rtni1st, tcgaBRCA_mRNA_data) ## ----eval=FALSE--------------------------------------------------------------- # # Compute regulon activity for the new samples # rtni1st_tcgasamples <- tni.gsea2(rtni1st_tcgasamples, regulatoryElements = risk.tfs) # tcga_regact <- tni.get(rtni1st_tcgasamples, what = "regulonActivity") ## ----eval=FALSE--------------------------------------------------------------- # # Get sample attributes from the 'rtni1st_tcgasamples' dataset # tcga_annot <- tni.get(rtni1st_tcgasamples, "colAnnotation") # # Adjust PAM50 attributes for pheatmap # tcga_annot <- within(tcga_annot,{ # 'LumA' = ifelse(subtype_BRCA_Subtype_PAM50%in%c("LumA"),1,0) # 'LumB' = ifelse(subtype_BRCA_Subtype_PAM50%in%c("LumB"),1,0) # 'Basal' = ifelse(subtype_BRCA_Subtype_PAM50%in%"Basal",1,0) # 'Her2' = ifelse(subtype_BRCA_Subtype_PAM50%in%c("Her2"),1,0) # 'Normal' = ifelse(subtype_BRCA_Subtype_PAM50%in%c("Normal"),1,0) # }) # attribs <- c("LumA","LumB","Basal","Her2","Normal") # tcga_annot <- tcga_annot[,attribs] ## ----eval=FALSE--------------------------------------------------------------- # # Plot regulon activity profiles # pheatmap(t(tcga_regact$dif), # main="TCGA-BRCA cohort subset (n=500 samples)", # annotation_col = tcga_annot, # show_colnames = FALSE, annotation_legend = FALSE, # clustering_method = "ward.D2", fontsize_row = 6, # clustering_distance_rows = "correlation", # clustering_distance_cols = "correlation") ## ----eval=FALSE, include=FALSE------------------------------------------------ # pheatmap(t(tcga_regact$dif), # main="TCGA-BRCA cohort subset (n=500 samples)", # treeheight_row = 40, treeheight_col = 20, # annotation_col = tcga_annot, # show_colnames = FALSE, annotation_legend = FALSE, # clustering_method = "ward.D2", # clustering_distance_rows = "correlation", # clustering_distance_cols = "correlation", # fontsize_row = 6, height = 4.5, width = 9, # filename = "fig5.pdf") ## ----label='Session information', eval=TRUE, echo=FALSE----------------------- sessionInfo()