Li's Bioinfo-Blog
  • |
  • 主页
  • 分类
  • 标签
  • 归档
  • 关于
  • 搜索
Home » 分类

📖 生信数据分析--分析流程,工具包等

WGCNA基因加权共表达网络分析

1、关于WGCNA原理 1.1 建立共表达网络 在基因共表达网络中,节点node代表基因,边edge代表两个基因间共表达关系。 若一个基因同时与多个基因存在相关性,称为hub基因。 若一群基因存在高度互相相关,称为module。 基因共表达网络的展示形式一般为 n × n 邻接矩阵adjacency matrix(n个基因) ...

Create:&nbsp;<span title='2022-04-25 00:00:00 +0000 UTC'>2022-04-25</span>&nbsp;|&nbsp;Update:&nbsp;2022-05-06&nbsp;|&nbsp;Words:&nbsp;5111&nbsp;|&nbsp;11 min&nbsp;|&nbsp;Lishensuo

GSVA包单样本富集分析

1、GSVA函数 1 2 3 4 5 6 7 # BiocManager::install("GSVA") library(GSVA) ?gsva gsva(expr = , #metrix格式表达矩阵(行--基因,列--样本) gset.idx.list = , #list格式基因集 method=c("gsva", "ssgsea", "zscore", "plage"), # defaul:gsva kcdf=c("Gaussian", "Poisson", "none")) # default:Gaussian 如上所示,gsva()函数有4个关键参数。 前两个参数分别提供表达矩阵与基因集,其中基因名需要为相同的格式。 method参数设置方法,默认为gsva kcdf参数设置分布模拟类型,默认为Gaussian,适合于array数据以及log标准化(CPM/RPKM/TPM等)的RNA-seq数据。Poisson分布适合RNA-seq的原始count矩阵。 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 ## (1) 表达矩阵 # BiocManager::install("GSVAdata") library(GSVAdata) data(commonPickrellHuang) ### array矩阵 array_mt = exprs(huangArrayRMAnoBatchCommon_eset) dim(array_mt) # [1] 11508 36 array_mt[1:4,1:4] # NA19099 NA18523 NA19144 NA19137 # 8567 8.370526 8.544890 8.270395 8.585984 # 23139 7.536228 7.232781 7.469420 7.801986 # 7580 6.325102 6.405961 6.510249 6.285510 # 55619 9.134495 9.049938 9.704679 9.285202 ### count矩阵 count_mt = exprs(pickrellCountsArgonneCQNcommon_eset) dim(count_mt) # [1] 11508 36 count_mt[1:4,1:4] # NA19099 NA18523 NA19144 NA19137 # 8567 326 209 318 343 # 23139 255 169 245 361 # 7580 72 69 124 76 # 55619 487 590 678 540 ## (2) 基因集 library(tidyverse) library(msigdbr) hm_H = msigdbr(species = "Homo sapiens", category = "H", subcategory = NULL) %>% as.data.frame() %>% dplyr::select(gs_cat, gs_subcat, gs_name, entrez_gene) hm_H.list = split(hm_H$entrez_gene, hm_H$gs_name) str(head(hm_H.list)) # List of 6 # $ HALLMARK_ADIPOGENESIS : int [1:210] 19 11194 10449 33 34 35 47 50 51 112 ... # $ HALLMARK_ALLOGRAFT_REJECTION: int [1:335] 16 6059 10006 43 92 207 322 567 567 586 ... # $ HALLMARK_ANDROGEN_RESPONSE : int [1:102] 10257 11057 2181 87 9510 11047 9590 207 220 56172 ... # $ HALLMARK_ANGIOGENESIS : int [1:36] 350 351 894 1281 1290 6372 2260 11167 3685 182 ... # $ HALLMARK_APICAL_JUNCTION : int [1:231] 58 60 70 71 72 87 88 89 81 81 ... # $ HALLMARK_APICAL_SURFACE : int [1:46] 102 79602 79602 84632 9465 351 50617 5205 2683 672 ... 3、GSVA分析 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 27 ## (1) array + gsva res_gsva = gsva(array_mt, hm_H.list) class(res_gsva) # [1] "matrix" "array" dim(res_gsva) # [1] 50 36 res_gsva[1:4,1:4] # NA19099 NA18523 NA19144 NA19137 # HALLMARK_ADIPOGENESIS -0.21147644 0.08157479 0.004672524 -0.1656925 # HALLMARK_ALLOGRAFT_REJECTION -0.32620137 -0.44105423 -0.452996166 -0.4288909 # HALLMARK_ANDROGEN_RESPONSE 0.02056377 0.09123573 -0.135386552 -0.2198959 # HALLMARK_ANGIOGENESIS -0.02435541 -0.18282244 0.281235763 0.5114103 ## (2) count + ssgsea res_gsva = gsva(count_mt, hm_H.list, method = "ssgsea", kcdf = "Poisson") class(res_gsva) # [1] "matrix" "array" dim(res_gsva) # [1] 50 36 res_gsva[1:4,1:4] # NA19099 NA18523 NA19144 NA19137 # HALLMARK_ADIPOGENESIS 0.3288768 0.32899026 0.2893928 0.3297140 # HALLMARK_ALLOGRAFT_REJECTION 0.5018248 0.49183093 0.4672342 0.4531009 # HALLMARK_ANDROGEN_RESPONSE 0.3210905 0.28289473 0.2818011 0.2845855 # HALLMARK_ANGIOGENESIS 0.1794361 0.07705857 0.2076328 0.1276449 4、limma差异分析 以上面array + gsva的分析结果为例,随机模拟分组进行limma差异分析 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 27 ## 模拟分组 meta = data.frame(sample=colnames(array_mt), group=rep(c("G1","G2"), each=18)) head(meta) # sample group # 1 NA19099 G1 # 2 NA18523 G1 # 3 NA19144 G1 # 4 NA19137 G1 # 5 NA18861 G1 # 6 NA19116 G1 ## limma差异分析 library(limma) mod <- model.matrix(~ factor(meta$group)) colnames(mod) <- c("G1", "G2") fit <- lmFit(res_gsva, mod) fit <- eBayes(fit) tt <- topTable(fit, coef=2, n=Inf) head(tt) # logFC AveExpr t P.Value adj.P.Val B # HALLMARK_COMPLEMENT 0.08684907 -0.10731154 2.495909 0.01676849 0.3324198 -3.059104 # HALLMARK_HEME_METABOLISM 0.07942485 -0.04250442 2.238926 0.03076698 0.3324198 -3.486635 # HALLMARK_IL2_STAT5_SIGNALING 0.07839683 -0.09276742 2.193525 0.03411795 0.3324198 -3.558694 # HALLMARK_MITOTIC_SPINDLE -0.10721385 -0.07656961 -2.057416 0.04617883 0.3324198 -3.768076 # HALLMARK_IL6_JAK_STAT3_SIGNALING 0.09342746 -0.20104143 2.026107 0.04943145 0.3324198 -3.814790 # HALLMARK_INFLAMMATORY_RESPONSE 0.10211588 -0.10571870 1.996939 0.05263925 0.3324198 -3.857810

Create:&nbsp;<span title='2023-01-15 00:00:00 +0000 UTC'>2023-01-15</span>&nbsp;|&nbsp;Update:&nbsp;2023-01-15&nbsp;|&nbsp;Words:&nbsp;872&nbsp;|&nbsp;2 min&nbsp;|&nbsp;Lishensuo

xCell与CIBERSORT等免疫浸润分析

xCell xCell包 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 devtools::install_github('dviraran/xCell') library(xCell) data("xCell.data") ##查看支持的64种细胞类型,同下图 colnames(xCell.data$spill$K) ##预测函数的关键参数解释 ?xCellAnalysis() # expr = 交代表达矩阵; ##如果是array,不需要额外标准化;如果是RNAseq,需要TPM/FPKM/TPM。 ##对于基因ID格式需要是symbol格式。 # rnaseq = TRUE 数据是否为RNAseq数据,如果是芯片数据设置为FALSE # cell.types.use = NULL 提供一个字符串,说明想要预测64种细胞中的哪些细胞类型 # parallel.sz = 4 调用的线程数,默认为4 NOTE: ...

Create:&nbsp;<span title='2022-04-25 00:00:00 +0000 UTC'>2022-04-25</span>&nbsp;|&nbsp;Update:&nbsp;2023-02-18&nbsp;|&nbsp;Words:&nbsp;2554&nbsp;|&nbsp;6 min&nbsp;|&nbsp;Lishensuo

ClusterGVis包绘制基因表达矩阵热图

ClusterGVis包是中国药科大学Jun Zhang博士开发的系列可视化工具包之一,可以基因表达矩阵进行高级的热图可视化分析。如下根据其github以及微信教程简单整理一下自己感兴趣的用法。 ...

Create:&nbsp;<span title='2023-03-19 00:00:00 +0000 UTC'>2023-03-19</span>&nbsp;|&nbsp;Update:&nbsp;2023-03-19&nbsp;|&nbsp;Words:&nbsp;1261&nbsp;|&nbsp;3 min&nbsp;|&nbsp;Lishensuo

survival包生存分析及glmnet包lasso回归

生存分析(survival analysis)的主要目的是发现与患者生存事件相关的指标因素,例如年龄性别、基因表达/突变等。如下学习相关基础知识及几种常见的生存分析方法。 ...

Create:&nbsp;<span title='2023-01-28 00:00:00 +0000 UTC'>2023-01-28</span>&nbsp;|&nbsp;Update:&nbsp;2023-01-29&nbsp;|&nbsp;Words:&nbsp;3794&nbsp;|&nbsp;8 min&nbsp;|&nbsp;Lishensuo

MuSiC包根据scRNAseq预测Bulk细胞组成

MuSiC(MUlti-Subject SIngle Cell deconvolution)是来自宾夕法尼亚大学Biostatistics, Epidemiology and Informatics系的Mingyao Li课题组于2019年发表于Nature Communication的一个工具R包,可根据单细胞转录组信息推测Bulk RNA-seq细胞组成。而后,该团队又于2022年在Briefing in bioinformatics发表了扩展版本MuSiC2,可以考虑更复杂的场景。 ...

Create:&nbsp;<span title='2023-03-26 00:00:00 +0000 UTC'>2023-03-26</span>&nbsp;|&nbsp;Update:&nbsp;2023-03-26&nbsp;|&nbsp;Words:&nbsp;1799&nbsp;|&nbsp;4 min&nbsp;|&nbsp;Lishensuo

从RNAseq的fastq.gz提取表达矩阵

RNA-seq数据比对流程主要分为三步(1)整理数据;(2)质控;(3)比对。其中每一步都涉及到若干软件的用法,如下简单整理基本的分析流程。 示例数据:GSE158623中6个样本的RNA-seq测序结果(human),对应SRR12720999~SRR12721004 ...

Create:&nbsp;<span title='2023-01-29 00:00:00 +0000 UTC'>2023-01-29</span>&nbsp;|&nbsp;Update:&nbsp;2023-01-29&nbsp;|&nbsp;Words:&nbsp;2347&nbsp;|&nbsp;5 min&nbsp;|&nbsp;Lishensuo

gseapy富集分析与可视化

参考 https://gseapy.readthedocs.io/en/latest/index.html 示例数据来自 https://github.com/zqfang/GSEApy/tree/master/tests/data 1 2 3 4 5 6 7 8 import os import pandas as pd import scanpy as sc import gseapy as gp from gseapy import barplot, dotplot import matplotlib.pyplot as plt 0. 基因集获取 Enrichr API supported organisms: Human, Mouse, Yeast, Fly, Fish, Worm 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 27 28 29 30 31 32 names = gp.get_library_name(organism = "Human") len(names) # 217 [name for name in names if "GO" in name] # ['GO_Biological_Process_2021', # 'GO_Biological_Process_2023', # 'GO_Biological_Process_2025', # 'GO_Cellular_Component_2021', # 'GO_Cellular_Component_2023', # 'GO_Cellular_Component_2025', # 'GO_Molecular_Function_2021', # 'GO_Molecular_Function_2023', # 'GO_Molecular_Function_2025', # 'SynGO_2022', # 'SynGO_2024'] [name for name in names if "KEGG" in name] # ['KEGG_2013', # 'KEGG_2015', # 'KEGG_2016', # 'KEGG_2019_Human', # 'KEGG_2019_Mouse', # 'KEGG_2021_Human'] [name for name in names if "MSigDB" in name] # ['MSigDB_Computational', 'MSigDB_Hallmark_2020', 'MSigDB_Oncogenic_Signatures'] kegg = gp.get_library(name='KEGG_2021_Human', organism='Human') type(kegg) # dict kegg['AMPK signaling pathway'][:4] # ['RAB2A', 'PPP2R1A', 'TSC2', 'TSC1'] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 # 还支持对MSigDB基因集数据库的获取 from gseapy import Msigdb msig = Msigdb() # all available versions msig.list_dbver() # all types of gene sets given a dbver msig.list_category(dbver="2025.1.Hs") # all gene sets of one type gmt = msig.get_gmt(category='h.all', dbver="2025.1.Hs") type(gmt) # dict gmt.keys() # hallmark pathway names gmt["HALLMARK_ADIPOGENESIS"] # gene list 1. ORA富集 1.1 分析 Online analysis 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 ## 1) 两个最重要的参数 (后同) # gene_list gene_list="tests/data/gene_list.txt" # file path, one gene per row gene_list= ["RAB2A","PRKAB2","PRKAA1"] # list of genes # gene_sets (from enrichr API) gene_sets='KEGG_2016' gene_sets=['KEGG_2021_Human','GO_Biological_Process_2025'] gene_list = pd.read_csv("tests/data/gene_list.txt", header=None)[0].tolist() gene_list[:4] # ['IGKV4-1', 'CD55', 'IGKC', 'PPFIBP1'] enr = gp.enrichr(gene_list=gene_list, # or "./tests/data/gene_list.txt", gene_sets=['GO_Biological_Process_2025'], organism='human', outdir=None, ) enr.results.head(5) # Index(['Gene_set', 'Term', 'Overlap', 'P-value', 'Adjusted P-value', # 'Old P-value', 'Old Adjusted P-value', 'Odds Ratio', 'Combined Score', # 'Genes'], # dtype='object') ## 2) 设置结果导出路径 enr = gp.enrichr(gene_list=gene_list, gene_sets=['GO_Biological_Process_2025'], organism='human', outdir="enrichr_kegg", ) os.listdir("enrichr_kegg") # ['GO_Biological_Process_2025.human.enrichr.reports.txt', # 'GO_Biological_Process_2025.human.enrichr.reports.pdf'] ## 3) 自定义Background genes (P值会存在一定差异) enr = gp.enrichr(gene_list=gene_list, gene_sets=['GO_Biological_Process_2025'], organism='human', background="tests/data/background.txt", # 自定义背景基因 / or list object outdir=None, ) Offline analysis (本地提供通路基因集数据) 1 2 3 4 5 6 7 8 9 10 11 from gseapy.parser import read_gmt gmt_dict = read_gmt("tests/data/c2.cp.kegg.v7.5.1.symbols.gmt") gmt_dict.keys() enr2 = gp.enrich(gene_list=gene_list, # or gene_list=glist gene_sets=["tests/data/c2.cp.kegg.v7.5.1.symbols.gmt"], background=None, outdir=None, verbose=True) enr2.results.head(5) 1.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 27 # ofname 设置图片的导出路径,若为None则不保存。也支持后续使用plt语法进一步调整图片 ax = dotplot(enr.results, column="Adjusted P-value", x='Gene_set', # set x axis, so you could do a multi-sample/library comparsion size=5, cutoff=0.05, top_term=10, figsize=(3, 5), title = "GO_BP", xticklabels_rot=45, # rotate xtick labels show_ring=True, # set to False to revmove outer ring marker='o', ofname='./plot/p1.png' ) ax = barplot(enr.results, column="Adjusted P-value", group='Gene_set', # set group, so you could do a multi-sample/library comparsion size=10, top_term=10, figsize=(3,5), #color=['darkred', 'darkblue'] # set colors for group color = {'GO_Biological_Process_2025': 'salmon'}, ofname='./plot/p2.png' ) 1 2 3 ax = dotplot(enr.res2d, title='Enrichr',cmap='viridis_r', size=10, figsize=(3,5), ofname='./plot/p3.png') ax = barplot(enr.res2d, title='Enrichr',cmap='viridis_r', size=4, figsize=(4, 5), color='darkred', ofname='./plot/p4.png') # 差别不大 2. Prerank富集 基于差异分析结果的GSEA ...

Create:&nbsp;<span title='2025-07-22 00:00:00 +0000 UTC'>2025-07-22</span>&nbsp;|&nbsp;Update:&nbsp;2025-07-22&nbsp;|&nbsp;Words:&nbsp;1367&nbsp;|&nbsp;3 min&nbsp;|&nbsp;Lishensuo

生信相关网站数据库集锦

1、HPA HPA:the Human Protein Atlas,由2003年来自瑞典的科研机构发起,旨在绘制综合性人类蛋白质图谱。 https://www.proteinatlas.org/ 蛋白质表达数据库,常见用途包括: (1)查看特定基因在不同组织、脑区,细胞类型,组织细胞类型,疾病(癌症),免疫细胞,肿瘤细胞系等表达情况。 (2)数据挖掘类文章常使用其进行比较基因在肿瘤部位与相应正常部位的蛋白水平表达差异。 2、ENCORI ENCORI,The Encyclopedia of RNA Interactomes,由中山大学生命科学学院屈良鹄团队开发,于2014年发表于Nucleic Acids Res。 https://starbase.sysu.edu.cn/index.php miRNA、lncRNA、RBP等多维相互作用网络,常见用途包括 miRNA/RBP的靶标查询(mRNA,lncRNA…) RNA interaction, ceRNA network TCGA肿瘤的差异表达,生存分析,相关性等 3、TCIA TCIA,The Cancer Immunome Database,由来自奥地利的因斯布鲁克大学医学院Zlatko Trajanoski团队开发,于2017年发表于Cell Reports https://tcia.at/home 20种solid cancer的免疫相关分析,例如 免疫基因表达、细胞浸润、肿瘤亚克隆等 亮点之一是提出Immunophenoscore指标用于预测免疫治疗响应 https://github.com/icbi-lab/Immunophenogram 样本临床信息也整理的较为完整 ...

Create:&nbsp;<span title='2023-04-02 00:00:00 +0000 UTC'>2023-04-02</span>&nbsp;|&nbsp;Update:&nbsp;2023-04-02&nbsp;|&nbsp;Words:&nbsp;832&nbsp;|&nbsp;2 min&nbsp;|&nbsp;Lishensuo

TCGA的SNV数据下载与maftools可视化

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")

Create:&nbsp;<span title='2023-04-09 00:00:00 +0000 UTC'>2023-04-09</span>&nbsp;|&nbsp;Update:&nbsp;2023-04-09&nbsp;|&nbsp;Words:&nbsp;608&nbsp;|&nbsp;2 min&nbsp;|&nbsp;Lishensuo
« Prev Page Next Page »
© 2025 Li's Bioinfo-Blog Powered by Hugo & PaperMod
您是本站第 位访问者,总浏览量为 次