前 · 言 第二单元第八讲:乳腺癌领域之PAM50分类 什么是PAM50 首次接触这个名词肯定很蒙,因为它是乳腺癌领域的分类名词,需要看许多资料才能了解,我也一样,看了一些推文、英文资料、文章,才做了一些总结 PAM50的意思是 使用的芯片是Agilent human 1Av2 microarrays or custom-designed Agilent human 22k arrays,数据在GSE10886,它研究了189个prototype samples,得到了一个50个分类基因与5个对照基因的RT-qPCR定量结果,得到4个gene expression–based “intrinsic” subtypes:Luminal A, Luminal B, HER2-enriched and Basal-like(详见:https://genome./pubsup/breastGEO/pam50_centroids.txt)。 关于这几种分子亚型的介绍:https://www./symptoms/types/molecular-subtypes
乳腺癌发育来自两种不同的细胞:基体细胞和管腔细胞,还有不同的激素类型(ER/PR、HER2受体),之前在临床上都是根据一些IHC marker来进行分类
乳腺癌是一个高度异质性的疾病,即使临床分期和病理分级相同,患者对治疗的反应和预后也是不同的。目前仍然是根据临床病理特点如HER2表达、雌激素受体状态、肿瘤大小、分级和淋巴结转移等选择辅助治疗,包括化疗,内分泌治疗,抗HER2治疗等。为了指导预后,常常采用TNM分期、临床病理指标,后来由于高通量数据的产生,多基因预测成为了一个新途径。
多基因检测有两项已经通过了FDA的批准:
另外前人的研究还有:
利用genefu包来熟悉PAM50分类器
用包需知 1自带了5个乳腺癌芯片数据集(breastCancerMAINZ=》GSE11121、breastCancerTRANSBIG=》GSE7390、breastCancerUPP=》GSE3494、breastCancerUNT=》GSE2990、breastCancerNKI=》数据集没有上传到GEO):https://vip./d/689-5
文章:The humoral immune system has a key prognostic impact in node-negative breast cancer. Cancer Res 2008 Jul 1;68(13):5405-13. Sci-hub: https:///10.1158/0008-5472.can-07-5206 方法:GPL96(HG-U133A) Affymetrix Human Genome U133A Array芯片,其中包含了200 tumors of patients who were not treated by systemic therapy after surgery using a discovery approach. 临床信息:biological process of proliferation、steroid hormone receptor expression、B cell and T cell infiltration
文章:Strong time dependence of the 76-gene prognostic signature for node-negative breast cancer patients in the TRANSBIG multicenter independent validation series. Clin Cancer Res 2007 Jun 1;13(11):3207-14. Sci-hub: https:///10.1158/1078-0432.ccr-06-2765 方法:GPL96 (HG-U133A) Affymetrix Human Genome U133A Array 芯片,frozen samples from 198 N- systemically untreated patients
文章:An expression signature for p53 status in human breast cancer predicts mutation status, transcriptional effects, and patient survival. Proc Natl Acad Sci U S A 2005 Sep 20;102(38):13550-5. Sci-hub:https:///10.2307/3376671 方法:GPL96 (HG-U133A) Affymetrix Human Genome U133A Array 芯片,freshly frozen breast tumors from a population-based cohort of 315 women representing 65% of all breast cancers resected in Uppsala County, Sweden, from January 1, 1987 to December 31, 1989.
文章:Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst 2006 Feb 15;98(4):262-72 Sci-hub:https:///10.1093/jnci/djj052 方法:GPL96 (HG-U133A) Affymetrix Human Genome U133A Array 芯片, 189 invasive breast carcinomas and from three published gene expression datasets from breast carcinomas.
用包需知 2这个R包除了包装了PAM50分类,还加入了其他许多分类标准,详见https:///bioc/genefu/man/,使用PAM50是因为它的引用量很高,认可度较高 开始用包# 加载数据rm(list = ls()) options(stringsAsFactors = F) load(file = '../input.Rdata') a[1:4,1:4] head(df) # 检查行名(基因名) > head(rownames(dat)) [1] "0610007P14Rik" "0610009B22Rik" "0610009L18Rik" "0610009O20Rik" [5] "0610010F05Rik" "0610010K14Rik" 除了很多不像常规基因名的基因以外,还有很多基因大小写不一致,这是因为这个数据是小鼠的,而小鼠的基因名与人类的不同在于:首字母大写,其余小写 首先就是将这里的 当然,最好直接使用小鼠的分类器,但是目前没有,因此只能使用人类的,不是很准确,但是这个分类是可以借鉴的 # 加载genefulibrary(genefu) # 可以看到会加载很多依赖包,包含机器学习、并行、分类法 Loading required package: limma Loading required package: biomaRt Loading required package: iC10 Loading required package: pamr Loading required package: cluster Loading required package: impute Loading required package: iC10TrainingData Loading required package: AIMS Loading required package: e1071 Loading required package: Biobase Loading required package: BiocGenerics Loading required package: parallel 这个包也需要转置后的表达矩阵(基因为列) > ddata=t(dat)> ddata[1:4,1:4] 0610007P14Rik 0610009B22Rik 0610009L18Rik 0610009O20Rik SS2_15_0048_A3 0.000000 0 0 0.000000 SS2_15_0048_A6 0.000000 0 0 0.000000 SS2_15_0048_A5 6.459884 0 0 2.544699 SS2_15_0048_A4 6.313884 0 0 3.025273 > dim(ddata) [1] 768 12198 > s=colnames(ddata);head(s);tail(s) ##获得基因名 [1] "0610007P14Rik" "0610009B22Rik" "0610009L18Rik" "0610009O20Rik" [5] "0610010F05Rik" "0610010K14Rik" [1] "ERCC-00160" "ERCC-00162" "ERCC-00163" "ERCC-00165" "ERCC-00170" [6] "ERCC-00171" ## 发现有的基因名是不符合常规认知的,因此需要进行基因名转换 # 看下人类这个基因注释包中都包含哪些,发现有org.Hs.egSYMBOL,应该就是需要的 ls("package:org.Hs.eg.db") # 这个注释信息是Bimap格式的,需要先转换成数据框,利用toTable函数 > class(org.Hs.egSYMBOL) [1] "AnnDbBimap" > s2g=toTable(org.Hs.egSYMBOL) # 求小鼠的基因与人类的基因的交集,利用match函数,返回位置信息(如果没有对应,就返回NA)。存在NA的原因就是:小鼠有的对应不上人类基因名,并且人类的基因也有未知的 > g=s2g[match(s,s2g$symbol),1] # 然后做成一个数据框 > dannot=data.frame(probe=s, "Gene.Symbol" =s, "EntrezGene.ID"=g) # 下面去掉ddata和dannot中NA的行 > ddata=ddata[,!is.na(dannot$EntrezGene.ID)] #ID转换 > dim(ddata) [1] 768 10487 # 相比之前大约去掉2000个基因 > dannot=dannot[!is.na(dannot$EntrezGene.ID),] # 看下去除NA后的基因注释和表达矩阵,必须保证注释的基因ID和表达矩阵的基因ID一一对应 > head(dannot) probe Gene.Symbol EntrezGene.ID 372 A4GALT A4GALT 53947 393 AAAS AAAS 8086 394 AACS AACS 65985 396 AAGAB AAGAB 79719 397 AAK1 AAK1 22848 398 AAMDC AAMDC 28971 > ddata[1:4,1:4] A4GALT AAAS AACS AAGAB SS2_15_0048_A3 8.516383 0 0.000000 0 SS2_15_0048_A6 7.111928 0 0.000000 0 SS2_15_0048_A5 3.415452 0 0.000000 0 SS2_15_0048_A4 6.848774 0 7.168196 0 可以进行genefu分析了,分型就是使用 annot=dannot,do.mapping=TRUE) # 结果就是将768个细胞 > table(s$subtype) Basal Her2 LumB LumA Normal 42 58 46 543 79 # 可以利用原始的样本信息数据框df进行clust分组与分子分型之间关系的探索 > df$subtypes=subtypes > table(df[,c(1,5)]) subtypes g Basal Her2 LumA LumB Normal 1 36 30 205 13 28 2 3 25 217 31 24 3 1 2 102 1 15 4 2 1 19 1 12
探索PAM50看一下pam50,它是一个列表 > str(pam50)List of 7 $ method.cor : chr "spearman" $ method.centroids: chr "mean" $ std : chr "none" $ rescale.q : num 0.05 $ mins : num 5 $ centroids : num [1:50, 1:5] 0.718 0.537 -0.575 -0.119 0.3 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:50] "ACTR3B" "ANLN" "BAG1" "BCL2" ... .. ..$ : chr [1:5] "Basal" "Her2" "LumA" "LumB" ... $ centroids.map :'data.frame':50 obs. of 3 variables: ..$ probe : chr [1:50] "ACTR3B" "ANLN" "BAG1" "BCL2" ... ..$ probe.centroids: chr [1:50] "ACTR3B" "ANLN" "BAG1" "BCL2" ... ..$ EntrezGene.ID : int [1:50] 57180 54443 573 596 332 644 891 898 991 990 ... 然后取出基因名,存储在 # 发现有的基因已经不是标准的symbol了,PAM50是2009年的基因名,因此需要进行修改 pam50genes[pam50genes$probe=='CDCA1',1]='NUF2' pam50genes[pam50genes$probe=='KNTC2',1]='NDC80' pam50genes[pam50genes$probe=='ORC6L',1]='ORC6' 以第一个基因为例:https://www./cgi-bin/carddisp.pl?gene=NUF2&keywords=NUF2 > x=dat # 找到pam50在原始表达矩阵行名中的基因,发现一共有38个 > pam50genes$probe[pam50genes$probe %in% rownames(x)] [1] "ANLN" "BAG1" "BCL2" "BIRC5" "BLVRA" "CCNB1" [7] "CCNE1" "CDC20" "CDC6" "NUF2" "CDH3" "CENPF" [13] "CEP55" "CXXC5" "EGFR" "ERBB2" "ESR1" "FOXC1" [19] "KIF2C" "NDC80" "MAPT" "MDM2" "MELK" "MIA" [25] "MKI67" "MLPH" "MMP11" "MYBL2" "MYC" "ORC6" [31] "PHGDH" "PTTG1" "RRM2" "SFRP1" "SLC39A6" "TYMS" [37] "UBE2C" "UBE2T" > x=x[pam50genes$probe[pam50genes$probe %in% rownames(x)] ,] 下面进行热图可视化 # 在原来group_list基础上,添加亚型信息,为了下面pheatmap中的anno_col设置tmp=data.frame(group=group_list, subtypes=subtypes) rownames(tmp)=colnames(x) # 画热图 library(pheatmap) pheatmap(x,show_rownames = T,show_colnames = F, annotation_col = tmp)
如果要继续归一化就是: x=t(scale(t(x)))x[x>1.6]=1.6 x[x< -1.6]= -1.6 pheatmap(x,show_rownames = T,show_colnames = F, annotation_col = tmp) |
|