前面的学徒作业系列有一个是《数据挖掘》学习班的学员提问:绘图本身很简单但是获取数据很难。本来呢,我是安排给了转录组讲师,希望她可以把这个解决方案制作好PPT给大家做一节公开课的。
但是奈何远程协作沟通是一个问题,所以目前仍然是处于悬而未决的状态。
好在我还有大把的实习生可以安排,群发出去后,七月的优秀实习生回复我了。
七月优秀实习生的答案
写在前面
jimmy老师的作业安排的很简略,
但是我在搜索这个习题的时候,看到了这样一段:生存分析速成指南,赶紧收藏!
点开之后,发现:
真的是气人,啊哈哈哈,果然学习没有捷径啊,我还是得自己独立去解决这个学徒任务。
1.下载数据
首先在 xena 上下载表达矩阵
2.整理数据
代码语言:javascript复制rm(list = ls())
options(stringsAsFactors = F)
if(F){
exp = read.table("data/HiSeqV2.gz",header = T,row.names = 1,check.names = F)
exp = as.matrix(exp)
dim(exp)
exp[1:4,1:4]
index=grep('^TP53$',rownames(exp))
TP53=exp[index,]
save(TP53,file = 'TP53.Rdata')
}
3.观察数据
看看 TP53 在正常样本和肿瘤样本的表达量
代码语言:javascript复制rm(list = ls())
options(stringsAsFactors = F)
library(ggstatsplot)
load('TP53.Rdata')
TP53
sample=names(TP53)
id=substring(sample,14,15)
table(id)
dat=data.frame(sample=sample,TP53=TP53)
dat$group=ifelse(id=='11','Solid tissue normal',
ifelse(id=='01','Primary Tumor','Metastatic'))
table(dat$group)
head(dat)
boxplot(TP53~group,data=dat)
colnames(dat)
###这个给列加因子
dat$group=factor(dat$group,levels = c('Solid tissue normal','Primary Tumor','Metastatic'))
#pdf('TP53-in-BRCA.pdf')
ggbetweenstats(data = dat,
x = group ,
y = TP53 )
dev.off()
挑选出匹配的样本
代码语言:javascript复制dat$pid=substring(dat$sample,1,12)
#首先是table(dat$pid)==2,也就是pid这一列中有多少样本是有2个
#然后挑出来pid这一类TURE
#最后挑出TURE的名字
pairedID=names(table(dat$pid)[table(dat$pid)==2])
##从pid中挑出有两个的配对样本,删掉只有一个的
dat=dat[dat$pid %in% pairedID,]
ggbetweenstats(data = dat ,
x = group ,
y = TP53 )
head(dat)
4.图图图
我抄袭了,哦不我借鉴了Jimmy老师的代码...我自己的太低级,是的我该进步了。
然后我发现接下来的这一步可能有一点点问题:
代码语言:javascript复制dat=dat[order(dat$pid),]
dat=cbind(dat[seq(1,nrow(dat),by=2),],dat[seq(2,nrow(dat),by=2),])
identical(dat[,4],dat[,8])
head(dat)
normal=dat[,2];tumor=dat[,6]
如果是按照pid这一列来排序,取的group这一列却不是匹配的,
pid排序了然而group没有排序,最后就变成了
到最后normal 和 tumor 就不能直接取 第 2 列和第 6 列
所以,我用我拙劣的手法略微的改动了一点点
代码语言:javascript复制index1=grep('normal',dat$group)
normal=dat[index1,]
index2=grep('Tumor',dat$group)
tumor=dat[index2,]
identical(tumor$pid,normal$pid)
tumor$pid = normal$pid[match(tumor$pid,normal$pid)]
## 嗯就记住,谁在括号外面,谁就在后面,x在外面,x就应该在后面...
tumor = na.omit(tumor)
p = identical(tumor$pid,normal$pid);p
if(!p) tumor = tumor[match(normal$pid,tumor$pid),]
这样就得到了两个很匹配的数据...
代码语言:javascript复制d <- data.frame(normal = normal, tumor = tumor)
colnames(d)
library(ggpubr)
ggpaired(d, cond1 = "normal.TP53", cond2 = "tumor.TP53",
fill = "condition", palette = "jco")
所以在肝癌数据中,正常组织和肿瘤组织中 TP53 的表达量并没有多大差异?
是的,基于数据说话。
TP53基因在有的肿瘤患者样本中高表达,在有的肿瘤患者样本中低表达,老师今天特地给我解释了这个问题啊
:happy:
叫辛普森悖论,大概是有高表达,有低表达,合并起来整体来看,却是没有多大差异。之后觉得悖论挺有意思,就去看了一下这本书《悖论简史》~
行吧,就这样了,我去进步了。