## ----setup, include=FALSE------------------------------------------------ knitr::opts_chunk$set(dpi = 300) knitr::opts_chunk$set(cache=FALSE) ## ----, echo = FALSE,hide=TRUE, message=FALSE,warning=FALSE--------------- devtools::load_all(".") ## ----, eval = FALSE------------------------------------------------------ ## source("http://bioconductor.org/biocLite.R") ## useDevel() ## biocLite("TCGAbiolinks") ## ----, eval = FALSE------------------------------------------------------ ## query <- TCGAquery(tumor = c("LGG","GBM"), level = 3, ## platform = c("HumanMethylation450","HumanMethylation27"), ## samples = c("TCGA-19-4065","TCGA-E1-5322-01A-01D-1467-05"), ## version = list(c("HumanMethylation450","LGG",1), ## c("HumanMethylation450","GBM",5))) ## ----, eval = TRUE------------------------------------------------------- query <- TCGAquery(tumor = "gbm") ## ----, eval = TRUE, echo = FALSE----------------------------------------- knitr::kable(disease.table, digits = 2, caption = "List of tumors",row.names = FALSE) ## ----, eval = TRUE------------------------------------------------------- query <- TCGAquery(tumor = "gbm", level = 3) query <- TCGAquery(tumor = "gbm", level = 2) query <- TCGAquery(tumor = "gbm", level = 1) query <- TCGAquery(tumor = "gbm", level = "mage-tab") ## ----, eval = TRUE------------------------------------------------------- query <- TCGAquery(tumor = "gbm", platform = "IlluminaHiSeq_RNASeqV2") ## ----, eval = TRUE, echo = FALSE----------------------------------------- knitr::kable(platform.table[,2:4], digits = 2, caption = "List of tumors",row.names = FALSE) ## ----, eval = TRUE------------------------------------------------------- query <- TCGAquery(tumor = "gbm", center = "mskcc.org") ## ----, eval = TRUE, echo = FALSE----------------------------------------- knitr::kable(center.table, digits = 2, caption = "List of tumors",row.names = FALSE) ## ----, eval = TRUE------------------------------------------------------- # You can define a list of samples to query and download providing relative TCGA barcodes. listSamples <- c("TCGA-E9-A1NG-11A-52R-A14M-07","TCGA-BH-A1FC-11A-32R-A13Q-07", "TCGA-A7-A13G-11A-51R-A13Q-07","TCGA-BH-A0DK-11A-13R-A089-07", "TCGA-E9-A1RH-11A-34R-A169-07","TCGA-BH-A0AU-01A-11R-A12P-07", "TCGA-C8-A1HJ-01A-11R-A13Q-07","TCGA-A7-A13D-01A-13R-A12P-07", "TCGA-A2-A0CV-01A-31R-A115-07","TCGA-AQ-A0Y5-01A-11R-A14M-07") # Query all available platforms with a list of barcode query <- TCGAquery(samples = listSamples) # Query with a partial barcode query <- TCGAquery(samples = "TCGA-61-1743-01A") ## ----, eval = FALSE------------------------------------------------------ ## library(TCGAbiolinks) ## ## BRCA_RNASeqV2_version <- TCGAquery_Version(tumor = "brca", ## platform = "illuminahiseq_rnaseqv2") ## ## ----, eval = TRUE, echo = FALSE----------------------------------------- library(TCGAbiolinks) knitr::opts_chunk$set(comment = "", warning = FALSE, message = FALSE, echo = TRUE, tidy = FALSE, size="small") knitr::kable(BRCA_RNASeqV2_version[,1:4], digits = 2, caption = "Table with date, version and number of samples of BRCA IlluminaHiSeq_RNASeqV2", row.names = FALSE) ## ----, eval = FALSE------------------------------------------------------ ## query <- TCGAquery(tumor = c("LGG","GBM"), platform = c("HumanMethylation450"), level = 3, ## version = list(c("HumanMethylation450","LGG",1))) ## ----, eval = FALSE------------------------------------------------------ ## # Get clinical data ## clinical_brca_data <- TCGAquery_clinic("brca","clinical_patient") ## clinical_uvm_data_bio <- TCGAquery_clinic("uvm","biospecimen_normal_control") ## clinical_brca_data_bio <- TCGAquery_clinic("brca","biospecimen_normal_control") ## clinical_brca_data <- TCGAquery_clinic("brca","clinical_patient") ## ----, eval = FALSE------------------------------------------------------ ## bar <- c("TCGA-G9-6378-02A-11R-1789-07", "TCGA-CH-5767-04A-11R-1789-07", ## "TCGA-G9-6332-60A-11R-1789-07", "TCGA-G9-6336-01A-11R-1789-07", ## "TCGA-G9-6336-11A-11R-1789-07", "TCGA-G9-7336-11A-11R-1789-07", ## "TCGA-G9-7336-04A-11R-1789-07", "TCGA-G9-7336-14A-11R-1789-07", ## "TCGA-G9-7036-04A-11R-1789-07", "TCGA-G9-7036-02A-11R-1789-07", ## "TCGA-G9-7036-11A-11R-1789-07", "TCGA-G9-7036-03A-11R-1789-07", ## "TCGA-G9-7036-10A-11R-1789-07", "TCGA-BH-A1ES-10A-11R-1789-07", ## "TCGA-BH-A1F0-10A-11R-1789-07", "TCGA-BH-A0BZ-02A-11R-1789-07", ## "TCGA-B6-A0WY-04A-11R-1789-07", "TCGA-BH-A1FG-04A-11R-1789-08", ## "TCGA-D8-A1JS-04A-11R-2089-08", "TCGA-AN-A0FN-11A-11R-8789-08", ## "TCGA-AR-A2LQ-12A-11R-8799-08", "TCGA-AR-A2LH-03A-11R-1789-07", ## "TCGA-BH-A1F8-04A-11R-5789-07", "TCGA-AR-A24T-04A-55R-1789-07", ## "TCGA-AO-A0J5-05A-11R-1789-07", "TCGA-BH-A0B4-11A-12R-1789-07", ## "TCGA-B6-A1KN-60A-13R-1789-07", "TCGA-AO-A0J5-01A-11R-1789-07", ## "TCGA-AO-A0J5-01A-11R-1789-07", "TCGA-G9-6336-11A-11R-1789-07", ## "TCGA-G9-6380-11A-11R-1789-07", "TCGA-G9-6380-01A-11R-1789-07", ## "TCGA-G9-6340-01A-11R-1789-07","TCGA-G9-6340-11A-11R-1789-07") ## ## S <- TCGAquery_SampleTypes(bar,"TP") ## S2 <- TCGAquery_SampleTypes(bar,"NB") ## ## # Retrieve multiple tissue types NOT FROM THE SAME PATIENTS ## SS <- TCGAquery_SampleTypes(bar,c("TP","NB")) ## ## # Retrieve multiple tissue types FROM THE SAME PATIENTS ## SSS <- TCGAquery_MatchedCoupledSampleTypes(bar,c("NT","TP")) ## ## # Get clinical data ## clinical_brca_data <- TCGAquery_clinic("brca","clinical_patient") ## female_erpos_herpos <- TCGAquery_clinicFilt(bar,clinical_brca_data, ## HER="Positive", ## gender="FEMALE", ## ER="Positive") ## ----, eval = TRUE, echo = FALSE----------------------------------------- bar <- c("TCGA-G9-6378-02A-11R-1789-07", "TCGA-CH-5767-04A-11R-1789-07", "TCGA-G9-6332-60A-11R-1789-07", "TCGA-G9-6336-01A-11R-1789-07", "TCGA-G9-6336-11A-11R-1789-07", "TCGA-G9-7336-11A-11R-1789-07", "TCGA-G9-7336-04A-11R-1789-07", "TCGA-G9-7336-14A-11R-1789-07", "TCGA-G9-7036-04A-11R-1789-07", "TCGA-G9-7036-02A-11R-1789-07", "TCGA-G9-7036-11A-11R-1789-07", "TCGA-G9-7036-03A-11R-1789-07", "TCGA-G9-7036-10A-11R-1789-07", "TCGA-BH-A1ES-10A-11R-1789-07", "TCGA-BH-A1F0-10A-11R-1789-07", "TCGA-BH-A0BZ-02A-11R-1789-07", "TCGA-D8-A1JS-04A-11R-2089-08", "TCGA-AN-A0FN-11A-11R-8789-08", "TCGA-AR-A2LQ-12A-11R-8799-08", "TCGA-AR-A2LH-03A-11R-1789-07", "TCGA-BH-A1F8-04A-11R-5789-07", "TCGA-AR-A24T-04A-55R-1789-07", "TCGA-AO-A0J5-05A-11R-1789-07", "TCGA-BH-A0B4-11A-12R-1789-07", "TCGA-B6-A1KN-60A-13R-1789-07", "TCGA-AO-A0J5-01A-11R-1789-07", "TCGA-AO-A0J5-01A-11R-1789-07", "TCGA-G9-6336-11A-11R-1789-07", "TCGA-G9-6380-11A-11R-1789-07", "TCGA-G9-6380-01A-11R-1789-07", "TCGA-G9-6340-01A-11R-1789-07","TCGA-G9-6340-11A-11R-1789-07") female_erpos_herpos <- TCGAquery_clinicFilt(bar, clinBRCA, HER = "Positive", gender = "FEMALE", ER = "Positive") print(female_erpos_herpos) ## ----, eval = FALSE------------------------------------------------------ ## # Check with subtypes from TCGAprepare and update examples ## GBM_path_subtypes <- TCGAquery_subtype(tumor = "gbm") ## ## LGG_path_subtypes <- TCGAquery_subtype(tumor = "lgg") ## ## LGG_clinic <- TCGAquery_clinic(tumor = "LGG", ## clinical_data_type = "clinical_patient") ## ----, eval = TRUE, echo = FALSE----------------------------------------- knitr::kable(lgg.gbm.subtype[1:10,c(1,2,3,4)], digits = 2, caption = "Table common samples among platforms from TCGAquery", row.names = TRUE) ## ----, eval = FALSE,echo=TRUE,message=FALSE,warning=FALSE---------------- ## query <- TCGAquery(tumor = "brca",level = 3) ## matSamples <- TCGAquery_integrate(query) ## matSamples[c(1,4,9),c(1,4,9)] ## ----, eval = TRUE, echo = FALSE----------------------------------------- query <- TCGAquery(tumor = "brca",level = 3) matSamples <- TCGAquery_integrate(query) knitr::kable(matSamples[c(1,4,9),c(1,4,9)], digits = 2, caption = "Table common samples among platforms from TCGAquery", row.names = TRUE) ## ----, eval = FALSE------------------------------------------------------ ## # First perform DEGs with TCGAanalyze ## # See previous section ## library(TCGAbiolinks) ## ## # Select only transcription factors (TFs) from DEGs ## TFs <- EAGenes[EAGenes$Family =="transcription regulator",] ## TFs_inDEGs <- intersect(TFs$Gene, dataDEGsFiltLevel$mRNA ) ## dataDEGsFiltLevelTFs <- dataDEGsFiltLevel[TFs_inDEGs,] ## ## # Order table DEGs TFs according to Delta decrease ## dataDEGsFiltLevelTFs <- dataDEGsFiltLevelTFs[order(dataDEGsFiltLevelTFs$Delta,decreasing = TRUE),] ## ## # Find Pubmed of TF studied related to cancer ## tabDEGsTFPubmed <- TCGAquery_investigate("breast", dataDEGsFiltLevelTFs, topgenes = 10) ## ## ----, eval = TRUE, echo = FALSE----------------------------------------- library(TCGAbiolinks) library(stringr) tabDEGsTFPubmed$PMID <- str_sub(tabDEGsTFPubmed$PMID,0,30) knitr::kable(tabDEGsTFPubmed, digits = 2, caption = "Table with most studied TF in pubmed related to a specific cancer",row.names = FALSE) #LatexPrintTableforPresentation(Table = tabDEGsTFPubmed,rowsForPage = nrow(tabDEGsTFPubmed), TableTitle = "tabDEGsTFPubmed", LabelTitle = "tabDEGsTFPubmed", withrows = T) ## ----, eval = FALSE------------------------------------------------------ ## library(TCGAbiolinks) ## ## # Define a list of package to find number of downloads ## listPackage <-c("limma","edgeR","survcomp") ## ## tabPackage <- TCGAquery_Social(siteToFind ="bioconductor.org",listPackage) ## ## ## # define a keyword to find in support.bioconductor.org returing a table with suggested packages ## tabPackageKey <- TCGAquery_Social(siteToFind ="support.bioconductor.org" ,KeyInfo = "tcga") ## ----, eval = TRUE, echo = FALSE----------------------------------------- library(TCGAbiolinks) listPackage <- c("limma","edgeR","survcomp") tabPackage <- TCGAquery_Social(siteToFind = "bioconductor.org", listPackage) knitr::kable(head(tabPackage), digits = 2, caption = "Table with number of downloads about a list of packages",row.names = FALSE) knitr::kable(tabPackage2[c(1,4,5,7),], digits = 2, caption = "Find most related question in support.bioconductor.org with keyword = tcga",row.names = FALSE) ## ----, eval = FALSE------------------------------------------------------ ## library(TCGAbiolinks) ## ## # Find most related question in biostar with TCGA ## tabPackage1 <- TCGAquery_Social(siteToFind ="biostars.org",KeyInfo = "TCGA") ## ## # Find most related question in biostar with package ## tabPackage2 <- TCGAquery_Social(siteToFind ="biostars.org",KeyInfo = "package") ## ## ----, eval = TRUE, echo = FALSE----------------------------------------- library(TCGAbiolinks) # Find most related question in biostar with TCGA #tabPackage <- TCGAquery_Social(siteToFind ="biostars.org",KeyInfo = "package") knitr::kable(head(tabPackage1), digits = 2, caption = "Find most related question in biostar with TCGA",row.names = FALSE) knitr::kable(tabPackage2[c(1,5,7),], digits = 2, caption = "Find most related question in biostar with package",row.names = FALSE) ## ----, eval = FALSE------------------------------------------------------ ## # get all samples from the query and save them in the TCGA folder ## # samples from IlluminaHiSeq_RNASeqV2 with type rsem.genes.results ## # samples to normalize later ## TCGAdownload(query, path = "data", type = "rsem.genes.results") ## ## TCGAdownload(query, path = "data", type = "rsem.isoforms.normalized_results") ## ## TCGAdownload(query, path = "dataBrca", type = "rsem.genes.results", ## samples = c("TCGA-E9-A1NG-11A-52R-A14M-07", ## "TCGA-BH-A1FC-11A-32R-A13Q-07")) ## ----, eval = FALSE------------------------------------------------------ ## # get all samples from the query and save them in the TCGA folder ## # samples from IlluminaHiSeq_RNASeqV2 with type rsem.genes.results ## # samples to normalize later ## data <- TCGAprepare(query, dir = "data", save = TRUE, filename = "myfile.rda") ## ----,eval=FALSE,echo=FALSE,message=FALSE,warning=FALSE------------------ ## library(TCGAbiolinks) ## library(SummarizedExperiment) ## query <- TCGAquery(tumor = "READ", platform = "IlluminaHiSeq_RNASeqV2", level = 3, ## samples = c("TCGA-DY-A1DE-01A-11R-A155-07", ## "TCGA-DY-A0XA-01A-11R-A155-07")) ## TCGAdownload(query,samples = c("TCGA-DY-A1DE-01A-11R-A155-07", ## "TCGA-DY-A0XA-01A-11R-A155-07"), ## type = "rsem.genes.normalized_results",path = "../dataREAD") ## dataREAD <- TCGAprepare(query,dir = "../dataREAD", type = "rsem.genes.normalized_results", ## samples = c("TCGA-DY-A1DE-01A-11R-A155-07", ## "TCGA-DY-A0XA-01A-11R-A155-07")) ## ----, eval = TRUE,message=FALSE,warning=FALSE--------------------------- library(TCGAbiolinks) library(SummarizedExperiment) head(assay(dataREAD,"normalized_count")) ## ----,eval=FALSE,echo=FALSE,message=FALSE,warning=FALSE------------------ ## library(SummarizedExperiment) ## library(TCGAbiolinks) ## query <- TCGAquery(tumor = "READ", platform = "IlluminaHiSeq_RNASeqV2", level = 3, ## samples = c("TCGA-DY-A1DE-01A-11R-A155-07", ## "TCGA-DY-A0XA-01A-11R-A155-07")) ## TCGAdownload(query,samples = c("TCGA-DY-A1DE-01A-11R-A155-07", ## "TCGA-DY-A0XA-01A-11R-A155-07"), ## type = "rsem.genes.normalized_results",path = "../dataREAD") ## ## dataREAD_df <- TCGAprepare(query,dir = "../dataREAD", ## type = "rsem.genes.normalized_results", ## samples = c("TCGA-DY-A1DE-01A-11R-A155-07", ## "TCGA-DY-A0XA-01A-11R-A155-07"), ## summarizedExperiment = FALSE) ## ----, eval = TRUE------------------------------------------------------- library(TCGAbiolinks) class(dataREAD_df) dim(dataREAD_df) head(dataREAD_df) ## ----, eval = FALSE------------------------------------------------------ ## ## # Define a list of samples to query and download providing relative TCGA barcodes. ## samplesList <- c("TCGA-02-0046-10A-01D-0182-01", ## "TCGA-02-0052-01A-01D-0182-01", ## "TCGA-02-0033-10A-01D-0182-01", ## "TCGA-02-0034-01A-01D-0182-01", ## "TCGA-02-0007-01A-01D-0182-01") ## ## # Query platform Genome_Wide_SNP_6 with a list of barcode ## query <- TCGAquery(tumor = "gbm", level = 3, platform = "Genome_Wide_SNP_6") ## ## # Download a list of barcodes with platform Genome_Wide_SNP_6 ## TCGAdownload(query, path = "samples") ## ## # Prepare matrix ## GBM_CNV <- TCGAprepare(query, dir = "samples", type = ".hg19.seg.txt") ## ----, eval = FALSE------------------------------------------------------ ## ############## Get tumor samples with TCGAbiolinks ## library(TCGAbiolinks) ## path <- "kirc" ## query <- TCGAquery(tumor = "KIRC", level = 3, platform = "HumanMethylation450") ## TCGAdownload(query, path =path) ## ## kirc.met <- TCGAprepare(query,dir = path, ## save = TRUE, ## filename = "metKirc.rda", ## summarizedExperiment = FALSE) ## ## kirc.met <- TCGAprepare_elmer(kirc.met, ## platform = "HumanMethylation450", ## save = TRUE, ## met.na.cut = 0.2) ## ## # Step 1.2 download expression data ## query.rna <- TCGAquery(tumor="KIRC",level=3, platform="IlluminaHiSeq_RNASeqV2") ## TCGAdownload(query.rna,path=path,type = "rsem.genes.normalized_results") ## ## kirc.exp <- TCGAprepare(query.rna, dir=path, save = TRUE, ## type = "rsem.genes.normalized_results", ## filename = "expKirc.rda", summarizedExperiment = FALSE) ## ## kirc.exp <- TCGAprepare_elmer(kirc.exp, ## save = TRUE, ## platform = "IlluminaHiSeq_RNASeqV2") ## ## # Step 2 prepare mee object ## library(ELMER) ## library(parallel) ## ## geneAnnot <- txs() ## geneAnnot$GENEID <- paste0("ID",geneAnnot$GENEID) ## geneInfo <- promoters(geneAnnot,upstream = 0, downstream = 0) ## probe <- get.feature.probe() ## mee <- fetch.mee(meth = kirc.met, exp = kirc.exp, TCGA = TRUE, ## probeInfo = probe, geneInfo = geneInfo) ## save(mee,file="case4mee.rda") ## ----, eval = FALSE------------------------------------------------------ ## ## # You can define a list of samples to query and download providing relative TCGA barcodes. ## ## listSamples <- c("TCGA-E9-A1NG-11A-52R-A14M-07","TCGA-BH-A1FC-11A-32R-A13Q-07", ## "TCGA-A7-A13G-11A-51R-A13Q-07","TCGA-BH-A0DK-11A-13R-A089-07", ## "TCGA-E9-A1RH-11A-34R-A169-07","TCGA-BH-A0AU-01A-11R-A12P-07", ## "TCGA-C8-A1HJ-01A-11R-A13Q-07","TCGA-A7-A13D-01A-13R-A12P-07", ## "TCGA-A2-A0CV-01A-31R-A115-07","TCGA-AQ-A0Y5-01A-11R-A14M-07") ## ## # Query platform IlluminaHiSeq_RNASeqV2 with a list of barcode ## query <- TCGAquery(tumor = "brca", samples = listSamples, ## platform = "IlluminaHiSeq_RNASeqV2", level = "3") ## ## # Download a list of barcodes with platform IlluminaHiSeq_RNASeqV2 ## TCGAdownload(query, path = "../dataBrca", type = "rsem.genes.results",samples = listSamples) ## ## # Prepare expression matrix with gene id in rows and samples (barcode) in columns ## # rsem.genes.results as values ## BRCARnaseq_assay <- TCGAprepare(query,"../dataBrca",type = "rsem.genes.results") ## ## BRCAMatrix <- assay(BRCARnaseq_assay,"raw_counts") ## ## # For gene expression if you need to see a boxplot correlation and AAIC plot ## # to define outliers you can run ## BRCARnaseq_CorOutliers <- TCGAanalyze_Preprocessing(BRCARnaseq_assay) ## ----, eval = TRUE, echo = FALSE,size = 8-------------------------------- library(TCGAbiolinks) dataGE <- dataBRCA[sample(rownames(dataBRCA),10),sample(colnames(dataBRCA),7)] knitr::kable(dataGE[1:10,2:3], digits = 2, caption = "Example of a matrix of gene expression (10 genes in rows and 2 samples in columns)", row.names = TRUE) ## ----, fig.width=6, fig.height=4, echo=FALSE, fig.align="center"--------- library(png) library(grid) img <- readPNG("PreprocessingOutput.png") grid.raster(img) ## ----, eval = TRUE------------------------------------------------------- # Downstream analysis using gene expression data # TCGA samples from IlluminaHiSeq_RNASeqV2 with type rsem.genes.results # save(dataBRCA, geneInfo , file = "dataGeneExpression.rda") library(TCGAbiolinks) # normalization of genes dataNorm <- TCGAanalyze_Normalization(tabDF = dataBRCA, geneInfo = geneInfo) # quantile filter of genes dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm, method = "quantile", qnt.cut = 0.25) # selection of normal samples "NT" samplesNT <- TCGAquery_SampleTypes(barcode = colnames(dataFilt), typesample = c("NT")) # selection of tumor samples "TP" samplesTP <- TCGAquery_SampleTypes(barcode = colnames(dataFilt), typesample = c("TP")) # Diff.expr.analysis (DEA) dataDEGs <- TCGAanalyze_DEA(mat1 = dataFilt[,samplesNT], mat2 = dataFilt[,samplesTP], Cond1type = "Normal", Cond2type = "Tumor", fdr.cut = 0.01 , logFC.cut = 1, method = "glmLRT") # DEGs table with expression values in normal and tumor samples dataDEGsFiltLevel <- TCGAanalyze_LevelTab(dataDEGs,"Tumor","Normal", dataFilt[,samplesTP],dataFilt[,samplesNT]) ## ----, eval = TRUE, echo = FALSE----------------------------------------- library(TCGAbiolinks) dataDEGsFiltLevel$FDR <- format(dataDEGsFiltLevel$FDR, scientific = TRUE) knitr::kable(dataDEGsFiltLevel[1:10,], digits = 2, caption = "Table DEGs after DEA",row.names = FALSE) ## ----, eval = FALSE------------------------------------------------------ ## library(TCGAbiolinks) ## # Enrichment Analysis EA ## # Gene Ontology (GO) and Pathway enrichment by DEGs list ## Genelist <- rownames(dataDEGsFiltLevel) ## ## system.time(ansEA <- TCGAanalyze_EAcomplete(TFname="DEA genes Normal Vs Tumor",Genelist)) ## ## # Enrichment Analysis EA (TCGAVisualize) ## # Gene Ontology (GO) and Pathway enrichment barPlot ## ## TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP), ## GOBPTab = ansEA$ResBP, ## GOCCTab = ansEA$ResCC, ## GOMFTab = ansEA$ResMF, ## PathTab = ansEA$ResPat, ## nRGTab = Genelist, ## nBar = 10) ## ## ----, fig.width=6, fig.height=4, echo=FALSE, fig.align="center"--------- library(png) library(grid) img <- readPNG("EAplot.png") grid.raster(img) ## ----, eval = FALSE------------------------------------------------------ ## clin.gbm <- TCGAquery_clinic("gbm", "clinical_patient") ## clin.lgg <- TCGAquery_clinic("lgg", "clinical_patient") ## ## TCGAanalyze_survival(plyr::rbind.fill(clin.lgg,clin.gbm), ## "radiation_therapy", ## main = "TCGA Set\nLGG and GBM",height = 10, width=10) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("survival.png") grid.raster(img) ## ----, eval = FALSE------------------------------------------------------ ## library(TCGAbiolinks) ## # Survival Analysis SA ## ## clinical_patient_Cancer <- TCGAquery_clinic("brca","clinical_patient") ## dataBRCAcomplete <- log2(BRCA_rnaseqv2) ## ## tokenStop<- 1 ## ## tabSurvKMcomplete <- NULL ## ## for( i in 1: round(nrow(dataBRCAcomplete)/100)){ ## message( paste( i, "of ", round(nrow(dataBRCAcomplete)/100))) ## tokenStart <- tokenStop ## tokenStop <-100*i ## tabSurvKM<-TCGAanalyze_SurvivalKM(clinical_patient_Cancer, ## dataBRCAcomplete, ## Genelist = rownames(dataBRCAcomplete)[tokenStart:tokenStop], ## Survresult = F, ## ThreshTop=0.67, ## ThreshDown=0.33) ## ## tabSurvKMcomplete <- rbind(tabSurvKMcomplete,tabSurvKM) ## } ## ## tabSurvKMcomplete <- tabSurvKMcomplete[tabSurvKMcomplete$pvalue < 0.01,] ## tabSurvKMcomplete <- tabSurvKMcomplete[!duplicated(tabSurvKMcomplete$mRNA),] ## rownames(tabSurvKMcomplete) <-tabSurvKMcomplete$mRNA ## tabSurvKMcomplete <- tabSurvKMcomplete[,-1] ## tabSurvKMcomplete <- tabSurvKMcomplete[order(tabSurvKMcomplete$pvalue, decreasing=F),] ## ## tabSurvKMcompleteDEGs <- tabSurvKMcomplete[ ## rownames(tabSurvKMcomplete) %in% dataDEGsFiltLevel$mRNA, ## ] ## ----, fig.width=6, fig.height=4, echo=FALSE, fig.align="center"--------- tabSurvKMcompleteDEGs$pvalue <- format(tabSurvKMcompleteDEGs$pvalue, scientific = TRUE) knitr::kable(tabSurvKMcompleteDEGs[1:10,1:4], digits = 2, caption = "Table KM-survival genes after SA", row.names = TRUE) knitr::kable(tabSurvKMcompleteDEGs[1:10,5:7], digits = 2, row.names = TRUE) ## ----, eval = FALSE------------------------------------------------------ ## data <- TCGAanalyze_DMR(data, groupCol = "cluster.meth",subgroupCol = "disease", ## group.legend = "Groups", subgroup.legend = "Tumor", ## print.pvalue = TRUE) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("figure5met.png") grid.raster(img) ## ----, eval = FALSE------------------------------------------------------ ## # normalization of genes ## dataNorm <- TCGAbiolinks::TCGAanalyze_Normalization(dataBRCA, geneInfo) ## ## # quantile filter of genes ## dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm, ## method = "quantile", ## qnt.cut = 0.25) ## ## # Principal Component Analysis plot for ntop selected DEGs ## TCGAvisualize_PCA(dataFilt,dataDEGsFiltLevel, ntopgenes = 200) ## ----, fig.width=6, fig.height=4, echo=FALSE, fig.align="center"--------- library(png) library(grid) img <- readPNG("PCAtop200DEGs.png") grid.raster(img) ## ----, eval = FALSE------------------------------------------------------ ## library(TCGAbiolinks) ## # Survival Analysis SA ## ## clinical_patient_Cancer <- TCGAquery_clinic("brca","clinical_patient") ## dataBRCAcomplete <- log2(BRCA_rnaseqv2) ## ## tokenStop<- 1 ## ## tabSurvKMcomplete <- NULL ## ## for( i in 1: round(nrow(dataBRCAcomplete)/100)){ ## message( paste( i, "of ", round(nrow(dataBRCAcomplete)/100))) ## tokenStart <- tokenStop ## tokenStop <-100*i ## tabSurvKM<-TCGAanalyze_SurvivalKM(clinical_patient_Cancer, ## dataBRCAcomplete, ## Genelist = rownames(dataBRCAcomplete)[tokenStart:tokenStop], ## Survresult = F,ThreshTop=0.67,ThreshDown=0.33) ## ## tabSurvKMcomplete <- rbind(tabSurvKMcomplete,tabSurvKM) ## } ## ## tabSurvKMcomplete <- tabSurvKMcomplete[tabSurvKMcomplete$pvalue < 0.01,] ## tabSurvKMcomplete <- tabSurvKMcomplete[!duplicated(tabSurvKMcomplete$mRNA),] ## rownames(tabSurvKMcomplete) <-tabSurvKMcomplete$mRNA ## tabSurvKMcomplete <- tabSurvKMcomplete[,-1] ## tabSurvKMcomplete <- tabSurvKMcomplete[order(tabSurvKMcomplete$pvalue, decreasing=F),] ## ## tabSurvKMcompleteDEGs <- tabSurvKMcomplete[rownames(tabSurvKMcomplete) %in% dataDEGsFiltLevel$mRNA,] ## ## tflist <- EAGenes[EAGenes$Family == "transcription regulator","Gene"] ## tabSurvKMcomplete_onlyTF <- tabSurvKMcomplete[rownames(tabSurvKMcomplete) %in% tflist,] ## ## TabCoxNet <- TCGAvisualize_SurvivalCoxNET(clinical_patient_Cancer,dataBRCAcomplete, ## Genelist = rownames(tabSurvKMcompleteDEGs), ## scoreConfidence = 700,titlePlot = "TCGAvisualize_SurvivalCoxNET Example") ## ----, fig.width=6, fig.height=4, echo=FALSE, fig.align="center"--------- library(png) library(grid) img <- readPNG("SurvivalCoxNETOutput.png") grid.raster(img) ## ----, eval = FALSE------------------------------------------------------ ## TCGAvisualize_meanMethylation(data,"group") ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("meanmet.png") grid.raster(img) ## ----, eval = FALSE------------------------------------------------------ ## starburst <- TCGAvisualize_starburst(coad.SummarizeExperiment, ## different.experssion.analysis.data, ## group1 = "CIMP.H", ## group2 = "CIMP.L", ## met.p.cut = 10^-5, ## exp.p.cut=10^-5, ## names = TRUE) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("figure5star.png") grid.raster(img) ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## PlatformCancer <- "IlluminaHiSeq_RNASeqV2" ## dataType <- "rsem.genes.results" ## pathGBM<- "../dataGBM" ## pathLGG <- "../dataLGG" ## ## library(BiocInstaller) ## useDevel() # we need Devel for SummarizedExperiment package ## library(SummarizedExperiment) ## library(TCGAbiolinks) ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## library(TCGAbiolinks) ## ## # defining common parameters ## cancer <- "BRCA" ## PlatformCancer <- "IlluminaHiSeq_RNASeqV2" ## dataType <- "rsem.genes.results" ## pathCancer <- paste0("../data",cancer) ## ## datQuery <- TCGAquery(tumor = cancer, platform = PlatformCancer, level = "3") ## lsSample <- TCGAquery_samplesfilter(query = datQuery) ## ## # get subtype information ## dataSubt <- TCGAquery_subtype(tumor = cancer) ## ## # Which samples are Primary Solid Tumor ## dataSmTP <- TCGAquery_SampleTypes(barcode = lsSample$IlluminaHiSeq_RNASeqV2, ## typesample = "TP") ## # Which samples are Solid Tissue Normal ## dataSmTN <- TCGAquery_SampleTypes(barcode = lsSample$IlluminaHiSeq_RNASeqV2, ## typesample ="NT") ## # get clinical data ## dataClin <- TCGAquery_clinic(tumor = cancer, ## clinical_data_type = "clinical_patient") ## ## TCGAdownload(data = datQuery, ## path = pathCancer, ## type = dataType, ## samples =c(dataSmTP,dataSmTN)) ## ## dataAssy <- TCGAprepare(query = datQuery, ## dir = pathCancer, ## type = dataType, ## save = TRUE, ## summarizedExperiment = TRUE, ## samples = c(dataSmTP,dataSmTN), ## filename = paste0(cancer,"_",PlatformCancer,".rda")) ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## ## dataPrep <- TCGAanalyze_Preprocessing(object = dataAssy, ## cor.cut = 0.6) ## ## dataNorm <- TCGAanalyze_Normalization(tabDF = dataPrep, ## geneInfo = geneInfo, ## method = "gcContent") ## ## dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm, ## method = "quantile", ## qnt.cut = 0.25) ## ## dataDEGs <- TCGAanalyze_DEA(mat1 = dataFilt[,dataSmTN], ## mat2 = dataFilt[,dataSmTP], ## Cond1type = "Normal", ## Cond2type = "Tumor", ## fdr.cut = 0.01 , ## logFC.cut = 1, ## method = "glmLRT") ## ## ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## ansEA <- TCGAanalyze_EAcomplete(TFname="DEA genes Normal Vs Tumor", ## RegulonList = rownames(dataDEGs)) ## ## TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP), ## GOBPTab = ansEA$ResBP, ## GOCCTab = ansEA$ResCC, ## GOMFTab = ansEA$ResMF, ## PathTab = ansEA$ResPat, ## nRGTab = rownames(dataDEGs), ## nBar = 20) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case1_EA.png") grid.raster(img) ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## dataSurv <- TCGAanalyze_SurvivalKM(clinical_patient = dataClin, ## dataGE = dataFilt, ## Genelist = rownames(dataDEGs), ## Survresult = FALSE, ## ThreshTop = 0.67, ## ThreshDown = 0.33, ## p.cut = 0.05) ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## ## require(dnet) # to change ## org.Hs.string <- dRDataLoader(RData = "org.Hs.string") ## ## TabCoxNet <- TCGAvisualize_SurvivalCoxNET(dataClin, ## dataFilt, ## Genelist = rownames(dataSurv), ## scoreConfidence = 700, ## org.Hs.string = org.Hs.string, ## titlePlot = "Case Study n.1 dnet") ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case1_dnet.png") grid.raster(img) ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## library(TCGAbiolinks) ## cancer <- "LGG" ## PlatformCancer <- "IlluminaHiSeq_RNASeqV2" ## dataType <- "rsem.genes.results" ## pathCancer <- paste0("../data",cancer) ## ## # Result....Function.....parameters...p1...pn..................................time execution ## ## datQuery <- TCGAquery(tumor = cancer, platform = PlatformCancer, level = "3") # time = 0.093s ## lsSample <- TCGAquery_samplesfilter(query = datQuery) ## dataSubt <- TCGAquery_subtype(tumor = cancer) ## dataSmTP <- TCGAquery_SampleTypes(barcode = lsSample$IlluminaHiSeq_RNASeqV2, ## typesample = "TP") ## dataClin <- TCGAquery_clinic(tumor = cancer, ## clinical_data_type = "clinical_patient") ## ## TCGAdownload(data = datQuery, path = pathCancer, type = dataType, ## samples = dataSmTP ) ## ## dataAssy <- TCGAprepare(query = datQuery, ## dir = pathCancer, ## type = dataType, ## save = TRUE, ## summarizedExperiment = TRUE, ## samples = dataSmTP, ## filename = paste0(cancer,"_",PlatformCancer,".rda")) ## ## dataPrep <- TCGAanalyze_Preprocessing(object = dataAssy,cor.cut = 0.6) # time = 13.028s ## dataNorm <- TCGAanalyze_Normalization(tabDF = dataPrep, ## geneInfo = geneInfo, ## method = "gcContent") # time = 165.577s ## ## datFilt1 <- TCGAanalyze_Filtering(tabDF = dataNorm,method = "varFilter") ## datFilt2 <- TCGAanalyze_Filtering(tabDF = datFilt1,method = "filter1") ## datFilt <- TCGAanalyze_Filtering(tabDF = datFilt2,method = "filter2") ## ## data_Hc1 <- TCGAanalyze_Clustering(tabDF = datFilt,method = "hclust", methodHC = "ward.D2") ## data_Hc2 <- TCGAanalyze_Clustering(tabDF = datFilt, ## method = "consensus", ## methodHC = "ward.D2") # time = 207.389 ## # deciding number of tree to cuts ## cut.tree <-4 ## paste0(c("EC"),(1:cut.tree)) ## ## ## consensusClusters contains barcodes for 4 groups ## ans <- hclust(ddist <- dist(datFilt), method = "ward.D2") ## hhc <- data_Hc2[[cut.tree]]$consensusTree ## consensusClusters<-data_Hc2[[cut.tree]]$consensusClass ## sampleOrder <- consensusClusters[hhc$order] ## ## consensusClusters <- as.factor(data_Hc2[[cut.tree]]$clrs[[1]]) ## names(consensusClusters) <- attr(ddist, "Labels") ## names(consensusClusters) <- substr(names(consensusClusters),1,12) ## ## # adding information about gropus from consensus Cluster in clinical data ## dataClin <- cbind(dataClin, groupsHC = matrix(0,nrow(dataClin),1)) ## rownames(dataClin) <- dataClin$bcr_patient_barcode ## ## for( i in 1: nrow(dataClin)){ ## currSmp <- dataClin$bcr_patient_barcode[i] ## dataClin[currSmp,"groupsHC"] <- as.character(consensusClusters[currSmp]) ## } ## ## # plotting survival for groups EC1, EC2, EC3, EC4 ## ## TCGAanalyze_survival(data = dataClin, ## clusterCol = "groupsHC", ## main = "TCGA kaplan meier survival plot from consensus cluster", ## height = 10, ## width=10, ## legend = "RNA Group", ## labels=paste0(c("EC"),(1:cut.tree)), ## color = as.character(levels(consensusClusters)), ## filename = "case2_surv.png") ## dev.off() ## TCGAvisualize_BarPlot(DFfilt = datFilt, ## DFclin = dataClin, ## DFsubt = dataSubt, ## data_Hc2 = data_Hc2, ## Subtype = "IDH.1p19q.Subtype", ## cbPalette = c("cyan","tomato","gold"), ## filename = "case2_Idh.png", ## height = 10, ## width=10, ## dpi =300) ## ## TCGAvisualize_BarPlot(DFfilt = datFilt, ## DFclin = dataClin, ## DFsubt = dataSubt, ## data_Hc2 = data_Hc2, ## Subtype = "MethylationCluster", ## cbPalette = c("black","orchid3","palegreen4","sienna3", "steelblue4"), ## filename = "case2_Met.png", ## height = 10, ## width=10, ## dpi =300) ## ## TCGAvisualize_BarPlot(DFfilt = datFilt, ## DFclin = dataClin, ## DFsubt = dataSubt, ## data_Hc2 = data_Hc2, ## Subtype = "AGE", ## cbPalette = c("yellow2","yellow3","yellow4"), ## filename = "case2_Age.png", ## height = 10, ## width=10, ## dpi =300) ## ## dev.off() ## pdf(file="case2_Heatmap2.pdf") ## TCGAvisualize_Heatmap(DFfilt = datFilt, ## DFclin = dataClin, ## DFsubt = dataSubt, ## data_Hc2= data_Hc2) ## dev.off() ## ## ## # Convert images from pdf to png. ## library(animation) ## ani.options(outdir = getwd()) ## ## im.convert("case2_Heatmap2.pdf", ## output = "case2_Heatmap2.png", ## extra.opts="-density 150") ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case2_surv.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case2_Idh.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case2_Met.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case2_Met.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case2_Age.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case2_Heatmap.png") grid.raster(img) ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## #----------------------------------- ## # STEP 1: Search, download, prepare | ## #----------------------------------- ## # 1.1 - Methylation ## # ---------------------------------- ## query.met <- TCGAquery(tumor = c("coad"), ## platform = c("HumanMethylation27", ## "HumanMethylation450"), ## level = 3) ## ## TCGAdownload(query.met, path = "/dados/ibm/comparing/biolinks/coad/") ## ## coad.met <- TCGAprepare(query = query.met, ## dir = "/dados/ibm/comparing/biolinks/coad/", ## save = TRUE, ## add.subtype= TRUE, ## filename = "metcoad.rda", ## reannotate = TRUE) ## ## #----------------------------------- ## # 1.2 - Expression ## # ---------------------------------- ## ## coad.query.exp <- TCGAquery(tumor = "coad", ## platform = "IlluminaGA_RNASeqV2", ## level = 3) ## ## TCGAdownload(coad.query.exp, ## path = "/dados/ibm/comparing/biolinks/RNA/", ## type = "rsem.genes.results") ## ## coad.exp <- TCGAprepare(query = coad.query.exp, ## dir = "/dados/ibm/comparing/biolinks/RNA/", ## type = "rsem.genes.results", ## add.subtype= TRUE, ## save = T, filename = "coadexp.rda") ## ## ## # removing the samples without classification ## coad.met <- subset(coad.met,select = !(colData(coad.met)$methylation_subtype %in% c(NA))) ## ## #-------------------------------------------- ## # STEP 2: Analysis ## #-------------------------------------------- ## # 2.1 - Mean methylation of samples ## # ------------------------------------------- ## TCGAvisualize_meanMethylation(coad.met, ## groupCol = "methylation_subtype", ## subgroupCol = "hypermutated", ## group.legend = "Groups", ## subgroup.legend = "hypermutated", ## filename = "coad_mean.png") ## ## #------------------------------------------------- ## # 2.2 - DMR - Methylation analysis - volcano plot ## # ------------------------------------------------ ## coad.aux <- subset(coad.met, ## select = colData(coad.met)$methylation_subtype %in% c("CIMP.L","CIMP.H")) ## ## ## # na.omit ## coad.aux <- subset(coad.aux,subset = (rowSums(is.na(assay(coad.aux))) == 0)) ## ## # Volcano plot ## coad.aux <- TCGAanalyze_DMR(coad.aux, groupCol = "methylation_subtype", ## group1 = "CIMP.H", ## group2="CIMP.L", ## p.cut = 10^-5, ## diffmean.cut = 0.25, ## legend = "State", ## plot.filename = "coad_CIMPHvsCIMPL_metvolcano.png") ## ## save(coad.aux,file = "coad_pvalue.rda") ## ## #------------------------------------------------- ## # 2.3 - DEA - Expression analysis - volcano plot ## # ------------------------------------------------ ## ## coad.exp.aux <- subset(coad.exp, ## select = colData(coad.exp)$methylation_subtype %in% c("CIMP.H","CIMP.L")) ## ## idx <- colData(coad.exp.aux)$methylation_subtype %in% c("CIMP.H") ## idx2 <- colData(coad.exp.aux)$methylation_subtype %in% c("CIMP.L") ## ## dataPrep <- TCGAanalyze_Preprocessing(object = coad.exp.aux, ## cor.cut = 0.6) ## ## dataNorm <- TCGAanalyze_Normalization(tabDF = dataPrep, ## geneInfo = geneInfo, ## method = "gcContent") ## ## dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm, ## qnt.cut = 0.25, ## method='quantile') ## ## dataDEGs <- TCGAanalyze_DEA(mat1 = dataFilt[,idx], ## mat2 = dataFilt[,idx2], ## Cond1type = "CIMP.H", ## Cond2type = "CIMP.l", ## method = "glmLRT") ## ## TCGAVisualize_volcano(dataDEGs$logFC,dataDEGs$FDR, ## filename = "Case3_volcanoexp.png", ## x.cut = 3, ## y.cut = 10^-4, ## names = rownames(dataDEGs), ## xlab = " Gene expression fold change (Log2)", ## legend = "State", ## title = "Volcano plot (CIMP.l vs CIMP.H)") ## ## #------------------------------------------ ## # 2.4 - Starburst plot ## # ----------------------------------------- ## starburst <- TCGAvisualize_starburst(coad.aux, dataDEGs,"CIMP.H","CIMP.L", ## filename = "starburst.png", ## met.p.cut = 10^-5, ## exp.p.cut = 10^-4, ## diffmean.cut = 0.25, ## logFC.cut = 3, ## names = TRUE) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("figure5exp.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("figure5met.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("figure5star.png") grid.raster(img) ## ----,eval=FALSE,echo=TRUE,message=FALSE,warning=FALSE------------------- ## #----------------------------------- ## # STEP 1: Search, download, prepare | ## #----------------------------------- ## # Step 1.1 download methylation data| ## # ----------------------------------- ## path <- "." ## query <- TCGAquery(tumor = "KIRC",level = 3, platform = "HumanMethylation450") ## TCGAdownload(query, path = path) ## kirc.met <- TCGAprepare(query,dir = path, ## save = TRUE, ## filename = "metKirc.rda", ## summarizedExperiment = FALSE) ## ## kirc.met <- TCGAprepare_elmer(kirc.met, ## platform = "HumanMethylation450", ## save = TRUE, ## met.na.cut = 0.2) ## ## # Step 1.2 download expression data ## query.rna <- TCGAquery(tumor="KIRC",level=3, platform="IlluminaHiSeq_RNASeqV2") ## TCGAdownload(query.rna,path=path,type = "rsem.genes.normalized_results") ## kirc.exp <- TCGAprepare(query.rna, dir=path, save = TRUE, ## type = "rsem.genes.normalized_results", ## filename = "expKirc.rda", summarizedExperiment = FALSE) ## ## kirc.exp <- TCGAprepare_elmer(kirc.exp, ## save = TRUE, ## platform = "IlluminaHiSeq_RNASeqV2") ## #----------------------------------- ## # STEP 2: ELMER integration | ## #----------------------------------- ## # Step 2.1 prepare mee object | ## # ----------------------------------- ## library(ELMER) ## library(parallel) ## ## geneAnnot <- txs() ## geneAnnot$GENEID <- paste0("ID",geneAnnot$GENEID) ## geneInfo <- promoters(geneAnnot,upstream = 0, downstream = 0) ## probe <- get.feature.probe() ## mee <- fetch.mee(meth = kirc.met, exp = kirc.exp, TCGA = TRUE, ## probeInfo = probe, geneInfo = geneInfo) ## ## #-------------------------------------- ## # STEP 3: Analysis | ## #-------------------------------------- ## # Step 3.1: Get diff methylated probes | ## #-------------------------------------- ## Sig.probes <- get.diff.meth(mee ,cores=detectCores(), ## dir.out ="kirc",diff.dir="hypo", ## pvalue = 0.01) ## ## #-------------------------------------------------- ## # Step 3.2: Identifying putative probe-gene pairs | ## #-------------------------------------------------- ## # Collect nearby 20 genes for Sig.probes ## nearGenes <- GetNearGenes(TRange=getProbeInfo(mee, probe=Sig.probes), ## cores=detectCores(), ## geneAnnot=getGeneInfo(mee)) ## ## # Identify significant probe-gene pairs ## Hypo.pair <- get.pair(mee=mee, ## probes=Sig.probes, ## nearGenes=nearGenes, ## permu.dir="./kirc/permu", ## dir.out="./kirc/", ## cores=detectCores(), ## label= "hypo", ## permu.size=10000, ## Pe = 0.001) ## ## Sig.probes.paired <- fetch.pair(pair=Hypo.pair, ## probeInfo = getProbeInfo(mee), ## geneInfo = getGeneInfo(mee)) ## ## #------------------------------------------------------------- ## # Step 3.3: Motif enrichment analysis on the selected probes | ## #------------------------------------------------------------- ## enriched.motif <- get.enriched.motif(probes=Sig.probes.paired, ## dir.out="kirc", label="hypo", ## background.probes = probe$name) ## ## #------------------------------------------------------------- ## # Step 3.4: Identifying regulatory TFs | ## #------------------------------------------------------------- ## TF <- get.TFs(mee=mee, ## enriched.motif=enriched.motif, ## dir.out="kirc", ## cores=detectCores(), label= "hypo") ## ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("case4_elmer.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("elmer1.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("elmer2.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("elmer3.png") grid.raster(img) ## ----, fig.width=6, fig.height=4, echo = FALSE, fig.align="center",hide=TRUE, message=FALSE,warning=FALSE---- library(png) library(grid) img <- readPNG("elmer4.png") grid.raster(img) ## ----sessionInfo--------------------------------------------------------- sessionInfo()