• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • (未测试)实战:TCGA数据差异分析三大R包及其结果对比

      (未测试)实战:TCGA数据差异分析三大R包及其结果对比

      • 发布者 一览
      • 分类 未分类
      • 日期 2020年2月6日
      • 评论 0评论

      参考网址: https://www.jianshu.com/p/b8a124501073

      原本还有第四个部分,小洁老师讲了另一个R包下载表达矩阵和临床信息的,
      TCGA-4.使用RTCGA包获取数据
      但是这个包有个缺点就是数据更新不及时,因此当时看到时候我就没有跟学了。直接跳到第五步TCGA-5.(转录组)差异分析三大R包及其结果对比
      但是呢,由于没跟学第四步这一步获取数据并做数据清洗的时候出了问题,一直没能完成,后来昨天花了点时间学了冰糖在菜鸟团的推文也是小洁老师的第四步教程相关的内容,对比来看一步步调试,再加上从技能树推文得到的小洁老师的画图函数后,终于完成了第五步的学习。
      还是很有收获的。

      1.提前准备安装和加载R包

      rm(list = ls())
      options(stringsAsFactors = F)
      if(!require(stringr))install.packages('stringr')
      if(!require(ggplotify))install.packages("ggplotify")
      if(!require(patchwork))install.packages("patchwork")
      if(!require(cowplot))install.packages("cowplot")
      if(!require(DESeq2))install.packages('DESeq2')
      if(!require(edgeR))install.packages('edgeR')
      if(!require(limma))install.packages('limma')

      2.准备数据

      本示例的数据是TCGA-KIRC的表达矩阵。tcga样本编号14-15位是隐藏分组信息的,详见:
      TCGA的样本id里藏着分组信息

      TCGA样本id,分组信息是在这个id的第14-15位,01-09是tumor,10-29是normal。

      #TCGA-KIRC
      library(TCGAbiolinks)
      #可以查看所有支持的癌症种类的缩写
      #TCGAbiolinks:::getGDCprojects()$project_id
      #还是选择之前的例子
      cancer_type="TCGA-KIRC"
      clinical <- GDCquery_clinic(project = cancer_type, type = "clinical")
      clinical[1:4,1:4]
      dim(clinical)
      
      query <- GDCquery(project = cancer_type, 
                        data.category = "Transcriptome Profiling", 
                        data.type = "miRNA Expression Quantification", 
                        workflow.type = "BCGSC miRNA Profiling")
      GDCdownload(query, method = "api", files.per.chunk = 50)
      expdat <- GDCprepare(query = query)
      expdat[1:3,1:3]
      library(tibble)
      rownames(expdat) <- NULL
      expdat <- column_to_rownames(expdat,var = "miRNA_ID")
      expdat[1:3,1:3]
      exp = t(expdat[,seq(1,ncol(expdat),3)])
      exp[1:4,1:4]
      expr=exp
      rowName <- str_split(rownames(exp),'_',simplify = T)[,3]
      expr<- apply(expr,2,as.numeric) 
      expr<- na.omit(expr)
      dim(expr)
      expr <- expr[,apply(expr, 2,function(x){sum(x>1)>10})]
      rownames(expr) <- rowName
      dim(expr)
      expr[1:4,1:4]
      save(expr,clinical,file = "tcga-kirc-download.Rdata")
      rm(list = ls())
      load("tcga-kirc-download.Rdata") #获取初步下载数据。
      meta <- clinical
      colnames(meta)
      meta <- meta[,c("submitter_id","vital_status",
                      "days_to_death","days_to_last_follow_up",
                      "race",
                      "age_at_diagnosis",
                      "gender" ,
                      "ajcc_pathologic_stage")]
      expr=t(expr)
      expr[1:4,1:4]
      group_list <- ifelse(as.numeric(str_sub(colnames(expr),14,15))<10,"tumor","normal")
      group_list <- factor(group_list,levels = c("normal","tumor"))
      
      table(group_list)
      # normal  tumor 
      # 71    545
      save(expr,group_list,file = "tcga-kirc-raw.Rdata")

      由于不知道小洁老师做了怎样的过滤,我得到的结果不同
      我觉得应该是在mata这个代码步骤后面选择一个指标过滤掉一些数据。
      先放着,这个代码在这个步骤中没有用到。以后应该会用到。
      由于不会自己写代码,后面的分析基本上就是走的小洁老师教程的内容。

      3.三大R包的差异分析

      #Deseq2
      library(DESeq2)
      colData <- data.frame(row.names =colnames(expr), 
                            condition=group_list)
      dds <- DESeqDataSetFromMatrix(
        countData = expr,
        colData = colData,
        design = ~ condition)
      #参考因子应该是对照组 dds$condition <- relevel(dds$condition, ref = "untrt")
      
      dds <- DESeq(dds)
      # 两两比较
      res <- results(dds, contrast = c("condition",rev(levels(group_list))))
      resOrdered <- res[order(res$pvalue),] # 按照P值排序
      DEG <- as.data.frame(resOrdered)
      head(DEG)
      # 去除NA值
      DEG <- na.omit(DEG)
      
      #添加change列标记基因上调下调
      #logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)) )
      logFC_cutoff <- 1
      DEG$change = as.factor(
        ifelse(DEG$pvalue < 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
               ifelse(DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')
      )
      head(DEG)
      
      DESeq2_DEG <- DEG
      
      #edgeR
      library(edgeR)
      
      dge <- DGEList(counts=expr,group=group_list)
      dge$samples$lib.size <- colSums(dge$counts)
      dge <- calcNormFactors(dge) 
      
      design <- model.matrix(~0+group_list)
      rownames(design)<-colnames(dge)
      colnames(design)<-levels(group_list)
      
      dge <- estimateGLMCommonDisp(dge,design)
      dge <- estimateGLMTrendedDisp(dge, design)
      dge <- estimateGLMTagwiseDisp(dge, design)
      
      fit <- glmFit(dge, design)
      fit2 <- glmLRT(fit, contrast=c(-1,1)) 
      
      DEG=topTags(fit2, n=nrow(expr))
      DEG=as.data.frame(DEG)
      logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )
      logFC_cutoff <- 1
      DEG$change = as.factor(
        ifelse(DEG$PValue < 0.05 & abs(DEG$logFC) > logFC_cutoff,
               ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
      )
      head(DEG)
      table(DEG$change)
      edgeR_DEG <- DEG
      
      #limma-voom
      library(limma)
      
      design <- model.matrix(~0+group_list)
      colnames(design)=levels(group_list)
      rownames(design)=colnames(expr)
      
      dge <- DGEList(counts=expr)
      dge <- calcNormFactors(dge)
      logCPM <- cpm(dge, log=TRUE, prior.count=3)
      
      v <- voom(dge,design, normalize="quantile")
      fit <- lmFit(v, design)
      
      constrasts = paste(rev(levels(group_list)),collapse = "-")
      cont.matrix <- makeContrasts(contrasts=constrasts,levels = design) 
      fit2=contrasts.fit(fit,cont.matrix)
      fit2=eBayes(fit2)
      
      DEG = topTable(fit2, coef=constrasts, n=Inf)
      DEG = na.omit(DEG)
      #logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )
      logFC_cutoff <- 1
      DEG$change = as.factor(
        ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
               ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
      )
      head(DEG)
      limma_voom_DEG <- DEG
      save(DESeq2_DEG,edgeR_DEG,limma_voom_DEG,group_list,file = "DEG.Rdata")
      
      #差异分析结果的可视化
      rm(list = ls())
      load("tcga-kirc-raw.Rdata")
      load("DEG.Rdata")
      source("3-plotfunction.R")
      logFC_cutoff <- 1
      expr[1:4,1:4]
      dat = log(expr+1)
      pca.plot = draw_pca(dat,group_list)
      
      cg1 = rownames(DESeq2_DEG)[DESeq2_DEG$change !="NOT"]
      cg2 = rownames(edgeR_DEG)[edgeR_DEG$change !="NOT"]
      cg3 = rownames(limma_voom_DEG)[limma_voom_DEG$change !="NOT"]
      
      h1 = draw_heatmap(expr[cg1,],group_list)
      h2 = draw_heatmap(expr[cg2,],group_list)
      h3 = draw_heatmap(expr[cg3,],group_list)
      
      v1 = draw_volcano(test = DESeq2_DEG[,c(2,5,7)],pkg = 1)
      v2 = draw_volcano(test = edgeR_DEG[,c(1,4,6)],pkg = 2)
      v3 = draw_volcano(test = limma_voom_DEG[,c(1,4,7)],pkg = 3)
      
      library(patchwork)
      (h1 + h2 + h3) / (v1 + v2 + v3) +plot_layout(guides = 'collect')
      
      #(v1 + v2 + v3) +plot_layout(guides = 'collect')
      ggsave("heat_volcano.png",width = 21,height = 9)
      #三大R包差异基因对比
      # 三大R包差异基因交集
      UP=function(df){
        rownames(df)[df$change=="UP"]
      }
      DOWN=function(df){
        rownames(df)[df$change=="DOWN"]
      }
      
      up = intersect(intersect(UP(DESeq2_DEG),UP(edgeR_DEG)),UP(limma_voom_DEG))
      down = intersect(intersect(DOWN(DESeq2_DEG),DOWN(edgeR_DEG)),DOWN(limma_voom_DEG))
      
      hp = draw_heatmap(expr[c(up,down),],group_list)
      
      #上调、下调基因分别画维恩图
      
      up.plot <- venn(UP(DESeq2_DEG),UP(edgeR_DEG),UP(limma_voom_DEG),
                      "UPgene"
      )
      down.plot <- venn(DOWN(DESeq2_DEG),DOWN(edgeR_DEG),DOWN(limma_voom_DEG),
                        "DOWNgene"
      )
      
      library(cowplot)
      library(ggplotify)
      up.plot = as.ggplot(as_grob(up.plot))
      down.plot = as.ggplot(as_grob(down.plot))
      library(patchwork)
      #up.plot + down.plot
      
      pca.plot + hp+up.plot +down.plot
      ggsave("deg.png",height = 10,width = 10)

      整个流程走完得到的结果如下:





      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      作者头像
      一览

      上一篇文章

      (未测试)一文重复一篇四分SCI------基于TCGA和geo的circRNA研究(1)
      2020年2月6日

      下一篇文章

      (未测试)R语言之生信,差异基因分析第一节
      2020年2月6日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?