哦别做梦了！

生信技能树

发布于 2020-12-03 14:47:18

4500

发布于 2020-12-03 14:47:18

文章被收录于专栏：生信技能树

前面的学徒作业系列有一个是《数据挖掘》学习班的学员提问：绘图本身很简单但是获取数据很难。本来呢，我是安排给了转录组讲师，希望她可以把这个解决方案制作好PPT给大家做一节公开课的。

但是奈何远程协作沟通是一个问题，所以目前仍然是处于悬而未决的状态。

好在我还有大把的实习生可以安排，群发出去后，七月的优秀实习生回复我了。

七月优秀实习生的答案

写在前面

jimmy老师的作业安排的很简略，

但是我在搜索这个习题的时候，看到了这样一段：生存分析速成指南，赶紧收藏！

点开之后，发现：

真的是气人，啊哈哈哈，果然学习没有捷径啊，我还是得自己独立去解决这个学徒任务。

1.下载数据

首先在 xena 上下载表达矩阵

2.整理数据

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 在正常样本和肿瘤样本的表达量

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()

挑选出匹配的样本

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老师的代码...我自己的太低级，是的我该进步了。

然后我发现接下来的这一步可能有一点点问题：

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 列

所以，我用我拙劣的手法略微的改动了一点点

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),]

这样就得到了两个很匹配的数据...

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:

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

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

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-11-21，如有侵权请联系 cloudcommunity@tencent.com 删除

pid

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

pid

登录后参与评论

0 条评论

热度