乳腺癌领域之PAM50分类

乳腺癌领域之PAM50分类
第⼆单元第⼋讲：乳腺癌领域之P A M50分类
⾸次接触这个名词肯定很蒙，因为它是乳腺癌领域的分类名词，需要看许多资料才能了解，我也⼀样，看了⼀些推⽂、英⽂资料、⽂章，才做了⼀些总结P A M 50的意思是P r e d i c t i o n A n a l y s i s o f M i c r o a r r a y 50 ，可以看到是芯⽚时代的产物了，它是2009年由P ar k er 提出的，原⽂在：ht t ps ://as c opubs.or g/doi/f ull/10.1200/J C O.2008.18.1370，⽬前接近3000引⽤量。

使⽤的芯⽚是A gilent hum an 1A v 2 m ic r oar r ay s or c us t om -des igned A gilent hum an 22k ar r ay s ，数据在G S E 10886，它研究了189个pr ot ot y pe s am ples ，得到了⼀个50个分类基因与5个对照基因的R T -qP C R 定量结果，得到4个gene ex pr es s ion–bas ed “int r ins ic ” s ubt y pes ：Lum inal A , Lum inal B , H E R 2-
enr ic hed and B as al-lik e(详见：ht t ps ://genom /pubs up/br eas t G E O /p am 50_c ent r oids.t x t )。

关于这⼏种分⼦亚型的介绍：ht t ps ://w w w.br eas t c anc er.or g/s y m pt om s /t y pes /m olec ular -s ubt y pes
前 · ⾔
什么是PAM50
L u mi n al A：hor m one-r ec ept or p o si t i ve(es t r ogen-r ec ept or and/or
pr oges t er one-r ec ept or pos it iv e),H E R2 n eg at i ve,l o w lev els of t he pr ot ein K i-67 =>gr ow s low ly and hav e t he bes t pr ognos is.
L u mi n al B：hor m one-r ec ept or p o si t i ve(es t r ogen-r ec ept or and/or
pr oges t er one-r ec ept or pos it iv e),eit her H E R2 p o si t i ve o r H E R2
n eg at i ve，h i g h lev els of K i-67 =>gr ow s light ly f as t er t han lum inal A&
pr ognos is is s light ly w or s e
T r i p l e-n eg at i ve/b asal-l i ke:hor m one-r ec ept or n eg at i ve(es t r ogen-r ec ept or and pr oges t er one-r ec ept or negat iv e),H E R2 n eg at i ve
Mo re c o mmo n w i th B R C A1 g e n e mu ta ti o n s a mo n g y o u n g e r
a n d A fri c a n-A me ri c a n w o me n..
H E R2-en r i ch ed:hor m one-r ec ept or n eg at i ve(es t r ogen-r ec ept or and
pr oges t er one-r ec ept or negat iv e),H E R2 p o si t i ve=>gr ow f as t er t han
lum inal c anc er s&w or s e pr ognos is
B U T o fte n s u c c e s s fu l l y tre a te d w i th ta rg e te d th e ra p ie s
a i me d a t th e H E R2 p ro te i n (e.g. H e rc e p ti n, P e rj e ta,
T y k e rb, N e rl y n x, K a d c y l a)
N o r mal-l i ke:si mi l ar t o lum inal A=>pr ognos is is s light ly w or s e t han
lum inal A but als o good
乳腺癌发育来⾃两种不同的细胞：基体细胞和管腔细胞，还有不同的激素类型
(E R/P R、H E R2受体)，之前在临床上都是根据⼀些I H C m ar k er来进⾏分类
T he m os t c om m on i mmu n o h i st o ch emi cal b r east
can cer pr ognos t ic and t her apeut ic mar ker s us ed inc lude:
es t r ogen r ec ept or,hum an epider m al gr ow t h f ac t or r ec ept or-2,
K i-67,pr oges t er one r ec ept or,and p53.(ht t ps://w w w.nc bi.nlm.
/pm c/ar t ic les/P M C4127609/)
乳腺癌是⼀个⾼度异质性的疾病，即使临床分期和病理分级相同，患者对治疗的反应和预后也是不同的。

⽬前仍然是根据临床病理特点如H E R2表达、雌激素受体状态、肿瘤⼤⼩、分级和淋巴结转移等选择辅助治疗，包括化疗，内分泌治疗，抗H E R2治疗等。

为了指导预后，常常采⽤T N M分期、临床病理指标，后来由于⾼通量数据的产⽣，多基因预测成为了⼀个新途径。

举个例⼦：可以看表达量，⽐如有50个基因，有10个特定基因⾼它们
就表⽰Lum inal A，有其他10个基因⾼就是Lum inal B，这就是⼀个模
式；我们只需要⽐较我们的表达矩阵和这个模式进⾏对应
多基因检测有两项已经通过了F D A的批准：
21-g en e O n co t yp eD x assay(G enom e H ealt h I nc,R edw ood C it y,C A)：r is k s t r at if y ear ly-s t age es t r ogen r ec ept or(E R)–pos it iv e br eas t c anc er
70-g en e M ammaP r i n t(A gendia,H unt ingt on B eac h,C A)：E R-pos it iv e and E R-negat iv e ear ly-s t age node-negat iv e br eas t c anc er.
另外前⼈的研究还有：
S i n g l e S amp l e P r ed i ct o r(S S P)：S S P2003 、S S P2006、P A M50
S u b t yp e C l assi f i cat i o n M o d el(S C M)：S C M O D1、S C M O D2 、s im ple
t hr ee-gene m odel (S C M G E N E）
利⽤genefu包来熟悉PAM50分类器
这个是B i o c o n d u c to r的包，使⽤正确的⽅式安装好官⽅教程在：h
ttp s://rd rr.i o/b i o c/g e n e fu/f/i n s t/d o c/g e n e fu.p d f
1
⾃带了5个乳腺癌芯⽚数据集（br eas t C anc er M A I N Z=》G S E11121、
br eas t C anc er T R A N S B I G=》G S E7390、br eas t C anc er U P P=》G S E3494、
br eas t C anc er U N T=》G S E2990、br eas t C anc er N K I=》数据集没有上传到
G E O）：ht t ps://v ip.biot r ainee.c om/d/689-5
b re a s tC a n
c e rMA IN Z=》GS E11121
⽂章：T he hum or al im m une s y s t em has a k ey pr ognos t ic im pac t in node-negat iv e br eas t c anc er.C anc er R es2008 J ul 1;68(13):5405-13.
S c i-hub:ht t ps://s c i-hub.t w/10.1158/0008-5472.c an-07-5206
⽅法：G P L96（H G-U133A）A f f y m et r ix H um an G enom e U133A A r r ay芯⽚，其中包含了200 t u mo r s o f p at i en t s w ho w er e not t r eat ed by s y s t em ic
t her apy af t er s ur ger y us ing a dis c ov er y appr oac h.
临床信息：biologic al pr oc es s of pr olif er at ion、s t er oid hor m one r ec ept or
ex pr es s ion、B c ell and T c ell inf ilt r at ion
b re a s tC a n
c e rT R A N S B IG=》GS E7390
⽂章：S t r ong t im e dependenc e of t he 76-g en e p r o g n o st i c si g n at u r e f or node-negat iv e br eas t c anc er pat ient s in t he T R A N S B I G m ult ic ent er independent v alidat ion s er ies.C lin C anc er R es2007 J un 1;13(11):3207-14.
S c i-hub:ht t ps://s c i-hub.t w/10.1158/1078-0432.c c r-06-2765
⽅法：G P L96（H G-U133A）A f f y m et r ix H um an G enom e U133A A r r ay芯⽚，f r oz en s am ples f r om198 N-syst emi cal l y u n t r eat ed p at i en t s
b re a s tC a n
c e rU P P=》GS E3494
⽂章：A n ex pr es s ion s ignat ur e f or p53 s t at us in hum an br eas t c anc er
pr edic t s m ut at ion s t at us,t r ans c r ipt ional ef f ec t s,and pat ient s ur v iv al.P r oc N at l A c ad S c i U S A2005 S ep 20;102(38):13550-5.
S c i-hub：ht t ps://s c i-hub.t w/10.2307/3376671
⽅法：G P L96（H G-U133A）A f f y m et r ix H um an G enom e U133A A r r ay芯⽚，f r es hly f r oz en br eas t t um or s f r om a populat ion-bas ed c ohor t of315
w o men r epr es ent ing 65% of all br eas t c anc er s r es ec t ed in U pps ala C ount y, S w eden,f r om J anuar y1,1987 t o D ec em ber31,1989.
b re a s tC a n
c e rU N T =》GS E2990
⽂章：G ene ex pr es s ion pr of iling in br eas t c anc er:under s t anding t he
m olec ular bas is of his t ologic gr ade t o im pr ov e pr ognos is.J N at l C anc er
I ns t2006 F eb 15;98(4):262-72
S c i-hub：ht t ps://s c i-hub.t w/10.1093/jnc i/djj052
⽅法：G P L96（H G-U133A）A f f y m et r ix H um an G enom e U133A A r r ay芯⽚，189inv as iv e br eas t c ar c inom as and f r om t h r ee p u b l i sh ed gene
ex pr es s ion dat as et s f r om br eas t c ar c inom as.
最后⼀个b re a s tC a n c e rN K I使⽤的是A g i l e n t公司芯⽚
2
这个R包除了包装了P A M50分类，还加⼊了其他许多分类标准，详见ht t ps://r dr r.io/ bioc/genef u/m an/，使⽤P A M50是因为它的引⽤量很⾼，认可度较⾼
#加载数据
r m(l i s t = l s())
o p t i o n s(s t r i n g s A s F a c t o r s = F)
l o a d(f i l e = '../i n p u t.R d a t a')
a[1:4,1:4]
h e a d(d f)
#检查⾏名(基因名)
> h e a d(r o w n a m e s(d a t))
[1] "0610007P14R i k" "0610009B22R i k" "0610009L18R i k" "0610009O20R i k" [5] "0610010F05R i k" "0610010K14R i k"
除了很多不像常规基因名的基因以外，还有很多基因⼤⼩写不⼀致，这是因为这个数据是⼩⿏的，⽽⼩⿏的基因名与⼈类的不同在于：⾸字母⼤写，其余⼩写
⾸先就是将这⾥的d a t基因名全变为⼤写
r o w n a m e s(d a t)=t o u p p e r(r o w n a m e s(d a t))
当然，最好直接使⽤⼩⿏的分类器，但是⽬前没有，因此只能使⽤⼈类的，不是很准确，但是这个分类是可以借鉴的
#加载g e n e f u
l i b r a r y(g e n e f u)
#可以看到会加载很多依赖包，包含机器学习、并⾏、分类法
L o a d i n g r e q u i r e d p a c k a g e: l i m m a
L o a d i n g r e q u i r e d p a c k a g e: b i o m a R t
L o a d i n g r e q u i r e d p a c k a g e: i C10
L o a d i n g r e q u i r e d p a c k a g e: p a m r
L o a d i n g r e q u i r e d p a c k a g e: c l u s t e r
L o a d i n g r e q u i r e d p a c k a g e: i m p u t e
L o a d i n g r e q u i r e d p a c k a g e: i C10T r a i n i n g D a t a
L o a d i n g r e q u i r e d p a c k a g e: A I M S
L o a d i n g r e q u i r e d p a c k a g e: e1071
L o a d i n g r e q u i r e d p a c k a g e: B i o b a s e
L o a d i n g r e q u i r e d p a c k a g e: B i o c G e n e r i c s
L o a d i n g r e q u i r e d p a c k a g e: p a r a l l e l
这个包也需要转置后的表达矩阵(基因为列)
>d d a t a=t(d a t)
>d d a t a[1:4,1:4]
0610007P14R i k 0610009B22R i k 0610009L18R i k 0610009O20R i k
S S2_15_0048_A30.000000000.000000
S S2_15_0048_A60.000000000.000000
S S2_15_0048_A5 6.45988400 2.544699
S S2_15_0048_A4 6.31388400 3.025273
> d i m(d d a t a)
[1]768 12198
>s=c o l n a m e s(d d a t a);h e a d(s);t a i l(s) ##获得基因名
[1] "0610007P14R i k" "0610009B22R i k" "0610009L18R i k" "0610009O20R i k" [5] "0610010F05R i k" "0610010K14R i k"
[1] "E R C C-00160" "E R C C-00162" "E R C C-00163" "E R C C-00165" "E R C C-00170"
[6] "E R C C-00171"
##发现有的基因名是不符合常规认知的，因此需要进⾏基因名转换
#看下⼈类这个基因注释包中都包含哪些，发现有o r g.H s.e g S Y M B O L，应该就是需要的l s("p a c k a g e:o r g.H s.e g.d b")
#这个注释信息是B i m a p格式的，需要先转换成数据框，利⽤t o T a b l e函数
> c l a s s(o r g.H s.e g S Y M B O L)
[1] "A n n D b B i m a p"
> s2g=t o T a b l e(o r g.H s.e g S Y M B O L)
#求⼩⿏的基因与⼈类的基因的交集，利⽤m a t c h函数，返回位置信息(如果没有对应，就返回N A)。

存在N A的原因就是：⼩⿏有的对应不上⼈类基因名，并且⼈类的基因也有未知的
> g=s2g[m a t c h(s,s2g$s y m b o l),1]
#然后做成⼀个数据框
> d a n n o t=d a t a.f r a m e(p r o b e=s,
"G e n e.S y m b o l" =s,
"E n t r e z G e n e.I D"=g)
#下⾯去掉d d a t a和d a n n o t中N A的⾏
>d d a t a=d d a t a[,!i s.n a(d a n n o t$E n t r e z G e n e.I D)] #I D转换
>d i m(d d a t a)
[1]768 10487 #相⽐之前⼤约去掉2000个基因
> d a n n o t=d a n n o t[!i s.n a(d a n n o t$E n t r e z G e n e.I D),]
#看下去除N A后的基因注释和表达矩阵，必须保证注释的基因I D和表达矩阵的基因I D⼀⼀对应
>h e a d(d a n n o t)
p r o b e G e n e.S y m b o l E n t r e z G e n e.I D
372 A4G A L T A4G A L T53947
393A A A S A A A S8086
394A A C S A A C S65985
396A A G A B A A G A B79719
397A A K1A A K122848
398A A M D C A A M D C28971
>d d a t a[1:4,1:4]
A4G A L T A A A S A A C S A A G A B
S S2_15_0048_A3 8.5163830 0.0000000
S S2_15_0048_A6 7.1119280 0.0000000
S S2_15_0048_A5 3.4154520 0.0000000
S S2_15_0048_A4 6.8487740 7.1681960
可以进⾏g e n e fu分析了，分型就是使⽤m o l e c u l a r.s u b t y p i n g函数
s<-m o l e c u l a r.s u b t y p i n g(s b t.m o d e l = "p a m50",d a t a=d d a t a,
a n n o t=d a n n o t,d o.m a p p i n g=T R U E)
#结果就是将768个细胞
>t a b l e(s$s u b t y p e)
B a s a l H e r2L u m B L u m A N o r m a l
42584654379
#可以利⽤原始的样本信息数据框d f进⾏c l u s t分组与分⼦分型之间关系的探索
> d f$s u b t y p e s=s u b t y p e s
> t a b l e(d f[,c(1,5)])
s u b t y p e s
g B a s a l H e r2 L u m A L u m B N o r m a l
136302051328
23252173124
312102115
42119112
注意：虽然这⾥可以实现分类，但是P A M50是针对乳腺癌患者进⾏分
类的，⽽我们这⾥是针对单细胞；⽽且细胞也不是癌细胞，是
C A F s(c anc er as s oc iat ed f iber blas t）不管是什么细胞，最后都能得
到⼀个表达矩阵，算法是不会考虑矩阵来源的，因此即便是正常细胞的矩阵，也可以分类成5种乳腺癌亚型，所以分类的前提还是⾃⼰熟悉数
据的⽣物学背景
50
看⼀下pam50，它是⼀个列表
> s t r(p a m50)
L i s t o f 7
$ m e t h o d.c o r: c h r "s p e a r m a n"
$ m e t h o d.c e n t r o i d s: c h r "m e a n"
$ s t d: c h r "n o n e"
$ r e s c a l e.q: n u m 0.05
$ m i n s: n u m 5
$ c e n t r o i d s: n u m[1:50, 1:5] 0.718 0.537 -0.575 -0.119 0.3 ...
..- a t t r(*, "d i m n a m e s")=L i s t o f 2
.. ..$ : c h r[1:50] "A C T R3B" "A N L N" "B A G1" "B C L2" ...
.. ..$ : c h r[1:5] "B a s a l" "H e r2" "L u m A" "L u m B" ...
$ c e n t r o i d s.m a p:'d a t a.f r a m e':50 o b s. o f3 v a r i a b l e s:
..$ p r o b e: c h r[1:50] "A C T R3B" "A N L N" "B A G1" "B C L2" ...
..$ p r o b e.c e n t r o i d s: c h r[1:50] "A C T R3B" "A N L N" "B A G1" "B C L2" ... ..$ E n t r e z G e n e.I D: i n t[1:50] 57180 54443 573
596 332 644 891 898 991 990 ...
然后取出基因名，存储在c e n t r o i d s中：
p a m50g e n e s=p a m50$c e n t r o i d s.m a p[c(1,3)]
#发现有的基因已经不是标准的s y m b o l了，P A M50是2009年的基因名，因此需要进⾏修改
p a m50g e n e s[p a m50g e n e s$p r o b e=='C D C A1',1]='N U F2'
p a m50g e n e s[p a m50g e n e s$p r o b e=='K N T C2',1]='N D C80'
p a m50g e n e s[p a m50g e n e s$p r o b e=='O R C6L',1]='O R C6'
以第⼀个基因为例：ht t ps://w w w.genec ar ds.or g/c gi-bin/c ar ddis p.pl?gene=N U F 2&k ey w or ds=N U F2
> x=d a t
#找到p a m50在原始表达矩阵⾏名中的基因，发现⼀共有38个
> p a m50g e n e s$p r o b e[p a m50g e n e s$p r o b e %i n% r o w n a m e s(x)]
[1] "A N L N""B A G1""B C L2""B I R C5""B L V R A""C C N B1"
[7] "C C N E1""C D C20""C D C6""N U F2""C D H3""C E N P F"
[13] "C E P55""C X X C5""E G F R""E R B B2""E S R1""F O X C1"
[19] "K I F2C""N D C80""M A P T""M D M2""M E L K""M I A"
[25] "M K I67""M L P H""M M P11""M Y B L2""M Y C""O R C6"
[31] "P H G D H""P T T G1""R R M2""S F R P1""S L C39A6" "T Y M S"
[37] "U B E2C""U B E2T"
> x=x[p a m50g e n e s$p r o b e[p a m50g e n e s$p r o b e %i n%
r o w n a m e s(x)],]
下⾯进⾏热图可视化
#在原来g r o u p_l i s t基础上，添加亚型信息，为了下⾯p h e a t m a p中的a n n o_c o l设置
t m p=d a t a.f r a m e(g r o u p=g r o u p_l i s t,
s u b t y p e s=s u b t y p e s)
r o w n a m e s(t m p)=c o l n a m e s(x)
#画热图
l i b r a r y(p h e a t m a p)
p h e a t m a p(x,s h o w_r o w n a m e s = T,s h o w_c o l n a m e s = F,
a n n o t a t i o n_c o l = t m p)
图⽚本⾝不重要，因为这⾥数据的使⽤是不合适的。

可以看到，⼤部分
基因都是lum inal A
如果要继续归⼀化就是：
x=t(s c a l e(t(x)))
x[x>1.6]=1.6
x[x< -1.6]= -1.6
p h e a t m a p(x,s h o w_r o w n a m e s = T,s h o w_c o l n a m e s = F,
a n n o t a t i o n_c o l = t m p)
“。