refgenie：参考基因组（阿拉丁）商店 http://refgenie.databio.org/ Here we provide a web interface and a RESTful API to access genome assets for popular reference genome assemblies. 该平台由位于弗吉尼亚大学公共卫生基因组学中心的计算生物学和生物信息学研究小组(Sheffield lab of computational biology)建立。上次修改/更新是2021年11月。 ...
1 2 3 4 5 if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("STRINGdb") library(STRINGdb) 1、定义要使用的STRING版本、物种，以及PPI阈值分数 1 2 3 4 string_db <- STRINGdb$new(version="11", species=9606, score_threshold=200, input_directory="") 2、示例基因 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 data(diff_exp_example1) genes = diff_exp_example1[1:50,] head(genes) # pvalue logFC gene # 1 0.0001018 3.333461 VSTM2L # 2 0.0001392 3.822383 TBC1D2 # 3 0.0001720 3.306056 LENG9 # 4 0.0001739 3.024605 TMEM27 # 5 0.0001990 3.854414 LOC100506014 # 6 0.0002393 3.082052 TSPAN1 ###基因名匹配protein ID #第一个参数是data.frame; 第二个参数是基因所在列的列名 genes_mapped <- string_db$map(genes, "gene" ) #Warning: we couldn't map to STRING 30% of your identifiers head(genes_mapped) # gene pvalue logFC STRING_id # 1 VSTM2L 0.0001018 3.333461 9606.ENSP00000362560 # 2 TBC1D2 0.0001392 3.822383 9606.ENSP00000481721 # 3 LENG9 0.0001720 3.306056 9606.ENSP00000479355 # 4 TMEM27 0.0001739 3.024605 9606.ENSP00000369699 # 40 LOC100506014 0.0001990 3.854414 <NA> # 5 TSPAN1 0.0002393 3.082052 9606.ENSP00000361072 #string_db$plot_network(genes_mapped$STRING_id) 3、下载这些基因间的互作关系 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ppi = string_db$get_interactions(genes_mapped$STRING_id) %>% distinct() ppi = ppi %>% dplyr::left_join(genes_mapped[,c(1,4)], by=c('from'='STRING_id')) %>% #列标序号根据具体情况而定 dplyr::rename(Gene1=gene) %>% ##列名根据具体情况而定 dplyr::left_join(genes_mapped[,c(1,4)], by=c('to'='STRING_id')) %>% dplyr::rename(Gene2=gene) %>% dplyr::select(Gene1, Gene2, combined_score) head(ppi) # Gene1 Gene2 combined_score # 1 C3 TYROBP 240 # 2 ABCA12 GRHL3 308 # 3 FAM189A1 TM4SF20 400 # 4 ABCA12 NIPAL4 824 # 5 GRHL3 NIPAL4 275 # 6 GRHL3 IGDCC4 238
1、TCGAbiolinks下载数据 使用TCGAbiolinks下载特定肿瘤类型的SNV数据 https://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/mutation.html 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 library(TCGAbiolinks) query <- GDCquery( project = "TCGA-CHOL", data.category = "Simple Nucleotide Variation", access = "open", legacy = FALSE, data.type = "Masked Somatic Mutation") GDCdownload(query) maf <- GDCprepare(query) dim(maf) # [1] 3764 141 ## (1) 因后续需要，修改Tumor_Sample_Barcode列 maf$long_Barcode = maf$Tumor_Sample_Barcode maf$Tumor_Sample_Barcode = substr(maf$Tumor_Sample_Barcode,1,12) length(unique(maf$Tumor_Sample_Barcode)) # 51 ## (2) 读取临床生存数据 clinical = readxl::read_xlsx("TCGA_Pan_Cancer_Clinical_Data_mmc1.xlsx") clinical_sle = clinical %>% dplyr::filter(type=="CHOL") %>% dplyr::select(bcr_patient_barcode, OS, OS.time, clinical_stage) %>% dplyr::rename(Tumor_Sample_Barcode=bcr_patient_barcode) dim(clinical_sle) # 45 2、maftools可视化 https://bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html https://www.jieandze1314.com/post/cnposts/237/ 1 2 3 4 5 6 7 8 9 library(maftools) maf_obj = read.maf(maf = maf, clinicalData = clinical_sle) #每个样本的突变情况统计 getSampleSummary(maf_obj) #每个基因的突变类型统计 getGeneSummary(maf_obj) 2.1概括图 1 2 3 plotmafSummary(maf = maf_obj, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE) 2.2 基因瀑布图 1 2 3 4 5 6 7 8 9 10 11 12 13 oncoplot(maf = maf_obj, top = 10) # 如下图 ## 添加临床注释 oncoplot(maf = maf_obj, top = 10, clinicalFeatures = c('clinical_stage',"OS"), draw_titv = TRUE) ## 选择特定基因集 set.seed(42) genes = sample(maf_obj@data$Hugo_Symbol,10) oncoplot(maf = maf_obj, genes = genes, clinicalFeatures = c('clinical_stage',"OS"), draw_titv = TRUE) 2.3 转换颠倒统计 1 2 3 4 5 # Transition 转换 ： 嘌呤(AG)或嘧啶(CT)内部之间转换 # Transversions 颠倒：嘌呤与嘧啶间互相转换 maf_obj.titv = titv(maf = maf_obj, plot = FALSE, useSyn = TRUE) plotTiTv(res = maf_obj.titv) 2.4 基因对突变统计 1 2 3 4 5 6 7 8 # green: co-occuring # yellow: mutually exclusive somaticInteractions(maf = maf_obj, top = 25, pvalue = c(0.05, 0.1)) set.seed(42) genes = sample(maf_obj@data$Hugo_Symbol,25) somaticInteractions(maf = maf_obj, genes = genes , pvalue = c(0.05, 0.1)) 2.5 生存分析 根据特定基因是否突变将病人分成WT与Mutant两组 1 2 3 4 5 6 7 8 9 10 11 mafSurvival(maf = maf_obj, genes = 'TP53', time = 'OS.time', Status = 'OS') # Group medianTime N # 1: Mutant 732 4 # 2: WT 650 41 ## 提取信息 # maf_obj@clinical.data %>% # dplyr::mutate(Group=ifelse(Tumor_Sample_Barcode %in% # subset(maf_obj@data, Hugo_Symbol=="TP53")$Barcode, # "Mutant","WT")) 2.6 基因对的生存相关性 1 2 3 4 5 6 7 8 9 10 prog_geneset = survGroup(maf = maf_obj, top = 200, geneSetSize = 2, time = "OS.time", Status = "OS", verbose = FALSE,minSamples = 3) prog_geneset # Gene_combination P_value hr WT Mutant # 1: PBRM1_PLXNA4 0.243 2.36e+00 42 3 # 2: PBRM1_PCLO 0.294 3.46e-01 42 3 # 3: PBRM1_TP53 0.320 3.64e-01 42 3 mafSurvGroup(maf = maf_obj, geneSet = c("PBRM1", "PLXNA4"), time = "OS.time", Status = "OS")
UCSCXenaShiny是基于集成了多种肿瘤数据库的UCSCXena平台，进行数据下载、分析、可视化的Shiny工具（以及同名R包），由上海科技大学王诗翔博士等共同开发；于2021年6月发表于Bioinformatics。下面主要学习其R包的相关函数，了解其核心功能。 ...
Github仓库：https://github.com/openbiox/UCSCXenaShiny Online App：https://shiny.zhoulab.ac.cn/UCSCXenaShiny/ ...