劉小澤寫于19.7.8-第二單元第八講：乳腺癌領(lǐng)域之PAM50分類
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
什么是PAM50
首次接觸這個(gè)名詞肯定很蒙，因?yàn)樗侨橄侔╊I(lǐng)域的分類名詞，需要看許多資料才能了解，我也一樣，看了一些推文、英文資料、文章，才做了一些總結(jié)
PAM50的意思是Prediction Analysis of Microarray 50 ，可以看到是芯片時(shí)代的產(chǎn)物了，它是2009年由Parker提出的，原文在：https://ascopubs.org/doi/full/10.1200/JCO.2008.18.1370，目前接近3000引用量。
使用的芯片是Agilent human 1Av2 microarrays or custom-designed Agilent human 22k arrays，數(shù)據(jù)在GSE10886，它研究了189個(gè)prototype samples，得到了一個(gè)50個(gè)分類基因與5個(gè)對(duì)照基因的RT-qPCR定量結(jié)果，得到4個(gè)gene expression–based “intrinsic” subtypes：Luminal A, Luminal B, HER2-enriched and Basal-like(詳見：https://genome.unc.edu/pubsup/breastGEO/pam50_centroids.txt)。
關(guān)于這幾種分子亞型的介紹：https://www.breastcancer.org/symptoms/types/molecular-subtypes
Luminal A：hormone-receptor positive (estrogen-receptor and/or progesterone-receptor positive), HER2 negative, low levels of the protein Ki-67 => grow slowly and have the best prognosis.
Luminal B：hormone-receptor positive (estrogen-receptor and/or progesterone-receptor positive), either HER2 positive or HER2 negative，high levels of Ki-67 => grow slightly faster than luminal A & prognosis is slightly worse
-
Triple-negative/basal-like: hormone-receptor negative (estrogen-receptor and progesterone-receptor negative) , HER2 negative
More common with BRCA1 gene mutations among younger and African-American women..
-
HER2-enriched: hormone-receptor negative (estrogen-receptor and progesterone-receptor negative), HER2 positive => grow faster than luminal cancers & worse prognosis
BUT often successfully treated with targeted therapies aimed at the HER2 protein (e.g. Herceptin, Perjeta, Tykerb, Nerlynx, Kadcyla)
Normal-like: similar to luminal A => prognosis is slightly worse than luminal A but also good
乳腺癌發(fā)育來(lái)自兩種不同的細(xì)胞：基體細(xì)胞和管腔細(xì)胞，還有不同的激素類型(ER/PR、HER2受體)，之前在臨床上都是根據(jù)一些IHC marker來(lái)進(jìn)行分類
The most common immunohistochemical breast cancerprognostic and therapeutic markers used include: estrogen receptor, human epidermal growth factor receptor-2, Ki-67, progesterone receptor, and p53. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4127609/)
乳腺癌是一個(gè)高度異質(zhì)性的疾病，即使臨床分期和病理分級(jí)相同，患者對(duì)治療的反應(yīng)和預(yù)后也是不同的。目前仍然是根據(jù)臨床病理特點(diǎn)如HER2表達(dá)、雌激素受體狀態(tài)、腫瘤大小、分級(jí)和淋巴結(jié)轉(zhuǎn)移等選擇輔助治療，包括化療，內(nèi)分泌治療，抗HER2治療等。為了指導(dǎo)預(yù)后，常常采用TNM分期、臨床病理指標(biāo)，后來(lái)由于高通量數(shù)據(jù)的產(chǎn)生，多基因預(yù)測(cè)成為了一個(gè)新途徑。
舉個(gè)例子：可以看表達(dá)量，比如有50個(gè)基因，有10個(gè)特定基因高它們就表示Luminal A，有其他10個(gè)基因高就是Luminal B，這就是一個(gè)模式；我們只需要比較我們的表達(dá)矩陣和這個(gè)模式進(jìn)行對(duì)應(yīng)
多基因檢測(cè)有兩項(xiàng)已經(jīng)通過(guò)了FDA的批準(zhǔn)：
- 21-gene OncotypeDx assay (Genome Health Inc, Redwood City, CA)：risk stratify early-stage estrogen receptor (ER) –positive breast cancer
- 70-gene MammaPrint (Agendia, Huntington Beach, CA)：ER-positive and ER-negative early-stage node-negative breast cancer.
另外前人的研究還有：
- Single Sample Predictor (SSP) ：SSP2003 、SSP2006、PAM50
- Subtype Classification Model (SCM)：SCMOD1、SCMOD2 、simple three-gene model (SCMGENE）
利用genefu包來(lái)熟悉PAM50分類器
這個(gè)是Bioconductor的包，使用正確的方式安裝好
官方教程在：https://rdrr.io/bioc/genefu/f/inst/doc/genefu.pdf
用包需知 1
自帶了5個(gè)乳腺癌芯片數(shù)據(jù)集（breastCancerMAINZ=》GSE11121、breastCancerTRANSBIG=》GSE7390、breastCancerUPP=》GSE3494、breastCancerUNT=》GSE2990、breastCancerNKI=》數(shù)據(jù)集沒(méi)有上傳到GEO）：https://vip.biotrainee.com/d/689-5
breastCancerMAINZ=》GSE11121
文章：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://sci-hub.tw/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
breastCancerTRANSBIG=》GSE7390
文章：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://sci-hub.tw/10.1158/1078-0432.ccr-06-2765
方法： GPL96 （HG-U133A） Affymetrix Human Genome U133A Array 芯片，frozen samples from 198 N- systemically untreated patients
breastCancerUPP=》GSE3494
文章：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://sci-hub.tw/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.
breastCancerUNT =》GSE2990
文章：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://sci-hub.tw/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.
最后一個(gè)breastCancerNKI使用的是Agilent公司芯片
用包需知 2
這個(gè)R包除了包裝了PAM50分類，還加入了其他許多分類標(biāo)準(zhǔn)，詳見https://rdrr.io/bioc/genefu/man/，使用PAM50是因?yàn)樗囊昧亢芨撸J(rèn)可度較高
開始用包
# 加載數(shù)據(jù)
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"
除了很多不像常規(guī)基因名的基因以外，還有很多基因大小寫不一致，這是因?yàn)檫@個(gè)數(shù)據(jù)是小鼠的，而小鼠的基因名與人類的不同在于：首字母大寫，其余小寫
首先就是將這里的dat基因名全變?yōu)榇髮?/p>
rownames(dat)=toupper(rownames(dat))
當(dāng)然，最好直接使用小鼠的分類器，但是目前沒(méi)有，因此只能使用人類的，不是很準(zhǔn)確，但是這個(gè)分類是可以借鑒的
# 加載genefu
library(genefu)
# 可以看到會(huì)加載很多依賴包，包含機(jī)器學(xué)習(xí)、并行、分類法
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
這個(gè)包也需要轉(zhuǎn)置后的表達(dá)矩陣(基因?yàn)榱?
> 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"
## 發(fā)現(xiàn)有的基因名是不符合常規(guī)認(rèn)知的，因此需要進(jìn)行基因名轉(zhuǎn)換
# 看下人類這個(gè)基因注釋包中都包含哪些，發(fā)現(xiàn)有org.Hs.egSYMBOL，應(yīng)該就是需要的
ls("package:org.Hs.eg.db")
# 這個(gè)注釋信息是Bimap格式的，需要先轉(zhuǎn)換成數(shù)據(jù)框，利用toTable函數(shù)
> class(org.Hs.egSYMBOL)
[1] "AnnDbBimap"
> s2g=toTable(org.Hs.egSYMBOL)
# 求小鼠的基因與人類的基因的交集，利用match函數(shù)，返回位置信息(如果沒(méi)有對(duì)應(yīng)，就返回NA)。存在NA的原因就是：小鼠有的對(duì)應(yīng)不上人類基因名，并且人類的基因也有未知的
> g=s2g[match(s,s2g$symbol),1]
# 然后做成一個(gè)數(shù)據(jù)框
> dannot=data.frame(probe=s,
"Gene.Symbol" =s,
"EntrezGene.ID"=g)
# 下面去掉ddata和dannot中NA的行
> ddata=ddata[,!is.na(dannot$EntrezGene.ID)] #ID轉(zhuǎn)換
> dim(ddata)
[1] 768 10487 # 相比之前大約去掉2000個(gè)基因
> dannot=dannot[!is.na(dannot$EntrezGene.ID),]
# 看下去除NA后的基因注釋和表達(dá)矩陣，必須保證注釋的基因ID和表達(dá)矩陣的基因ID一一對(duì)應(yīng)
> 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
可以進(jìn)行g(shù)enefu分析了，分型就是使用molecular.subtyping函數(shù)
s<-molecular.subtyping(sbt.model = "pam50",data=ddata,
annot=dannot,do.mapping=TRUE)
# 結(jié)果就是將768個(gè)細(xì)胞
> table(s$subtype)
Basal Her2 LumB LumA Normal
42 58 46 543 79
# 可以利用原始的樣本信息數(shù)據(jù)框df進(jìn)行clust分組與分子分型之間關(guān)系的探索
> 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
注意：雖然這里可以實(shí)現(xiàn)分類，但是PAM50是針對(duì)乳腺癌患者進(jìn)行分類的，而我們這里是針對(duì)單細(xì)胞；而且細(xì)胞也不是癌細(xì)胞，是CAFs(cancer associated fiberblast）
不管是什么細(xì)胞，最后都能得到一個(gè)表達(dá)矩陣，算法是不會(huì)考慮矩陣來(lái)源的，因此即便是正常細(xì)胞的矩陣，也可以分類成5種乳腺癌亞型，所以分類的前提還是自己熟悉數(shù)據(jù)的生物學(xué)背景
探索PAM50
看一下pam50，它是一個(gè)列表
> 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 ...
然后取出基因名，存儲(chǔ)在centroids中：
pam50genes=pam50$centroids.map[c(1,3)]
# 發(fā)現(xiàn)有的基因已經(jīng)不是標(biāo)準(zhǔn)的symbol了，PAM50是2009年的基因名，因此需要進(jìn)行修改
pam50genes[pam50genes$probe=='CDCA1',1]='NUF2'
pam50genes[pam50genes$probe=='KNTC2',1]='NDC80'
pam50genes[pam50genes$probe=='ORC6L',1]='ORC6'
以第一個(gè)基因?yàn)槔?a target="_blank">https://www.genecards.org/cgi-bin/carddisp.pl?gene=NUF2&keywords=NUF2
> x=dat
# 找到pam50在原始表達(dá)矩陣行名中的基因，發(fā)現(xiàn)一共有38個(gè)
> 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)] ,]
下面進(jìn)行熱圖可視化
# 在原來(lái)group_list基礎(chǔ)上，添加亞型信息，為了下面pheatmap中的anno_col設(shè)置
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)
圖片本身不重要，因?yàn)檫@里數(shù)據(jù)的使用是不合適的。可以看到，大部分基因都是luminal A
如果要繼續(xù)歸一化就是：
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)
