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： ...
ClusterGVis包是中国药科大学Jun Zhang博士开发的系列可视化工具包之一，可以基因表达矩阵进行高级的热图可视化分析。如下根据其github以及微信教程简单整理一下自己感兴趣的用法。 ...
MuSiC(MUlti-Subject SIngle Cell deconvolution)是来自宾夕法尼亚大学Biostatistics, Epidemiology and Informatics系的Mingyao Li课题组于2019年发表于Nature Communication的一个工具R包，可根据单细胞转录组信息推测Bulk RNA-seq细胞组成。而后，该团队又于2022年在Briefing in bioinformatics发表了扩展版本MuSiC2，可以考虑更复杂的场景。 ...
RNA-seq数据比对流程主要分为三步（1）整理数据；（2）质控；（3）比对。其中每一步都涉及到若干软件的用法，如下简单整理基本的分析流程。 示例数据：GSE158623中6个样本的RNA-seq测序结果(human)，对应SRR12720999~SRR12721004 ...
参考 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 ...
ConsensusClusterPlus包是肿瘤分型研究的常用工具，其于2010年发表于Bioinformatics。 Paper：https://academic.oup.com/bioinformatics/article/26/12/1572/281699 Tutorial：https://bioconductor.org/packages/release/bioc/vignettes/ConsensusClusterPlus/inst/doc/ConsensusClusterPlus.pdf 1 2 # BiocManager::install("ConsensusClusterPlus") library(ConsensusClusterPlus) 1、示例表达矩阵 行名为基因，列名为样本的基因表达矩阵；根据需要，进行标准化处理 基因集的选择是关键的一步，可根据统计学（高变基因）或者生物学（特定功能相关基因）进行选择 如下是参考教程的数据 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 library(ALL) data(ALL) d=exprs(ALL) dim(d) # [1] 12625 128 mads=apply(d,1,mad) d=d[rev(order(mads))[1:5000],] # 筛选高方差基因 dim(d) # [1] 5000 128 ## 标准化处理：每行（基因）减去每行的中位数 d = sweep(d,1, apply(d,1,median,na.rm=T)) d[1:4,1:4] # 01005 01010 03002 04006 # 36638_at 1.556121 0.9521271 -0.05018082 4.780378 # 39318_at 1.191353 2.5013225 -2.38793537 -1.199521 # 38514_at 1.020716 3.2785671 1.55949145 -3.345919 # 266_s_at 1.829260 0.3624327 1.54913247 -1.286294 2、样本亚型鉴定 实际就是一个函数即可：ConsensusClusterPlus() 其中涉及到较多参数的选择，具体如下 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 ## 默认参数 ConsensusClusterPlus(d, maxK=3, reps=10, pItem=0.8, pFeature=1, title="untitled_consensus_cluster", clusterAlg="hc", distance="pearson", plot=NULL, writeTable=FALSE, seed=42) ## 参数含义 # maxK：考虑的最大聚类数，建议取10~20 # reps：抽样次数，建议1000 # pItem与pFeature：分别表示对样本与基因的抽样比例 # title：图片或文件的保存路径 # clusterAlg：聚类方式 c("hc","pam","km") # distance: 距离计算方式 c("pearson","spearman","euclidean","binay","maximum","canberra","minkowski") # plot : 是否绘图，以及图形类型 c(NULL, "pdf", "png", "pngBMP") # writeTable: 是否保存文件 # seed：随机种子 results = ConsensusClusterPlus(d, maxK=10, reps=1000, pItem=0.8, pFeature=1, title="./tmp/", clusterAlg="hc", distance="pearson", plot="pdf", writeTable=FALSE, seed=42) R对象结果：list格式 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ## 查看聚类数为2的结果信息 names(results[[3]]) # [1] "consensusMatrix" "consensusTree" "consensusClass" "ml" "clrs" table(results[[3]]$consensusClass) # 1 2 3 # 69 28 31 dim(results[[3]]$consensusMatrix) # [1] 128 128 results[[3]]$consensusMatrix[1:4,1:4] # [,1] [,2] [,3] [,4] # [1,] 1.0000000 0.3408360 0.7996870 0.3515249 # [2,] 0.3408360 1.0000000 0.1224806 1.0000000 # [3,] 0.7996870 0.1224806 1.0000000 0.1190108 # [4,] 0.3515249 1.0000000 0.1190108 1.0000000 ## consensus values 0 (never clustered together) ## consensus values 1 (always clustered together) 图形输出结果，主要两类图 （1）特定聚类结果的热图 ...
https://bioconductor.org/packages/release/bioc/html/IsoformSwitchAnalyzeR.html 1、背景知识 The usage of Alternative Transcription Start sites (aTSS可选择转录起始位点), Alternative Splicing (AS可选择剪切位点) and alternative Transcription Termination Sites (aTTS可选择终止位点) are collectively collectively results in the production of different isoforms. Alternative isoforms are widely used as recently demonstrated by The ENCODE Consortium, which found that on average, 6.3 different transcripts are generated per gene; a number which may vary considerably per gene. ...
肿瘤预后类文章的常规步骤之一是绘制诺模图，并进行相关分析。以下总结了相关基础绘制工具。 示例数据集 1 2 3 4 5 6 library(survival) head(lung) # inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss # 1 3 306 2 74 1 1 90 100 1175 NA # 2 3 455 2 68 1 0 90 90 1225 15 # 3 3 1010 1 56 1 0 90 90 NA 15 1、rms包 参考用法：https://atm.amegroups.com/article/view/14736/15089 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 library(rms) # var.labels = c(age="Age in Years", # lac="lactate", # sex="Sex of the participant", # shock="shock", # y="outcome", # Y="ordinal") # label(data) = lapply(names(var.labels), # function(x) label(data[,x]) = var.labels[x]) mod.cox <- cph(Surv(time,status) ~ ph.ecog+sex+age,lung, surv=TRUE) ddist <- datadist(lung) options(datadist='ddist') surv.cox <- Survival(mod.cox) # 定义函数 med <- Quantile(mod.cox) surv <- Survival(mod.cox) nom.cox <- nomogram(mod.cox, # 根据total point进行特定函数计算 fun=list(function(x) surv.cox(365, x), # 一年生存率 function(x) med(lp=x, q=0.5)), # 中位生存时间 funlabel=c("200-Day Sur. Prob.", "Median Survival Time"), lp=F, # 不展示Linear Predictor conf.int=c(0.1,0.5) # 两个置信区间 ) plot(nom.cox, col.conf=c('red','green'), # 置信区间的颜色 col.grid = c("grey30","grey") # 网格的颜色 ) # f = cph(Surv(time, status) ~ age + sex + ph.karno, data = lung, # x = T, y = T, sur = T) # pred_score = apply(lung, 1, function(x){ # pred=Predict(f, age=x["age"], sex=x["sex"], ph.karno=x["ph.karno"]) # return(pred$yhat) # }) %>% unlist() # summary(pred_score) 校准曲线 ...
EnhancedVolcano包可根据差异分析结果，基于ggplot2绘图结构，方便地绘制美观的火山图，下面根据自己的理解小结下基本用法。 官方全面的教程：https://github.com/kevinblighe/EnhancedVolcano 示例差异基因数据 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 library(airway) library(magrittr) data('airway') airway$dex %<>% relevel('untrt') ens <- rownames(airway) library(org.Hs.eg.db) symbols <- mapIds(org.Hs.eg.db, keys = ens, column = c('SYMBOL'), keytype = 'ENSEMBL') symbols <- symbols[!is.na(symbols)] symbols <- symbols[match(rownames(airway), names(symbols))] rownames(airway) <- symbols keep <- !is.na(rownames(airway)) airway <- airway[keep,] library('DESeq2') dds <- DESeqDataSet(airway, design = ~ cell + dex) dds <- DESeq(dds, betaPrior=FALSE) res <- results(dds, contrast = c('dex','trt','untrt')) res <- lfcShrink(dds, contrast = c('dex','trt','untrt'), res=res, type = 'normal') res <- as.data.frame(res) head(res) # baseMean log2FoldChange lfcSE stat pvalue padj # TSPAN6 710.0931707 -0.37807189 0.09851236 -3.8404448 0.0001228116 0.0009522932 # TNMD 0.0000000 NA NA NA NA NA # DPM1 521.2572396 0.19826365 0.10931684 1.8155169 0.0694445184 0.1910397405 # SCYL3 237.6068046 0.03234467 0.13821470 0.2371917 0.8125081096 0.9118161375 # C1orf112 58.0358739 -0.08835419 0.25056704 -0.3194810 0.7493618190 0.8773885438 # FGR 0.3194343 -0.08459224 0.15186225 -0.3948862 0.6929268648 NA 如上，只要包含包含基因名、差异倍数、P值三部分信息的差异结果就可以用于绘制火山图。 ...
IOBR包集signature打分与免疫浸润分析为一体的肿瘤数据分析工具，由南方医科大学南方医院廖旺军教授，曾东强博士等人于2021年7月发表于Frontiers in Immunology，引用数以超过100余次。现根据其github教程学习其部分功能用法。 ...