生不了啦GEO SRA ArrayExpress ENA 数据挖掘

关于实验中要选择的细胞系,你是课题组有啥用啥么?

2019-06-29  本文已影响55人  9d760c7ce737

做实验的时候,如果想说清楚靶分子的重要性,常规操作就是在表达量低的细胞系中过表达基因,在表达量高的细胞系中敲减靶分子,如果观察到相反的现象,那么事情就讲清楚了。

接下里就是确定要研究的细胞系,一般情况下就是收集齐细胞系,然后做一个qRT-PCR就行了,如果是编码基因,最好用western blot。

那现在问题有两个:

  1. 第一,极大的可能,收集不全要研究的细胞系
  2. 第二,如果大家做过qRT-PCR,就知道他的条件极其苛刻,波动性很大。

有没有数据,可以让我们方便地查看基因在多个细胞系中的表达量呢?
有的,有人把所有肿瘤细胞系都测了RNA-seq,而且把数据好上传了。数据存放在CCLE


点进去界面很简洁,我们可以输入一个基因XBP1, 测试一下。

这时候,他给我们返回了这个基因在所有组织分类中的表达,比如,我们发现他在乳腺癌的60个细胞系中,表达算高的。如果我想查看这个基因在所有乳腺癌细胞系中的表达,怎么办呢?
好像不可以。

1.那我们就下载数据自己搞

点击data


里面有很多文件,我们需要的是两个,一个是转录本水平的TPM文件。

一个是样本的注释文件,往下拉就有了


2.读入表达量数据

CCLE_data <- data.table::fread("CCLE_RNAseq_rsem_transcripts_tpm_20180929.txt",data.table =F )

这个数据很大,包含了1000多个细胞系的结果。


选择breast乳腺癌细胞系所在的位置

index <- grep("BREAST",names(CCLE_data))

接下来要做的,就是把这个数据变成能作图的清洁数据(这是R语言中最重要的能力)

library(tidyr)
library(tibble)
library(dplyr)
data <- CCLE_data %>% 
  select(c(2,index)) %>% 
  separate(transcript_id,into = "transcript_id",sep = "\\.") %>% 
  column_to_rownames("transcript_id") %>% 
  t() %>% 
  data.frame() %>% 
  rownames_to_column("celline") %>% 
  separate(celline,into = "celline",sep = "\\_")

这串代码很流畅,也很抽象,必须自己运行才能理解,我们直接看现在的结果吧。
第一列是细胞系的名称,其余列都是一个个转录本,只是用ensemble ID来表示了。


现在还有个问题,就是细胞系的名称,不规范,比如BT474,应该为BT-474, 所以我们需要注释文件来处理

3.修改细胞系的名称

anno_CCLE <- data.table::fread("anno_CCLE.txt")
names(anno_CCLE)[5] <- "Site_Primary"
names(anno_CCLE)[1] <- "CCLE_name"
names(anno_CCLE)[2] <- "Cell_line_primary_name"
library(dplyr)
breast_ccle_index <- anno_CCLE %>% 
  filter(Site_Primary=="breast") %>% 
  select(CCLE_name,Cell_line_primary_name) %>%
  separate(CCLE_name,into = "celline",sep = "\\_") %>% 
  column_to_rownames("celline")
data <- cbind(cell_name=breast_ccle_index[data$celline,1],data)
save(data,file = "CCLE_BREAST_PLOT.Rdata")

现在数据就改过来了


4.选择转录本

首先要确定XBP1的转录本,大部分人在这里就崩溃了。
一个基因会被转录成多个转录本,这些转录本的CDS可能不一样,所以,过表达一个基因,实际上过表达选定转录本的CDS区。这个可以从ensemble数据库得到。


这时候我们会看到,XBP1有很多转录本,而且不同转录本对应的编码蛋白大小还不一样,所以这时候是要读文献确定的。我们这里随便选一个ENST00000216037

5.细胞系作图

library(ggplot2)
gene_id <- "ENST00000216037"
test <- data[,c("cell_name",gene_id)]
names(test)[2] <- "gene_id"
ggplot(data = test,aes(x=forcats::fct_reorder(cell_name,desc(gene_id)),y=gene_id))+
  geom_bar(stat='identity', aes(fill=cell_name))+
  theme_bw()+
  ylab(paste0(gene_id," (TPM)"))+
  xlab(NULL)+
  theme(legend.position = "none")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

50几个细胞系的结果一目了然。
当然,如果先细胞系太多了,自己组上没有那么多,可以选择部分来做图。

library(ggplot2)
gene_id <- "ENST00000216037"
test <- data[,c("cell_name",gene_id)]
rownames(test) <- test$cell_name
test <- test[c("T-47D","SK-BR-3","BT-549","ZR-75-1","MDA-MB-231","MDA-MB-468","MCF7","BT-474","MDA-MB-436"),]
names(test)[2] <- "gene_id"
ggplot(data = test,aes(x=forcats::fct_reorder(cell_name,desc(gene_id)),y=gene_id))+
  geom_bar(stat='identity', aes(fill=cell_name))+
  theme_bw()+
  ylab(paste0(gene_id," expression (TPM)"))+
  xlab(NULL)+
  theme(legend.position = "none")+ 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

此时我们选择"ZR-75-1"来敲减,选择"MDA-MB-231"来过表达就没有人来说什么了吧,如果你只想在一个细胞系中来做实验,选择中等表达的,可以敲减可以过表达,那么MCF-7就是个好选择。

总结

今天的技能,属于那种没有人教,需要自己学会的技能。在科研工作中有很多高光时刻都属于这一类,在R语言处理数据的时候,也是这样。
学一个知识,处理上百种情况, 这在成为了医学生后几乎不可能,所有的知识都需要背,自己推理的结果往往跟实际情况不一样。但是在学习了R语言之后,情况好了很多,在科研中有时候用一下R语言会很有成就感。如果你也想拥有这种能力,可以看看下面这个实诚的帖子:
久等了,果子学生信最新课程报名通知

上一篇下一篇

猜你喜欢

热点阅读