Posted on

GTEx(Genotype-Tissue Expression,基因型-组织表达)数据库,研究从来自449名生前健康的人类捐献者的7000多份尸检样本,涵盖44个组织(42种不同的组织类型),包括31个实体器官组织、10个闹分区、全血、2个来自捐献者血液和皮肤的细胞系,作者利用这些样本研究基因表达在不同组织和个体中有何差异。

GTEx对几乎所有转录基因的基因表达模式进行了观察,从而能够确定基因组中影响基因表达的特定区域。

此外,合并GTEx与TCGA数据库数据能够有效解决TCGA数据库中正常组织样本量不足的缺陷,从而提高比较的准确性。

1. 数据来源

(数据比较大,如果下载困难可以留言)

2. 注释来自 TCGA GTEx 的样本

library(stringr)
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(data.table)


#################======= step1: clean GTEx pheno data   =======#################
gtex <- read.table("samplepair.txt",header=T,sep='\t')

tcga_ref <- gtex[,1:2]

gtex$type <- paste0(gtex$TCGA,"_normal_GTEx")
gtex$sample_type <-"normal"
gtex <- gtex[,c("TCGA","GTEx","type","sample_type")]
names(gtex)[1:2] <- c("tissue","X_primary_site")


gp <- read.delim(file="GTEX_phenotype.gz",header=T,as.is = T)
gtex2tcga <- merge(gtex,gp,by="X_primary_site")
gtex_data <- gtex2tcga[,c(5,2:4)]
names(gtex_data)[1] <- "sample"
#write.table(gtex_data,"GTEx_pheno.txt",row.names=F,quote=F,sep='\t')

#################======= step2: clean a TCGA pheno data   =======#################
tcga <- read.delim(file="TCGA_phenotype_denseDataOnlyDownload.tsv.gz",header=T,as.is = T)
tcga <- merge(tcga_ref,tcga,by.y="X_primary_disease",by.x="Detail",all.y = T)
tcga <- tcga[tcga$sample_type %in% c("Primary Tumor","Solid Tissue Normal"),]

tcga$type <- ifelse(tcga$sample_type=='Solid Tissue Normal',
                    paste(tcga$TCGA,"normal_TCGA",sep="_"),paste(tcga$TCGA,"tumor_TCGA",sep="_"))
tcga$sample_type <- ifelse(tcga$sample_type=='Solid Tissue Normal',"normal","tumor")
tcga<-tcga[,c(3,2,6,5)]
names(tcga)[2] <- "tissue"
#write.table(tcga,"tcga_pheno.txt",row.names = F,quote=F,sep='\t')


#################======= step3: remove samples without tpm data =======############
gtex_exp <-  fread("gtex_RSEM_gene_tpm.gz",data.table = F)
gtexS <- gtex_data[ gtex_data$sample%in%colnames(gtex_exp)[-1],]

tcga_exp <- fread("tcga_RSEM_gene_tpm.gz",data.table = F)
tcgaS <- tcga[tcga$sample %in%colnames(tcga_exp)[-1],]
tcga_gtex <- rbind(tcgaS,gtexS)
write.table(tcga_gtex,"tcga_gtex_sample.txt",row.names = F,quote=F,sep='\t')

3. 提取感兴趣的基因

library(stringr)
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(data.table)
library(tibble)

rm(list = ls())
options(stringsAsFactors = FALSE) 

target <- "YTHDC2"

idmap <- read.delim("gencode.v23.annotation.gene.probemap",as.is=T)
tcga_exp <- fread("tcga_RSEM_gene_tpm.gz",data.table = F)
gtex_exp <- fread("gtex_RSEM_gene_tpm.gz",data.table=F)
tcga_gtex <- read.table("tcga_gtex_sample.txt",sep='\t',header = T)

id <- idmap$id[which(idmap$gene==target)]
tcga_data <- t(tcga_exp[tcga_exp$sample==id,colnames(tcga_exp)%in%c("sample",tcga_gtex$sample)])
tcga_data <- data.frame(tcga_data[-1,])
tcga_data <- rownames_to_column(tcga_data,"sample")
names(tcga_data)[2] <- "tpm"

gtex_data <- t(gtex_exp[gtex_exp$sample==id,colnames(gtex_exp)%in%c("sample",tcga_gtex$sample)])
gtex_data <- data.frame(gtex_data[-1,])
gtex_data <- rownames_to_column(gtex_data,"sample")
names(gtex_data)[2] <- "tpm"

tmp  <- rbind(tcga_data,gtex_data)
exp <- merge(tmp,tcga_gtex,by="sample",all.x=T)
exp <- exp[,c("tissue","sample_type","tpm")]
exp <- arrange(exp,tissue)
write.table(exp,"Merge gene expression/YTHDC2 expression.txt",row.names = F,quote=F,sep='\t')

4. 可视化基因表达

library(ggplot2)
library(ggpubr)
library(RColorBrewer)

rm(list = ls())
options(stringsAsFactors = FALSE) 

exp <- read.table("Merge gene expression/YTHDC2 expression.txt",header=T,sep='\t')

ylabname <- paste("YTHDC2", "expression")
colnames(exp) <- c("Tissue", "Group", "Gene")

p1 <- ggboxplot(exp, x = "Tissue", y = "Gene", fill = 'Group',
                ylab = ylabname,
                color = "Group", 
                palette = c("#00AFBB",  "#FC4E07"),
                ggtheme = theme_minimal())

##计算每种肿瘤正常和肿瘤组织的样本量
count_N<-exp %>% group_by(Tissue, Group) %>% tally
count_N$n <- paste("n =",count_N$n)
##添加N = 到图中
p1 <-p1+geom_text(data=count_N, aes(label=n, y=-9,color=Group), position=position_dodge2(0.9),size = 3,angle=90, hjust = 0)+
  theme(axis.text.x = element_text(angle = 45,hjust = 1.2))

#计算t检验显著性
comp<- compare_means(Gene ~ Group, group.by = "Tissue", data = exp,
                         method = "t.test", symnum.args = list(cutpoints = c(0,0.001, 0.01, 0.05, 1), symbols = c( "***", "**", "*", "ns")),
                         p.adjust.method = "holm")
#添加显著性标记
p2 <- p1 + stat_pvalue_manual(comp, x = "Tissue", y.position = 7.5,
                     label = "p.signif", position = position_dodge(0.8))
p2
#dev.off()

##保存图片
### pdf version
ggsave("figure/pancancer_Plot.pdf", width = 14, height = 5)

### png version
#png("figure/pancancer_Plot.png", width = 465, height = 225, units='mm', res = 300)

代码参考GitHub:https://github.com/cmutd/TCGA_GTEx

71 Replies to “GTEx联合TCGA数据库差异分析(更新)”

  1. 进哥你好,我想咨询一下做弥漫大b淋巴瘤,GTEx里应该提取那个组织的表达组数据啊,应该是淋巴结吧,我没有找到。。

    1. 额,这种血液肿瘤就直接以blood做对照吧,最好是以同种组织的正常组织作为对照,没有的情况下找相近的

  2. 进哥哥,您最后用了t.test,可以认为tcga和gtex两组间是满足正态方差齐吗?

    1. 您好,终于有人发现了这个统计上的漏洞。事实上,大部分癌症类型是满足的,对于部分样本少,或者有缺失数据的可能就不满足正态方差齐,这时候用Wilcoxon秩和检验更加合适,选择t.test只是图个方便。t检验和Wilcoxon秩和检验

  3. 博主您好,我在把YTHDC2换了一个目的基因,出的图会变形,没办法调整Y轴

    1. 您好,你指的变形是样本量的位置是吗?这一句调整位置,改一下y=:geom_text(data=count_N, aes(label=n, y=-9,color=Group)
      如果其它问题可以详细描述一下 或者加我微信

    1. 你想要一次分析几个基因在33种癌症中的表达吗?要做当然可以,关键看你最终将以什么形式呈现,可以加微信讨论

      1. gtex_exp <- fread("gtex_RSEM_gene_tpm.gz",data.table = F)文件过大,无法打开,更改data.table=T勉强可以,想问进哥哥有啥办法只提取这个大数据里我需要的数据即可,不必完全打开。

        1. 你好 之前评论回复里面有分割的GTEx和Pancancer数据,你可以只下载你需要的组织类型数据,然后进行合并分析,如果还是需要多种癌症,写个循环即可

          1. 是这样的,我用的gz文件和进哥哥展示的还不太一样,我使用的是transcrip expression的,所以我是不是得自己分割了,,,,

  4. 进哥,我还有一个问题,THYM正常对照可以用blood吗?虽然胸腺里面是免疫细胞,但是和全血还是差别很大的

    1. 您好 这个我不清楚 讲道理应该不可以,你看文献怎么处理的,如果实在没有正常组织,可以考虑分析不同分期的差异

      1. tcga<-tcga[,c(3,2,6,5)]
        Error in `[.data.frame`(tcga, , c(3, 2, 6, 5)) : 选择了未定义的列

        进哥,我按照你的代码整理tcga表型文件的时候提取3,2,6,5例时显示有一例不存在,我看了tcga table只有5列,detail ,TCGA, sample, sample_type_id, sample_type.这该怎么解决?

        1. tcga$type <- ifelse(tcga$sample_type=='Solid Tissue Normal', paste(tcga$TCGA,"normal_TCGA",sep="_"),paste(tcga$TCGA,"tumor_TCGA",sep="_")) 这一行运行了吗 还是运行报错了?因为你没有type这一列

  5. 进哥,像UVM这种联合GTEX都拿不到正常组织数据的癌种,怎么分析基因差异呢?

  6. 请问您的TCGAphenotype文件是什么样的来源?是进行过处理的吗?还是从xena直接下载的呢?如果进行过处理请问是怎么样处理的呢?

    1. 您好 这个文件是xena直接下载的,保存的是癌组织来源,没有进行处理,直接与表达数据合并

  7. 张博士,你好,请问下DLBCL中GTEx的数据用的是全血细胞吗?可以把DLBCL的TCGA数据和GTEx数据发我吗?

  8. 您好,请问可以把GTEx的数据和GEO的数据合并分析吗?需要下载什么样的数据?研究的病理亚型在TCGA 样本量太少了

    1. 合并的话就是保证使用的TCGA和GTEx数据集是同一种标准化方式,比如这篇文章里用的是log2(TPM+0.001)。
      但是您说所研究的病理亚型在TCGA样本量太少,那你合并GTEx也没有帮助哇,GTEx都是正常组织。如果TCGA无法满足,可以尝试搜索相关GEO数据库或一些其他肿瘤数据库,比如ICGC等

  9. 您好,进哥。我想咨询一下,您的这些数据有正常脑组织样本数据吗?与 MiRNA 有关的。期待并感谢您的解答

    1. 正常脑组织数据有的,但是GTEx数据库只有genotype和gene expression数据,应该是没有mirna expression数据

  10. 您好!想要您 “tcga_RSEM_gene_tpm.gz ”这个数据,因为文件很大,下载不下来,非常感谢!

    1. 你好,我这边有下载好的,可是一样很大,网盘发你下载也会很慢很慢。另外我有根据组织类型分割好的,前面评论有链接

      1. 进哥,如果方便,恳请发一份卵巢癌及子宫内膜癌的TCGA+GTEx整合的表达矩阵文件给我吧,谢谢!

      2. 谢谢您!想问一下分割好的数据,即来自TCGA和GTEx两个都是log2(TPM)后的吗?

  11. 博主你好,求一份卵巢癌TCGA+GTEx整合的表达矩阵文件,谢谢!

      1. 我跑上面的代码出现这一行错误:
        > write.table(exp,”Merge gene expression/YTHDC2 expression.txt”,row.names = F,quote=F,sep=’\t’)
        Error in file(file, ifelse(append, “a”, “w”)) :
        cannot open the connection
        In addition: Warning message:
        In file(file, ifelse(append, “a”, “w”)) :
        cannot open file ‘Merge gene expression/YTHDC2 expression.txt’: No such file or directory

        1. 因为你路径下不存在”Merge gene expression/YTHDC2 expression.txt”,你需要改成你自己电脑上保存的路径及文件名

          1. 把YTHDC2换了一个目的基因,出的图会变形,需要调整哪个参数?

  12. 博主,我可以求一个卵巢相关的GTEx的数据吗,GTEx的数据我试了很多种方法都下不了。非常感谢!!!!!!

  13. 博主您好,请问我想要看看Xp上的某些基因在不同组织中XCI ratio,要怎么做啊,谢谢

      1. 之前读取数据的时候内存不够,所以我用了memory.limit(newLimit)修改了内存,但是现想把这个再降回来,现在R语言 用了太多内存。怎么给清空呢。
        used (Mb) gc trigger (Mb) max used (Mb)
        Ncells 6880145 367.5 11096502 592.7 11096502 592.7
        Vcells 14379954 109.8 409974109 3127.9 903763272 6895.2

  14. 请问一下,单个基因进行差异分析,tpm数据需要经过处理吗,有些是负值。

    1. 你好 这里面用的数据是经过log2(tpm+0.001)处理的,对于有些表达量低的,转化之后是负数。所以 这些数据不需要进行其他处理,直接可以使用

  15. 进哥GTEx与TCGA连用的时候是不是去除批效应,GEPIA那个网站去除了吗?

  16. 博主您好,我下不了GTEx的数据,请问您可以提供下吗,谢谢

    1. Xena数据库国内访问没应该有问题,是下载太慢吗?我有下载好的,但是一样很大,放网盘再给你下载一样也很慢,我这边也有按照肿瘤类型分割好的数据,如果不做泛癌我可以发你需要的癌症数据

      1. 我需要BRCA(乳腺癌)跟THCA(甲状腺癌)两个类型的数据,可以发我邮箱吗?谢谢

          1. 我需要TCGA+GTex联合的数据,根据肿瘤类型分割好的数据,就只要BRCA跟THCA 这两个瘤种数据

          2. 链接:https://pan.baidu.com/s/17blyTZb-Kni8u9yqIsweyA?pwd=eqsp
            提取码:eqsp
            在家里电脑又分割了一下,我把GTEX和pancancer的所有分割好的数据发你,你自己下载下来合并分析,在R里面用merge根据Ensemble ID合并即可

  17. 王博士您好!我想请问一下计算t检验显著性时会报错Error in `mutate()`:
    ! Problem while computing `p = purrr::map(…)`.
    Caused by error in `if (stderr < 10 * .Machine$double.eps * max(abs(mx), abs(my))) …`:
    ! missing value where TRUE/FALSE needed
    Run `rlang::last_error()` to see where the error occurred.
    是什么原因?谢谢!

    1. 您好,这个问题应该是你的表达量那一列非数值型,你查看一下基因表达量那一列是否有NA,也看一下这一列的数据类型是否为numberic,如果有,可以Excel删除这些NA值再导入R,或者在R中用as.numberic()把表达量那一列转换一下,若有问题,可以加微信讨论

  18. 请问一下,我看到在XENA网站上下载的GTEx和TCGA数据的TPM计算方式不一样(GTEX是log2(TPM+0.001),TCGA是log2(TPM+1)),可以直接合并进行分析吗?不用转化成统一的格式吗?还有我看gepia网站的差异分析中是没有负值的,但是我做的TPM中是有负值的,log2转化后就是nan,该怎么办?

  19. 博主您好,我想问一下在GTEX中下载的这个gtex_RSEM_gene_tpm.gz 正常样本中,来自血液的样本编号是什么呢? 或者是骨髓来源的 样本编号是什么? 网上找了半天也没找到,恳请您的答复。

    1. 您好,这个样本注释信息在GTEX_phenotpye.gz 文件中有,有关于样本来源组织的信息(包括whole blood)
      这个文件文中有链接下载

  20. 写的很棒,学到了很多,但是下面这两个一起跑会显示:Error: cannot allocate vector of size 523 Kb,始终没有找到解决的办法?求大神给解答一下。谢谢
    gtex_exp <- fread("gtex_RSEM_gene_tpm.gz",data.table = F)
    gtexS <- gtex_data[ gtex_data$sample%in%colnames(gtex_exp)[-1],]

    tcga_exp <- fread("tcga_RSEM_gene_tpm.gz",data.table = F)
    tcgaS <- tcga[tcga$sample %in%colnames(tcga_exp)[-1],]

    1. 网上找的教程,可以试试:
      cannot allocate vector就是典型的数据太大读不了
      一、升级硬件
      二、改进算法
      三、修改操作系统分配给R的内存上限, memory.size(T)查看已分配内存

      memory.size(F)查看已使用内存

      memory.limit()查看内存上限

      object.size()看每个变量占多大内存。
      memory.size()查看现在的work space的内存使用
      memory.limit()查看系统规定的内存使用上限。

      如果现在的内存上限不够用,可以通过memory.limit(newLimit)更改到一个新的上限。注意,在32位的R中,封顶上限为4G,无法在一个程序上使用超过4G (数位上限)。这种时候,可以考虑使用64位的版本。

      1. 谢谢师兄,我自己试了一下,跑不出结果,哭了
        1.注释第三步的时候出现报错
        > tcga_exp source(“GeomSplitViolin.R”)
        Error in file(filename, “r”, encoding = encoding) : 无法打开链结
        此外: Warning message:
        In file(filename, “r”, encoding = encoding) :
        无法打开文件’GeomSplitViolin.R’: No such file or directory
        这个我完全看不懂,后面几乎都报错,也有“无法配置矢量”的报错
        谢谢师兄指点!

        1. source(“GeomSplitViolin.R”)这一句多余的,原本这个文件里面是一个分裂小提琴图的函数,后面并没有用到。后面不应该会报错,你可以加我微信讨论,在我的简历里面有手机号和二维码

发表评论

邮箱地址不会被公开。 必填项已用*标注