专栏首页生信技能树TCGA数据库LUSC亚型批量差异分析

TCGA数据库LUSC亚型批量差异分析

前些天我布置了一个学徒作业:这一个图背后是12个差异分析的综合

作业参考的文献:Integrated analysis reveals five potential ceRNA biomarkers in human lung adenocarcinoma

所以我设置的学徒作业是:下载TCGA数据库中LUSC的转录组信号值矩阵,LUSC病人分成了4类T1-4亚型分别与Normal组做差异分析,就是3*4=12个表达矩阵,12次差异分析,画PCA图,热图,火山图,以及用于差异分析结果比较的Venn图。

下面让我们一起看看一个优秀学徒的表演,该学徒很久以前在我们这里分享过他跨专业进入生信学习圈子的感悟:在华大工作五年还不如生信技能树3天?

下载数据

紧跟群主的TCGA视频课程,从UCSC的XENA下载LUSC表达矩阵,临床信息,探针注释GMT文件!

视频课程见:4个小时TCGA肿瘤数据库知识图谱视频教程又有学习笔记啦

rm(list = ls())
pd=read.table("TCGA-LUSC.GDC_phenotype.tsv.gz",header=T,sep = "\t", quote = "\"",dec = ".", fill = TRUE, comment.char = "",row.names = 1)
gset_mRNA=read.table("TCGA-LUSC.htseq_counts.tsv.gz",header=T,sep = "\t", quote = "\"",dec = ".", fill = TRUE, comment.char = "",row.names = 1)
gset_miRNA=read.table("TCGA-LUSC.mirna.tsv.gz",header=T,sep = "\t",row.names = 1)

#由于导入数据列名自动把“-”变为了“.”,可用gsub替换回来。

colnames(gset_miRNA)=gsub("\\.","-", colnames(gset_miRNA))
colnames(gset_mRNA)=gsub("\\.","-", colnames(gset_mRNA))

# 去掉mRNA和miRNA表达矩阵不一致的样本

table(colnames(gset_mRNA) %in% colnames(gset_miRNA))
gset_mRNA=gset_mRNA[,colnames(gset_mRNA) %in% colnames(gset_miRNA)]
gset_miRNA=gset_miRNA[,colnames(gset_miRNA) %in% colnames(gset_mRNA)]

table(substring(colnames(gset_mRNA),14,16))
table(pd$pathologic_T)  

# TCGA数据库的样本分组规律(感谢Jimmy大神提供):Tumor types range from 01 - 09, normal types from 10 - 19 and control samples from 20 - 29.

# 分期里T1分为了T1a,T1b,T1;T2分为了T2a,T2b,T2;这里就不细分了,分别合在一起组成T1和T2期亚群

group_list_mRNA=ifelse(substring(colnames(gset_mRNA),14,15) == "11","Normal",ifelse(colnames(gset_mRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T1",]),"T1",ifelse(colnames(gset_mRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T1",]),"T1",ifelse(colnames(gset_mRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T2",]),"T2",ifelse(colnames(gset_mRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T3",]),"T3",ifelse(colnames(gset_mRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T4",]),"T4","q"))))))
group_list_miRNA=ifelse(substring(colnames(gset_miRNA),14,15) == "11","Normal",ifelse(colnames(gset_miRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T1",]),"T1",ifelse(colnames(gset_miRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T1",]),"T1",ifelse(colnames(gset_miRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T2",]),"T2",ifelse(colnames(gset_miRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T3",]),"T3",ifelse(colnames(gset_miRNA) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == "T4",]),"T4","q"))))))


table(group_list_mRNA)                                 
table(group_list_miRNA) 

# 查下TCGA官网数据的说明,mRNA表达矩阵取了log2(fpkm+1),miRNA表达矩阵取了log2(RPM+1),所以这里要还原回去

head(gset_miRNA)
head(gset_mRNA)
gset_miRNA=(2^gset_miRNA -1)*1000  #差异分析用的RPKM,所以要乘以1000
gset_mRNA=2^gset_mRNA -1

# mRNA表达矩阵包括codingRNA和lncRNA,需要拆分下

library(rtracklayer)
library(SummarizedExperiment)
gtf1 <- rtracklayer::import('gencode.v22.annotation.gtf/gencode.v22.annotation.gtf')
gtf_df <- as.data.frame(gtf1)
head(gtf_df)
ensem2symbol <- gtf_df[gtf_df$type == 'gene',c('gene_id', 'gene_type', 'gene_name', 'source')]
rownames(ensem2symbol) <- ensem2symbol$gene_id
head(sort(table(ensem2symbol$gene_type),decreasing = T))

gset_cdRNA=gset_mRNA[rownames(gset_mRNA) %in% rownames(ensem2symbol[ensem2symbol$gene_type == "protein_coding",]),]
gset_lncRNA=gset_mRNA[rownames(gset_mRNA) %in% rownames(ensem2symbol[ensem2symbol$gene_type == "lincRNA",]),]


#保存表达矩阵和分组信息
save(pd,ensem2symbol,gset_miRNA,gset_cdRNA,gset_lncRNA,group_list_mRNA,group_list_miRNA,file = "step1_output.Rdata")
  • gsub批量整体替换字符很方便
  • ifelse函数筛选T1-T4的样本ID,得到表达矩阵及分组信息
  • 用基因探针GMT文件注释拆分mRNA表达矩阵成cdRNA(编码蛋白的基因)和lncRNA表达矩阵
  • 注意TCGA上对表达矩阵的格式说明,DESeq2差异分析是对count值表达矩阵,必要需要还原!

DESeq2差异分析

1.比较LUSC患者T1-4分型与正常样本的差异基因或miRNA RNA表达矩阵

1.1 检查数据
## 全部肿瘤样本及正常样本表达矩阵的PCA图,热图

rm(list = ls())
load(file = "step1_output.Rdata")
dat=gset_cdRNA
group_list=ifelse(substring(colnames(gset_cdRNA),14,15) == "11","Normal","Tumor")

### 去掉低质量的探针行
dat=dat[apply(dat,1, function(x) sum(x>1) > 50),]

### 检查数据
table(group_list)
source("run_check_h_pca.R")
pro = "cdRNA_Tumor_vs_Normal"
run_check_h_pca(pro)
  • 样本分组 GroupNormalT1T2T3T4样本个数381062796921
  • 全部Tumor样本和Normal组的热图和PCA图可以看出,Tumor组样本大都与Normal组有显著差异,从而可进行下一步差异分析。
1.2 差异分析
## T1-T4亚型组与正常组表达矩阵分别差异分析

### 去掉低质量的探针行
gset=gset_cdRNA
gset=gset[apply(gset,1, function(x) sum(x>1) > 50),]

### for循环批量差异分析
source("function.R")
for (Typel in c("T1","T2","T3","T4")) {
  pro=paste0("cdRNA_",Typel,"_vs_Normal_2")
  run_filter_check_DESeq2(Typel,pro)
}
1.3 比较差异分析结果
## 比较T1-4分别与正常样本差异基因情况
rm(list=ls())
load("cdRNA_T1_vs_Normal_deseq2_DEG.Rdata")
p0=p1 # 一个赋值小错误,由于包装函数里是p1,所以导致赋值时p1和p4一样,所以暂且改成p0赋值了
DEG_T1=DEG
load("cdRNA_T2_vs_Normal_deseq2_DEG.Rdata")
p2=p1
DEG_T2=DEG
load("cdRNA_T3_vs_Normal_deseq2_DEG.Rdata")
p3=p1
DEG_T3=DEG
load("cdRNA_T4_vs_Normal_deseq2_DEG.Rdata")
p4=p1
DEG_T4=DEG

library(ggplot2)
library(gridExtra)
plots <- list(p0,p2,p3,p4)
p= grid.arrange(grobs = plots, ncol = 2,
                top = "compare")
ggsave(p,filename = "cdRNA_compare_volcano.png",width = 20,height = 15)


venn <- function(x,y,z,a,name){
  if(!require(VennDiagram))install.packages('VennDiagram')
  library (VennDiagram)
  venn.diagram(x= list(T1_vs_Normal = x,T2_vs_Normal = y,T3_vs_Normal = z,T4_vs_Normal = a),filename = "cdRNA_DEGgene.png", 
               height = 1000, width = 1000,resolution =300, imagetype="tiff", col="transparent",fill=c("cornflowerblue","green","yellow","darkorchid1"),alpha = 0.5, cex=0.45, cat.cex=0.45)
}

DEG_T1$change!="NOT"

DEGs=function(df){
  rownames(df)[df$change!="NOT"]
}

library(VennDiagram)
venn(DEGs(DEG_T1),DEGs(DEG_T2),DEGs(DEG_T3),DEGs(DEG_T4),
     "DEGgene"
)

DEGsl = intersect(intersect(intersect(DEGs(DEG_T1),DEGs(DEG_T2)),DEGs(DEG_T3)),DEGs(DEG_T4))

length(DEGsl)

2. lncRNA和miRNA表达矩阵一样批量分析

  • 这里就直接上文献中类似的venn图结果
  • T1-4期患者样本分别与正常样本差异分析的阈值:log2FC=1,FDR=0.01
  • T1-4期患者样本分别与正常样本差异分析结果
  • cdRNA:19814个基因里有5573个共同差异基因
  • lncRNA:7656个基因里有1571个共同的差异lncRNA基因
  • miRNA:1881个miRNA里有164个共同的差异miRNA

总结

  • 首先特别感谢Jimmy 大神,以上函数代码均引自Jimmy大神及生信技能树
  • TCGA数据库的样本编码规律:Tumor types range from 01 - 09, normal types from 10 - 19 and control samples from 20 - 29,方法对样本进行分组。
  • 基因注释GMT文件把mRNA矩阵注释拆分成了coding RNA和lncRNA表达矩阵。
  • DESeq2包分析的表达矩阵必须是count矩阵,TCGA等数据库下载的表达矩阵需查看格式说明,如是否取log,如有需要转换回来。
  • 包装函数for循环方便根据临床表型信息批量筛选表达矩阵,绘制热图、PCA图、差异分析并得到火山图、venn图。
  • 模仿文献分析方法挖掘数据需要仔细阅读文献,查看表达矩阵的过滤条件和差异分析阈值(FC和log2FC有区别)。

函数代码

检查数据函数
run_check_h_pca <- function(pro = "T1_vs_Normal"){
  rm(list = ls())  ## 魔幻操作,一键清空~
  options(stringsAsFactors = F)
  table(group_list)
  # 每次都要检测数据
  dat[1:4,1:4]

  ## 下面是画PCA的必须操作,需要看说明书。
  dat=t(dat)#画PCA图时要求是行名时样本名,列名时探针名,因此此时需要转换
  dat=as.data.frame(dat)#将matrix转换为data.frame
  dat=cbind(dat,group_list) #cbind横向追加,即将分组信息追加到最后一列
  library("FactoMineR")#画主成分分析图需要加载这两个包
  library("factoextra") 
  # The variable group_list (index = 54676) is removed
  # before PCA analysis
  dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#现在dat最后一列是group_list,需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
  fviz_pca_ind(dat.pca,
               geom.ind = "point", # show points only (nbut not "text")
               col.ind = dat$group_list, # color by groups
               # palette = c("#00AFBB", "#E7B800"),
               addEllipses = TRUE, # Concentration ellipses
               legend.title = "Groups"
  )
  ggsave(filename = paste0(pro,'_all_samples_PCA.png'))

  rm(list = ls())  ## 魔幻操作,一键清空~
  load(file = 'step1_output.Rdata') 
  dat[1:4,1:4] 

  cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行('1'是按行取,'2'是按列取)取每一行的方差,从小到大排序,取最大的1000个
  library(pheatmap)
  #pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出,组合起来为一个新的表达矩阵
  n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
  n[n>2]=2 
  n[n< -2]= -2
  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(g=group_list)
  rownames(ac)=colnames(n) #把ac的行名给到n的列名,即对每一个探针标记上分组信息(是'noTNBC'还是'TNBC')
  ## 可以看到TNBC具有一定的异质性,拿它来区分乳腺癌亚型指导临床治疗还是略显粗糙。
  pheatmap(n,show_colnames =F,show_rownames = F,
           annotation_col=ac,filename = paste0(pro,"_heatmap_top1000_sd.png"))


  rm(list = ls()) 
  load(file = 'step1_output.Rdata') 
  dat[1:4,1:4] 
  exprSet=dat
  pheatmap::pheatmap(cor(exprSet)) 
  # 组内的样本的相似性应该是要高于组间的!
  colD=data.frame(group_list=group_list)
  rownames(colD)=colnames(exprSet)
  pheatmap::pheatmap(cor(exprSet),
                     annotation_col = colD,
                     show_rownames = F,
                     filename = paste0(pro,"_cor_all.png"))
  dim(exprSet)
  exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
  dim(exprSet)

  exprSet=log(edgeR::cpm(exprSet)+1)
  dim(exprSet)
  exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
  dim(exprSet)
  M=cor(log2(exprSet+1)) 
  pheatmap::pheatmap(M,
                     show_rownames = F,
                     annotation_col = colD,
                     filename = paste0(pro,"_cor_top500.png"))
}
DESeq2差异分析函数
run_filter_check_DESeq2 <- function(Typel,pro){
  #按照T1-T4分类肺癌组织样本的得到四个表达矩阵
  RNA_Normal=gset[,substring(colnames(gset),14,15) == "11"]
  RNA_LUSC=gset[,substring(colnames(gset),14,15) == "01"]

  colnames(RNA_Normal)
  colnames(RNA_LUSC)

  # T1-4期亚型的表达矩阵
  RNA_T=RNA_LUSC[,colnames(RNA_LUSC) %in% rownames(pd[substring(pd[,"pathologic_T"],1,2) == Typel,])]
  dat=cbind(RNA_Normal,RNA_T)
  colnames(dat)
  group_list=c(rep("Normal",38),rep(Typel,ncol(RNA_T)))
  table(group_list)

  ## 差异分析 DESeq2

  if(!require(DESeq2))BiocManager::install("DESeq2")
  library(DESeq2)
  #需修改results()的contrast参数
  #输入:表达矩阵和分组信息
  #输出:差异分析结果、火山图
  #构建colData (condition存在于colData中,是表示分组的因子型变量)
  countData <- floor(dat)
  colData <- data.frame(row.names =colnames(countData), 
                        condition=group_list)
  dds <- DESeqDataSetFromMatrix(
    countData = countData,
    colData = colData,
    design = ~ condition)

  dds <- DESeq(dds)
  # 两两比较
  res <- results(dds, contrast = c("condition",Typel,"Normal"))
  resOrdered <- res[order(res$padj),] # 按照P值排序
  DEG <- as.data.frame(resOrdered)
  DEG <- na.omit(DEG)

  #添加change列标记基因上调下调,在DESeq里FDR就是pdaj,所以要把pvalue改成pdaj
  #logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)) )
  logFC_cutoff <- 1
  DEG$change = as.factor(
    ifelse(DEG$padj < 0.01 & abs(DEG$log2FoldChange) > logFC_cutoff,
           ifelse(DEG$log2FoldChange >= logFC_cutoff ,'UP','DOWN'),'NOT')
  )
  head(DEG)
  boxplot(as.numeric(countData[rownames(DEG)[1],])~group_list)
  # Heatmap热图
  choose_gene=head(rownames(DEG),100) ## 选定部分差异基因
  choose_matrix=countData[choose_gene,]
  pheatmap::pheatmap(choose_matrix) 

  annotation_col = data.frame(
    group = factor(group_list))
  rownames(annotation_col) = colnames(countData)

  pheatmap::pheatmap(choose_matrix,
                     scale = "row",
                     annotation_col = annotation_col)

  # 火山图
  library(ggplot2)

  this_tile <- paste0('Cutoff for log2FoldChange is ',round(logFC_cutoff,3),
                      '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                      '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
  )
  p1 = ggplot(data=DEG, 
              aes(x=log2FoldChange, y=-log10(padj),color=change)) +
    geom_point(alpha=0.4, size=3.5) +
    theme_set(theme_set(theme_bw(base_size=20)))+
    xlab("log2 fold change") + ylab("-log10 padj") +
    geom_vline(xintercept=c(-logFC_cutoff,logFC_cutoff),lty=4,col="black",lwd=0.8) +
    geom_hline(yintercept = -log10(0.05),lty=4,col="black",lwd=0.8) +
    ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
    scale_colour_manual(values = c('blue','grey','red')) ## corresponding to the levels(res$change)
  p1
  ggsave(p1, filename = paste0(pro,"deseq2_vocalno.png"), width = 10, height = 7 )
  save(DEG,p1,file = paste0(pro,"_deseq2_DEG.Rdata"))

}文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI,多一点数据认知,让他们的科研上一个台阶:

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-03-09

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 需要5个步骤来说明你想研究的基因的重要性

    根据你生物学故事的大小来说,有不同的论证方法,当然,普罗大众喜闻乐见的当然是公共数据库网页工具不停的点点点,就出来了一堆图表证明自己的基因有研究的意义。证据链有...

    生信技能树
  • 构建shell脚本一文就够

    非常多的朋友在看我们公众号过往转录组,WES,等流程分享的时候发现很难理解我们的代码,其实就是缺乏shell脚本知识,那么这篇教程你就不容错过。 内容 使用多个...

    生信技能树
  • 每月一生信流程之rnaseqGene

    但是RNA-seq的分析肯定远不止那些啦,拿到基于基因的表达矩阵固然可以根据转录组经典表达量矩阵下游分析大全 里面的R包和代码进行统计可视化,但是表达矩阵并不是...

    生信技能树
  • Hades开源白盒审计系统V1.0.0

    为什么要开发这么一个白盒审计系统呢?其实,自己开发白盒审计系统的想法已经在我的脑海中存在很久了。一方面,之前大学毕业之后在白盒方面花了挺多时间的,不过由于各方面...

    周俊辉
  • 各种基本算法实现小结(四)—— 图及其遍历

    ====================================================================

    阳光岛主
  • 基于 HTML5 WebGL 的楼宇智能化集成系统(三)

    2018年7月,信息化部印发了《工业互联网平台建设及推广指南》和《工业互联网平台评价方法》,掀起了 工业互联网 的浪潮,并成为热词写入了报告中。同为...

    HT for Web
  • 深入正则表达式(0):正则表达式概述

    正则表达式(regular expression,在代码中常简写为regex、regexp或RE),又称正规表示式、正規表示法、正規運算式、規則運算式、常規表示...

    周陆军
  • Java 8中集合优雅快速的处理方式

    相信现在大多数的伙伴们,都在使用Java 8了,而 Java 8相比以前的版本,是作出了革命性的改变。Java8的特性大致可总结为,开发速度更快,代码更少,增加...

    攻城狮的那点事
  • leetcode381. Insert Delete GetRandom O(1) - Duplicates allowed

    设计一个数据结构,支持能够在O(1)的时间内完成对数字的插入,删除和获取随机数的操作,允许插入重复的数字,同时要求每个数字被随机获取的概率和该数字当前在数据结构...

    眯眯眼的猫头鹰
  • 框架篇-Django博客应用-权限控制

    然后让 BlogPublishView 和 BlogEditView 类继承 AdminRequiredMixin即可。

    小团子

扫码关注云+社区

领取腾讯云代金券