要重复的图
image-20190817123355782
老大给的提示：本世纪初Perou和Sorlie 等人作为分子分型(MC)的先驱，通过全基因表达谱确立了浸润性乳腺癌（IBCS）的5个独特亚型： Luminal A、Luminal B、正常乳腺样，HER2基因过表达、基底细胞样，每个亚型各自具有独特的发病率、生存率和疗效。
目前基因表达谱在日常实践中既不经济也不实用，因此许多学者研究采用免疫组织化学发（IHC）作为确定浸润性乳腺癌分子分型的替代方法。
本文讨论已具有临床意义并现实可行的浸润性乳腺癌分子分型 IHC替代法的持续研究成果。
所以需要确认，TCGA的BRCA的 PAM50分型和IHC分型的重合情况。
首先要知道是如何分型的
通过单细胞天地的一篇文章，找到乳腺癌亚型分型的分类
原网址（https://www.breastcancer.org/symptoms/types/molecular-subtypes）显示内容如下
- Luminal A breast cancer is hormone-receptor positive (estrogen-receptor and/or progesterone-receptor positive), HER2 negative, and has low levels of the protein Ki-67, which helps control how fast cancer cells grow. Luminal A cancers are low-grade, tend to grow slowly and have the best prognosis.
- Luminal B breast cancer is hormone-receptor positive (estrogen-receptor and/or progesterone-receptor positive), and either HER2 positive or HER2 negative with high levels of Ki-67. Luminal B cancers generally grow slightly faster than luminal A cancers and their prognosis is slightly worse.
- Triple-negative/basal-like breast cancer is hormone-receptor negative (estrogen-receptor and progesterone-receptor negative) and HER2 negative. This type of cancer is more common in women with BRCA1 gene mutations. Researchers aren’t sure why, but this type of cancer also is more common among younger and African-American women.
- HER2-enriched breast cancer is hormone-receptor negative (estrogen-receptor and progesterone-receptor negative) and HER2 positive. HER2-enriched cancers tend to grow faster than luminal cancers and can have a worse prognosis, but they are often successfully treated with targeted therapies aimed at the HER2 protein, such as Herceptin (chemical name: trastuzumab), Perjeta (chemical name: pertuzumab), Tykerb (chemical name: lapatinib), Nerlynx (chemical name: neratinib), and Kadcyla (chemical name: T-DM1 or ado-trastuzumab emtansine).
- Normal-like breast cancer is similar to luminal A disease: hormone-receptor positive (estrogen-receptor and/or progesterone-receptor positive), HER2 negative, and has low levels of the protein Ki-67, which helps control how fast cancer cells grow. Still, while normal-like breast cancer has a good prognosis, its prognosis is slightly worse than luminal A cancer’s prognosis.
根据图片要求，需要从临床分型中筛选出含有ER/PR(可以把ER和PR合起来理解是HR)，HER2的受体状态的列，及含有PAM50信息的列都找出来。
1.从UCSC xena上下载TCGA BRCA上的数据
image-20190817130141493
2.读取数据
mutype_file <- read.table('BRCA_clinicalMatrix.gz',header = T,sep = '\t',quote = '',stringsAsFactors = F)
colnames(mutype_file)
phetype<-(mutype_file[,c(1,65,70,114,20)])#挑选出来ER\PR、HER2、PAM50分型
phe3<-phetype
phe3
3.处理数据
由于读取进来的数据框中的字符串的长度为0,所以需要去掉，由于不知道如何去掉这些，所以我先把这些但凡字符串长度为0的都变成NA，然后再去掉NA就可以了。
library(stringr)
for (i in 1:nrow(phe3)) {
for (j in 1:ncol(phe3)) {
if(str_count(phe3[i,j])==0){
phe3[i,j]<-NA}
}
}
dim(ph3)
[1] 1247 5
phe4<-na.omit(phe3)
dim(phe4)
[1] 799 5
phe3
4.看一下ER/PR、HER2的阴阳性
#用这个lapply函数，可以把第2到4列中都table下，不用一个一个table。学会这个function(x)代表的意思非常重要，有点不好描述，但是找到感觉，就会用对。lapply是可以用在向量上的，所以我的phe[,]
lapply(2:4, function(x){
table(phe4[,x])
})
[[1]]
Indeterminate Negative Positive
2 176 621
[[2]]
Indeterminate Negative Positive
4 256 539
[[3]]
Equivocal Indeterminate Negative Positive
138 8 499 154
从上面得知，除了Positive、Negative以为，还有Indeterminate还有Equivocal（模凌两可的），需要去掉。
由于Equivocal只存在于第4列，也就是HER2这列，就先去掉第4列中含有Equivocal的所在的样本。
phe5<-phe4[-(grep('Equivocal',phe4[,4])),]
dim(phe5)
[1] 661 5
接下来去掉第2、3、4列中含有Inderterminate所在的样本。
phe5<-phe5[-(grep('Indeterminate',phe5[,2])),]
phe5<-phe5[-(grep('Indeterminate',phe5[,3])),]
phe5<-phe5[-(grep('Indeterminate',phe5[,4])),]
dim(phe5)
[1] 647 5
5.去掉ER和PR不相同的行
phe5
从上一步得到的phe5可以看到，有的ER和PR的阳性和阴性是不同的，而从最上面引用的分类可以知道，HR+是ER和PR共同阳性，HR-是ER和PR共同阴性，是这样的一个分类，至于如果ER+但是PR-的类型，如何判断HR是算阳性还是阴性，我不清楚，也另当别论。所以，现在的目的是去掉那些ER和PR阴性阳性不一致的样本。
phe6<-phe5[,2:3]
d<- do.call(rbind,apply(phe6, 1, function(x){#对phe6的每一行进行操作
a=rep(0,length(nrow(phe6))) #建立一个无字符的值
if(x[1]==x[2]) #进行if判断，如果这一行的第一个字符与第二个字符相等
a=x #就把这一行的内容给a
})) #然乎接着判断第二行。处来的是列表，都do.call,r.bind结合起来
view(d)
d
6.将得到的d数据框痛痛phe5对应起来，因为需要的有HER2列和PAM50列在一起的表达矩阵
其实我是完全可以在phe5中做的，只不过当时怕破坏，才在phe5中提取了ER和PR这两列为phe6
好了，现在可以用之前学过的match函数，从phe5中按照刚刚从d表达矩阵中筛选出的样本的坐标顺序，提取出为一个新的表达矩阵，正好练习下match函数。
e=phe5[match((row.names(d)),row.names(phe5)),]#非常喜欢的match函数，要多练习
e
得到的e和上面的d是一样的，因为就是从phe5里按照表达矩阵d来提取的嘛，当然一样啦！
接下来，新建一列HR,HR也就是ER和PR的阴阳性相同的，所以提取ER这一列，并新建HR一列
e$HR=e[,2]
e
#进行改名，并且将HR、HER2、PAM50这三列提取为一个新列，
colnames(e)[4]='HER2'
colnames(e)[5]='PAM50'
f<-e[,c('HR','HER2','PAM50')]#新表达矩阵f
g<-f #为了不破坏f，从新复制一个f
g
7.将表达矩阵中的向量进行改名
g[,1][as.data.frame(g[,1])=='Positive']<-'HR+'
g[,1][as.data.frame(g[,1])=='Negative']<-'HR-'
g[,2][as.data.frame(g[,2])=='Positive']<-'HER2+'
g[,2][as.data.frame(g[,2])=='Negative']<-'HER2-'
g
插曲：这步尝试了比较多的次数
第一种：用gsub、paste函数
第二种：apply循环
- 第一种
h<-gsub('Positive','HR+',g[c('HR')])
h<-gsub('Negative','HR-',h)
h
i<-gsub('Positive','HER2+',g[c('HER2')])
i<-gsub('Negative','HER2-',i)
i
>h
h
>i
i
h和i都是长度为1的向量。并且还含有\n我用"，"拆分以后，想在用paste再连接起来，一直无法一对一的连接，就是第一个HR连接第一个HER2,第二个HR连接第二个HER2.
paste函数
所以这个我就放弃了
- 第二种
####下面这个是把HR列和HER2都变成了HR-，这个循环不行
apply(g, 2, function(x){gsub('Negative','HR-',x)})
最后，找到了可以再data.frame中修改向量的方法，参考的是网上搜索到的这样一个方法：
拥有数据框架，如何替换所有行和列的所有特定值。说例如我想用NA替换所有空记录(不输入位置)：
df <- data.frame(list(A=c("","xyz","jkl"), B=c(12,"",100)))
A B
1 12
2 xyz
3 jkl 100
##想用NA替换上面的null的地方
> df[df==""]<-NA
> df
A B
1 <NA> 12
2 xyz <NA>
3 jkl 100
用的是df[df==""]<-NA这样一个方法，我的修改后适应自己的要求,改了如下。
g[,1][as.data.frame(g[,1])=='Positive']<-'HR+'
有一点，就是要改四次。
8.把HR和HER2两列合并为一列
在网上搜索到的都是如何把两列数据框合在一起，比如cbind，merge函数。但是可以看最前面要画的图图是HR+/HER+、HR+/HER-、HR-/HER+、HR-/HER-（TNBC）这么组合的，所以现在的目的是在R中合并某一数据框的两列数据，因此cbind、merge函数统统排不上用场了，用tidyr包中的unite函数，解决问题。
library(tidyr)
k<-unite(g, "HR/HER2", HR, HER2, sep = "/", remove = FALSE)
#函数本身一目了然，是tidyr包里面的，其中remove=FALSE是说新建一列，而不是在原有的基础上改动
k-最后处理完的表达矩阵
9.画图
插曲：为什么要把数据处理成上面k那样子？
是因为我看来ggstatsplot包
里面的一个example是这样的
ggstatsplot::ggpiestats(
data = datasets::mtcars,
main = vs,
condition = cyl,
bf.message = TRUE,
nboot = 10,
factor.levels = c("0 = V-shaped", "1 = straight"),
legend.title = "Engine"
)
datasets::mtcars
example出图
library(ggstatsplot)
k$subtype <- factor(k$PAM50,levels = c("Normal","Her2","LumB","LumA","Basal"))
ggpiestats(data = k,
main = subtype,
condition =HR/HER2,
palette = "Set1" )
但是，我现在的是报错。
报错
还没有发现问题所在。
网友评论