pbcmc包的介绍（根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理，感谢Jimmy老师！）

  • 1.使用BioConductor的BreastCancerXXX数据集
  • 2.使用任何微阵列R数据包
  • 3.用PAM50 centroids作为验证

pbcmc: Permutation-Based Confidence for Molecular Classification
Bioconductor link: https://bioconductor.riken.jp/packages/3.3/bioc/html/pbcmc.html
亮点：Pbcmc包的特点是基于排列测试（置换检验，permutation test）对基因表达分类器（也称为molecular signatures，ms）进行不确定性评估。为了实现这一目标，通过排列基因labels以构成模拟对象。 然后，针对相应的亚型分类器对每个模拟对象进行测试，以构建零分布。 因此，可以为每个受试者提供分类置信度评估报告，以帮助医生选择治疗方案。 目前，它只适用于Genefu包中的PAM50实现，但可以很容易地扩展到其他molecular signatures中。

（注：运行到第三步用PAM50 centroids作为验证的object<-filtrate(object, verbose=TRUE)时提示出错：'pam50' is not an exported object from 'namespace:genefu'，宠粉曾老师解决后还写了一篇推文，大家可以移步观看！使用R包的内置数据不能通过两个冒号吗？）

1.使用BioConductor的BreastCancerXXX数据集

为了使用PAM50 MS，用户必须加载BioConductor的BreastCancerXXX数据集，其中XXX代表UPP、NKI、VDX、TRANSBIG、Mainz或UNT。例如，我们可以加载NKI数据库，前提是安装了所需的库，使用以下代码：

# install.packages("pbcmc_1.6.0.tar.gz",repos = NULL, type="source")    #用BiocManager安装失败，进行了本地安装。
#
# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("BiocParallel")      #并行计算
# BiocManager::install("breastCancerNKI")   #数据集
# BiocManager::install("genefu")            #PAM50
# #BiocManager::install("biomaRt")          #可能需要加载的依赖包
# #BiocManager::install("vctrs")            #可能需要加载的依赖包library("biomaRt")
library("pbcmc")
library("BiocParallel")
library("breastCancerNKI")
library("genefu")object<-loadBCDataset(Class=PAM50, libname="nki", verbose=TRUE)
object
# A PAM50 molecular permutation classifier object
# Dimensions:
#             nrow ncol
# exprs      24481  337
# annotation 24481   10
# targets      337   21

该对象是一个PAM50的例子，它包含exprs矩阵，该矩阵带有targets data.frame（目标数据框）中的基因表达值、相关注释和临床数据。另一方面，用户可以借用PAM50的构造函数使用他/她自己的数据创建对象，或者使用 .PAM50(MAList_Object) 函数将 Limma MAList 对象转换为PAM50。在第一种情况下，用户只需要：

a）M gene expression object，基因表达对象M，行为基因、列为样本。
b）annotation data.frame，注释数据框，必须包含：“Probe”、“NCBI.gene.Symbol”和“EntrezGene.ID”这三列。
c）targets data.frame，目标数据框，这是一个可选slot。如果要准备它，行数应该与M中的样本数一样，列数与这些样本可用的临床或实验数据数一样。以“nki”为例，有21列。

2.使用任何微阵列R数据包

可以使用pbcmc直接加载包并将数据提取到PAM50对象所需的expression M、annotation和targets(可选)slot中来构建与上一节相同的示例。为了简单起见，我们将使用前五个样本，它也适用于单个样本。

library("breastCancerNKI")
data("nki")##The expression
M<-exprs(nki)[, 1:5, drop=FALSE]
dim(M)
head(M)
#                 NKI_4  NKI_6  NKI_7  NKI_8  NKI_9
# Contig45645_RC -0.215  0.071  0.182 -0.343 -0.134
# Contig44916_RC -0.207  0.055  0.077  0.302  0.051
# D25272         -0.158 -0.010  0.059  0.169 -0.007
# J00129         -0.819 -0.391 -0.624 -0.528 -0.811
# Contig29982_RC -0.267 -0.310 -0.120 -0.447 -0.536
# Contig26811     0.229  0.157  0.120  0.283 -0.112##The annotation
genes<-fData(nki)[, c("probe", "NCBI.gene.symbol", "EntrezGene.ID")]  #"probe", "NCBI.gene.symbol", "EntrezGene.ID"
head(genes)
#                         probe NCBI.gene.symbol EntrezGene.ID
# Contig45645_RC Contig45645_RC            GREM2         64388
# Contig44916_RC Contig44916_RC            SUHW2        140883
# D25272                 D25272             <NA>            NA
# J00129                 J00129              FGB          2244
# Contig29982_RC Contig29982_RC           SCARA5        286133
# Contig26811       Contig26811             <NA>            NA##Additional information 附加信息(可选的)
targets<-pData(nki)[1:5, ,drop=FALSE]
head(targets)
#       samplename dataset series id filename size age er grade pgr her2
# NKI_4      NKI_4     NKI    NKI  4       NA  2.0  41  1     3  NA   NA
# NKI_6      NKI_6     NKI    NKI  6       NA  1.3  49  1     2  NA   NA
# NKI_7      NKI_7     NKI    NKI  7       NA  2.0  46  0     1  NA   NA
# NKI_8      NKI_8     NKI    NKI  8       NA  2.8  48  0     3  NA   NA
# NKI_9      NKI_9     NKI    NKI  9       NA  1.5  48  1     3  NA   NA
#       brca.mutation e.dmfs t.dmfs node t.rfs e.rfs treatment tissue
# NKI_4             0      0   4747    0  4747     0         0      1
# NKI_6             0      0   4075    0  4075     0         0      1
# NKI_7             0      0   3703    0  3703     0         0      1
# NKI_8             0      0   3215    0  3215     0         0      1
# NKI_9             0      0   3760    0  3760     0         0      1
#       t.os e.os
# NKI_4 4744    0
# NKI_6 4072    0
# NKI_7 3700    0
# NKI_8 3213    0
# NKI_9 3757    0

现在我们准备进行下一节的工作流程，即：

1.通过PAM50对基因进行筛选，只保留所需的50个基因。
2.使用genefu的PAM50算法对样本进行分类，并选择所需的基因归一化方式：none, scale, robust 或者 median。对于单个样本，使用“none”，但我们建议将“median”用于群体（population approaches）。
3.排列基因分类标签以构建零分布，并生成 pbcmc 中提出的不确定性估计。
并用summary、subjectReport和databaseReport进一步探索得到的结果。

3.用PAM50 centroids作为验证

例如，我们可以使用genefu的PAM50 centroids来检查我们的实现，我们事先知道每个样本的真实分类：

library("genefu")
data(pam50)
M<-pam50$centroids  ##The expression
genes<-pam50$centroids.map  ##The annotation
names(genes)<-c("probe", "NCBI.gene.symbol", "EntrezGene.ID")
object<-PAM50(exprs=M, annotation=genes)
object
# A PAM50 molecular permutation classifier object
# Dimensions:
#            nrow ncol
# exprs        50    5
# annotation   50    3
# targets       0    0

请注意，对于上面的输出，targets slot是空的，即nrow=0和nol=0。此外，只包括这50个基因和5个IGP（intrinsic genes profiles）的表达，以及它们对应的三个注释（“probe”,“NCBI.gene.symbol” 和“EntrezGene.ID”）。

建议查看slot（exprs、annotation和targets）以确定它们已正确加载：

head(exprs(object))   ##每个受试者的基因表达
#             Basal       Her2        LumA        LumB     Normal
# ACTR3B  0.7183319 -0.4816657  0.00998107 -0.19055133  0.4657229
# ANLN    0.5373723  0.2669316 -0.57924572  0.09880418 -0.8369396
# BAG1   -0.5745069 -0.4760729  0.75822116 -0.40545862  0.3165530
# BCL2   -0.1187604 -0.1579140  0.28748744 -0.44133950  0.5339789
# BIRC5   0.3004886  0.4057331 -0.88143437  0.60385078 -0.8766364
# BLVRA  -0.6426775  0.3353360  0.04204202  0.69120496 -0.1634128head(annotation(object))  ##必须保留的注释部分
#         probe NCBI.gene.symbol EntrezGene.ID
# ACTR3B ACTR3B           ACTR3B         57180
# ANLN     ANLN             ANLN         54443
# BAG1     BAG1             BAG1           573
# BCL2     BCL2             BCL2           596
# BIRC5   BIRC5            BIRC5           332
# BLVRA   BLVRA            BLVRA           644head(targets(object))  ##临床数据（如果有的话）
#data frame with 0 columns and 0 rows  ##这里为空

正如所预期的那样，五个centroids被加载到Exprs slot中，其对应的“Probe”、“NCBI.gene.Symbol”和“EntrezGene.ID”编号在annotation slot中，并且没有targets的可用数据。现在，可以按照筛选、分类和排列的建议工作流程处理数据：

object<-filtrate(object, verbose=TRUE)   #筛选
object<-classify(object, std="none", verbose=TRUE)   #分类
object<-permutate(object, nPerm=10000, pCutoff=0.01, where="fdr",corCutoff=0.1, keep=TRUE, seed=1234567890, verbose=TRUE,BPPARAM=bpparam())   #排列

筛选的目的是只保留将在分类中起作用的基因。在这个例子中，它不会对原来的exprs slot产生任何改变，只是提取所需的50个基因。但是，如果提供了一个完整的微阵列，那么不是给IGP（intrinsic genes profiles）编码的探针将被移除。此外，相同基因的编码(重复或类似的注释)的探针将按照标准化(STD)参数所述进行处理。

筛选基因过后，就可以使用原始的PAM50算法对它们进行分类。但这时，建议使用至少β=10000次排列来获得亚型分配置信度，对调整后的p值使用I类误差α=0.01，相关性差异阈值CorCutoff=0.1。

事实上，这个过程是计算密集型的，所以可以使用BiocParallel包（BPPARAM=bpparam()）利用所有可用的计算核心。此外，用户可以通过包含verbose=TRUE选项来跟进排列进度。如果我们现在看一下object：

object
# A PAM50 molecular permutation classifier object
# Dimensions:
#            nrow ncol
# exprs        50    5
# annotation   50    3
# targets       0    0
#
# Classification:
#            nrow ncol
# probability   5    5
# correlation   5    5# $subtype
# Basal Her2 LumA LumB Normal
#     1    1    1    1      1# Permutations test ran with following parameters:
#   Permutations=10000, fdr<0.01, corCutoff>0.1, keep=TRUE
# Permutation:
#   correlation available: TRUE
#            nrow ncol
# pvalues       5    5
# fdr           5    5
# subtype       5    5

我们可以看到它已经更新了。首先，classification slot包含两个数据集：一个包含亚型的概率P(IGPi)，以及每个样本与五个IGP（intrinsic genes profiles）的相关性。$subtype项显示了样本对应的IGP频率表。此外，还用pvalues, fdr 和 subtypes显示了所使用的排列参数（permutation parameters）。请注意，在这种情况下，使用了keep=TRUE选项，因此可以使用模拟的相关零分布数据点(simulated correlation null distribution data points ，ρH0IGP)。

在本例中，我们使用了genefu的PAM50 centroids，因此，对于对象输出中的每个IGP只有一个样本。这在pbcmc包确定的分类和原始亚型之间的summary(object)矩阵中也能看出来。此外，这个示例仅显示分配给原始 PAM50亚型的样本(assigned subjects，A)，而未分配 (not assigned，NA) 的行/列仅包含零 (0)。如果发现有不明确的(ambiguous，AMB)样本，则Classes列将包括有争议的亚型的附加行（例如，“LumA, Normal” 或“Her2, LumB”等)。

summary(object)
#            Subtype
# Classes      Basal  Her2 LumA LumB Normal Not Assigned
# Basal            1     0    0    0      0            0
# Her2             0     1    0    0      0            0
# LumA             0     0    1    0      0            0
# LumB             0     0    0    1      0            0
# Normal           0     0    0    0      1            0
# Not Assigned     0     0    0    0      0            0# 最后，我们可以检查单个受试者的报告了解分类的情况，以便为医生的治疗提供适当的建议：
subjectReport(object, subject=1)

下图是一个 grid.arrange 对象，它基本上由三个主要部分组成：

（1）tableGrob：包含以下部分的汇总表

$Summary: 样本名和通过 PAM50.Subtype 或建议的方法 (Permuted.Subtype) 获得亚型。

$Fields ：对于五种 PAM50 亚型：
①相关性：第 i 个 PAM50 centroid与样本表达量的相关性ρ(PGP, IGPi)。
②p 值：使用模拟数据 pIGP 获得的排列p值。
③FDR：使用错误发现率调整的p值。

（2）facet wrap：用ggplot2绘制的样本表达值与PAM50 centroid的散点图和线性回归拟合线(蓝色)。如果样本具有唯一的亚型，则图将显示为红色。此外，如果使用keep=TRUE选项运行模拟排列，则将相应的未排列相关性在零分布箱线图上绘制为大圆点。

（3）textGrob：模拟中使用的排列参数（permutation parameter slot ）。

pbcmc还能使用databaseReport函数获得整个数据库的pdf报告。报告第一页是数据库的总体摘要，即对照原始PAM50亚型结果的置换检验分类的summary统计列联表。如图所示是输出的样本报告。

图例：genefu得到的“Basal”内在基因谱(intrinsic gene profile，IGP)的PAM50排序样本报告。最上面的表格总结了第i个样本的结果，例如有每个IGP的相关性、p值和错误发现率(FDR)。此外，还有样本基因表达值与 IGP 的散点图和线性回归线（蓝色）。红色表示分配的亚型。最后还有每个IGP（intrinsic genes profiles）的零分布箱线图（带有代表未排列相关性的大圆点）。

pbcmc包的介绍（根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理）相关推荐

  1. CancerSubtypes包的介绍（根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理）

    CancerSubtypes包的介绍(根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理,感谢Jimmy老师!) 1. 引言 2. 数据处理 2.1 基本处理 2.1.1 通过检查数据分布来分 ...

  2. ProTICS包的介绍（根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理）

    ProTICS包的介绍(根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理,感谢Jimmy老师!) 1.设置环境 2.Part1的结果 3.Part2的结果 4.Part3的结果 5.相关函数 ...

  3. 生信技能树 电脑配置linux,生信技能树--Jimmy - Linux 20题

    一.在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列. 二.在创建好的文件夹下面,比如我的是 /Users/jimmy/tmp/1/2/3/4/5/6/7/8/9 ,里面 ...

  4. R：生信技能树学习笔记一

    生信技能树小破站:R应该这样学1-4 1.查看已经安装的包的地址 .libPaths() 2.怎么查看函数用法 #在RStudio的右下角窗口的help可以看到 ?函数名 3.三个有用的函数 1.he ...

  5. R：生信技能树学习笔记二

    生信技能树小破站:R应该这样学5-7 1.热图 rm(list=ls()) library(pheatmap) a1=rnorm(100) dim(a1)=c(5,20) #设置维度 pheatmap ...

  6. 生信技能树【代码大全搜录】

    生信技能术代码大全: rm(list = ls()) options()$repos options()$BioC_mirror #options(BioC_mirror="https:// ...

  7. python表情包多样化聊天室_Python | 信不信我分分钟批量做你大堆的表情包？

    作者 | 依然很拉风 作为一个数据分析师,应该信奉一句话----"一图胜千言".不过这里要说的并不是数据可视化,而是一款全民向的产品形态----表情包!!!! 表情包不仅仅是一种符 ...

  8. 生信技能树linux虚拟机,科学网—Windows10安装Linux子系统Ubuntu 20.04LTS，轻松使用生信软件，效率秒杀虚拟机 - 刘永鑫的博文...

    很多优秀的生物信息学软件,如QIIME.QIIME 2.LEfSe等没有Windows版,而使用VirutalBox虚拟机不仅效率低,而且挂载外部硬盘和使用中也经常遇到各种问题,配置和使用详见 - 扩 ...

  9. 生信技能树课程记录笔记（七）20220531

    一.数据框排序 法一:sort函数 默认升序. 例:sort(test$Sepal.Length) 法二:order函数 默认升序,返回数据下标组成的数组. 可以给向量排序,也可以给数据框排序 例:t ...

  10. TCGA学习笔记一(生信技能树概述版)

    1.背景介绍 重要数据 外显子数据 表达数据 小RNA测序数据 拷贝数芯片 甲基化数据 蛋白质组学数据 临床信息 癌症背景知识 网页工具大全 GDC cbioportal:按照paper来分类的 UC ...

最新文章

  1. Java 的zip压缩和解压缩
  2. 运维工作钱少、事多而且杂？年轻人，你这个思想很危险吶
  3. 不认识java代码_程序员进阶：优雅的代码对于一个架构师的重要性
  4. python列表切片规则_Python 列表切片
  5. Flash 插件又被曝出新漏洞，让攻击者可以控制 Mac
  6. linux中vim中文显示乱码
  7. Kotlin 的优点
  8. 你真的懂JavaScript基础类型吗
  9. html5查看ies文件,5千+ IES光域网文件 5312 IES Files + IES 预览
  10. 网关串口+EM310
  11. get请求中文乱码问题解决
  12. lcx端口转发及远程终端问题
  13. Android MVP 实践之路（理解篇）
  14. Android之提示Could not find com.android.support:appcompat-v7:25.3.1
  15. 灰色预测的MATLAB程序
  16. 上升了百分之几怎么算_上涨百分之多少怎么算
  17. 国医大师王绵之：汤药煎服经验谈
  18. Web开发技术的演变
  19. 剑指 Offer 46. 把数字翻译成字符串（javascript）
  20. seo教程电子书（SEO搜索引擎优化基础教程）

热门文章

  1. PyQt5最全27 绘图之drawLine绘制不同类型的直线
  2. 关于Win10系统下VIA HD AUDIO威盛声卡没声音问题 - 有效解决办法
  3. 电脑没有声音怎么安装声卡驱动?驱动人生声卡驱动安装失败原因
  4. LiquidCrystal_I2C 显示不正常 只显示第一个首字符！
  5. 用java编写猜数字游戏
  6. font: 0.5rem/1 tahoma, arial, 'Microsoft YaHei', simsun;
  7. 圣经闪卡 - Holy Bible Flash Cards
  8. 2019汤家凤高等数学强化班讲义
  9. java web 和js区别_jsp和javascript之间有什么区别？
  10. 计算机打字键盘亮怎么设置,键盘指示灯亮着却不能打字的解决方法