哦别做梦了!

2020-12-03 14:47:18 浏览数 (1)

前面的学徒作业系列有一个是《数据挖掘》学习班的学员提问:绘图本身很简单但是获取数据很难。本来呢,我是安排给了转录组讲师,希望她可以把这个解决方案制作好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:

叫辛普森悖论,大概是有高表达,有低表达,合并起来整体来看,却是没有多大差异。之后觉得悖论挺有意思,就去看了一下这本书《悖论简史》~

行吧,就这样了,我去进步了。

pid

0 人点赞