• 主页
  • 课程

    关于课程

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

    同等学历教学

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

      关于课程

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

      同等学历教学

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

      未分类

      • 首页
      • 博客
      • 未分类
      • 已测试:TCGA差异表达分析和生存分析

      已测试:TCGA差异表达分析和生存分析

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

      TCGA数据库介绍请自行查阅资料,推荐链接:https://www.jianshu.com/p/dd066035b399

      上图可以看到,TCGA的数据包括了基因测序,基因拷贝数分析,甲基化分析,转录组分析,miRNA分析,临床数据。

      那么,我们可以分析哪些数据呢?

      基因测序数据,我们可以看出基因突变;基因拷贝数,我们可以看到基因组中哪些基因拷贝数多;甲基化数据可以找出差异的甲基化位点或者区域;转录组可以知道哪些基因高表达,哪些低表达;miRNA分析,可以找出差异的miRNA,或者找出跟癌症相关的miRNA。它还有临床的数据,基因的数据可以和临床数据结合,比如不同的分组,是否突变不同,对生存期是否不同。想知道某个基因的表达,对病人的预后是否影响,可以把转录组和预后进行分析。分子水平也可以联合分析,如转录组可以和miRNA联合起来,可以看出哪些miRNA会和基因的表达有关,哪个miRNA会对某个基因有调控的作用。甲基化也可以和转录组联合,甲基化的水平对基因表达的影响。

      上图是样品的名称,sample中,01-09就是肿瘤组织,10-19是正常组织,20-29是正常人。

      进入TCGA官网

      https://portal.gdc.cancer.gov/

      左上角的选择case, 可以选择不同组织,根据左侧列表进行筛选。我们选择前列腺癌,如下图,右边显示的是2326个数据。

      下面program里面选择TCGA。如果我们要选择多个数据库,也可以点选。

      case里面选择完了,我们继续选择file里面的数据。

      我们选择转录本的数据。

      加入购物车以后,点击购物车,则可以看到如下可以下载的页面,我们可以点击各个下载的地方,如download里面,有一个manifest文件,是需要直接下载的,此文件是解释文件,解释我们即将下载的文件。cart是基因表达量的文件。clinical是临床数据。biospecimen人群数据。filemetadata是名单库,需要和我们的下载的名单对应起来。

      以上文件依次下载,但是cart文件可能经常下载不了,网络问题等导致下载失败,我们就需要用到命令行下载。

      需要先进行工具下载,TCGA主页右上角有一个GDC apps,下载data transfer tool。根据自己的系统下载,不需要安装,因为软件是二进制的。

      找到manifest文件,放在同一个目录下,打开dos环境,运行gdc-client.exe download -m MANIFEST.txt

      一定下载最新版本的GDC工具,不然会报错,TCGA经常更新。

      上图显示正在下载。

      下一步则进行数据整理,整理成矩阵以后,则进行id转为gene symbol,整理生存数据。

      把所有的文件都放在一个目录下。这个可以手动,也可以通过代码。

      用perl移动文件至同一个目录,这里要注意,不能把miRNA和mRNA的数据下载到同一个目录下,不然都给融合到一个文件夹了。具体的文件可下载附件,或者留言索取。如下图所示。

      进入cmd,进入目录文件夹,输入以下命令,回车,就可以了。

      C:\Users\Administrator\Desktop\111>perl putFilesToOneDir.pl

      就会出现以下一个文件夹,就是完全整理到一个文件夹了。

      然后全部解压。

      然后把mRNA的数据整合到一起,用mRNAmerge脚本和metadata文件,前面已经下载过,复制到上面解压后文件的文件夹。同样按照perl的操作方式,进入到目录以后,运行perl,空格,后面跟着mRNAmerge脚本,再加上metadata文件。

      C:\Users\Administrator\Desktop\111\files>perl mRNA_merge.pl metadata.cart.2020-02-20.json

      回车后,这个运行需要一定的时间。

      运行完了以后,系统会告诉你,normal有多少,tumor有多少。这个数字要记录好,后面会用得到。

      运行完毕,会出现一个mRNAmatrix.txt文件,就是我们后续分析需要用到的。

      打开这个文件,我们可以看到,前51个就是normal。关键的字段是11,而不是1到9。

      下一步,我们要把ensemble ID转化为基因名。

      对应文件需要到ensemble的官网去下载。http://asia.ensembl.org/index.html

      download,下载。

      下载好以后解压,然后运行perl ensemblToSymbol.pl Homo_sapiens.GRCh38.99.chr.gtf mRNAmatrix.txt mRNA.symbol.txt(这里是我们自己命名的输出文件名字。)

      运行后,则基因的名字转换完毕。

      foldChange=1   #设置差异倍数
      padj=0.05   #设置P值
      library("edgeR")
      rt=read.table("mRNA.symbol.txt",sep="\t",header=T,check.names=F)  #改成自己的文件名
      rt=as.matrix(rt)
      rownames(rt)=rt[,1]
      exp=rt[,2:ncol(rt)]
      dimnames=list(rownames(exp),colnames(exp))
      data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
      data=avereps(data)
      data=data[rowMeans(data)>1,]
      group=c(rep("normal",51),rep("tumor",489))                         #按照癌症和正常样品数目修改
      design <- model.matrix(~group)
      y <- DGEList(counts=data,group=group)
      y <- calcNormFactors(y)
      y <- estimateCommonDisp(y)
      y <- estimateTagwiseDisp(y)
      et <- exactTest(y,pair = c("normal","tumor"))
      topTags(et)
      ordered_tags <- topTags(et, n=100000)
      allDiff=ordered_tags$table
      allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]
      diff=allDiff
      newData=y$pseudo.counts   #差异分析
      
      write.table(diff,file="edgerOut.xls",sep="\t",quote=F)
      diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]
      write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)
      diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]
      write.table(diffUp, file="up.xls",sep="\t",quote=F)
      diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]
      write.table(diffDown, file="down.xls",sep="\t",quote=F)
      
      normalizeExp=rbind(id=colnames(newData),newData)
      write.table(normalizeExp,file="normalizeExp.txt",sep="\t",quote=F,col.names=F)   #输出所有基因校正后的表达值(normalizeExp.txt)
      diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),])
      write.table(diffExp,file="diffmRNAExp.txt",sep="\t",quote=F,col.names=F)         #输出差异基因校正后的表达值(diffmRNAExp.txt)
      
      heatmapData <- newData[rownames(diffSig),]
      #volcano火山图
      pdf(file="vol.pdf")
      xMax=max(-log10(allDiff$FDR))+1 #这一步很关键,如果allDiff里面有0,那xMax是无穷大,所以会报错,让你输入有限数值的定义,这里我们可以看一下我们的数据里面,FDR的最小值是多少,计算出来,直接输入该值,即xMax=300,就像下面的yMax。
      yMax=12
      plot(-log10(allDiff$FDR), allDiff$logFC, xlab="-log10(FDR)",ylab="logFC",
           main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)
      diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC>foldChange,]
      points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="red",cex=0.4)
      diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC<(-foldChange),]
      points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="green",cex=0.4)
      abline(h=0,lty=2,lwd=3)
      dev.off()
      #heatmap热图
      hmExp=log10(heatmapData+0.001)
      library('gplots')
      hmMat=as.matrix(hmExp)
      pdf(file="heatmap.pdf",width=60,height=90)
      par(oma=c(10,3,3,7))
      heatmap.2(hmMat,col='greenred',trace="none")
      dev.off()

      结果如下,我们产生了相关文件,这些文件都是在后续画图中有用的文件。

      热图和火山图如下

      下面我们进一步下载miRNA的数据,与上述不同的是,这次我们选择miRNAseq,而不是mRNAseq。分析完了以后看看关联情况。与上述下载方式类似。

      把所有需要的文件下载下来,进一步进行分析。需要用到menifest文件。

      cart文件需要代码下载。如上图所示

      那么下载下来以后,我们继续把miRNA的文件融合到一起,这里不需要放到同一个文件夹,因为现在已经是txt文件,直接用perl脚本整合到一个文件夹。整合完了以后,已经是标准的miRNA ID,不再需要转换。下面进行差异分析和热图火山图的绘制。融合的文件如下,miRNAmatrix.txt,内容如下图。

      上图显示,normal52,后面分组需要用到,请记录。

      foldChange=1
      padj=0.05
      
      library("edgeR")
      rt=read.table("miRNAmatrix.txt",sep="\t",header=T,check.names=F)  #改成自己的文件名
      rt=as.matrix(rt)
      rownames(rt)=rt[,1]
      exp=rt[,2:ncol(rt)]
      dimnames=list(rownames(exp),colnames(exp))
      data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
      data=avereps(data)
      data=data[rowMeans(data)>0,]
      
      #group=c("normal","tumor","tumor","normal","tumor")
      group=c(rep("normal",4),rep("tumor",179))                         #按照癌症和正常样品数目修改
      design <- model.matrix(~group)
      y <- DGEList(counts=data,group=group)
      y <- calcNormFactors(y)
      y <- estimateCommonDisp(y)
      y <- estimateTagwiseDisp(y)
      et <- exactTest(y,pair = c("normal","tumor"))
      topTags(et)
      ordered_tags <- topTags(et, n=100000)
      
      allDiff=ordered_tags$table
      allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]
      diff=allDiff
      newData=y$pseudo.counts
      
      write.table(diff,file="edgerOut.xls",sep="\t",quote=F)
      diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]
      write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)
      diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]
      write.table(diffUp, file="up.xls",sep="\t",quote=F)
      diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]
      write.table(diffDown, file="down.xls",sep="\t",quote=F)
      
      normalizeExp=rbind(id=colnames(newData),newData)
      write.table(normalizeExp,file="normalizeExp.txt",sep="\t",quote=F,col.names=F)   #输出所有miRNA校正后的表达值(normalizeExp.txt)
      diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),])
      write.table(diffExp,file="diffmiRNAExp.txt",sep="\t",quote=F,col.names=F)         #输出差异miRNA校正后的表达值(diffmiRNAExp.txt)
      
      heatmapData <- newData[rownames(diffSig),]
      
      #volcano
      pdf(file="vol.pdf")
      xMax=max(-log10(allDiff$FDR))+1
      yMax=12
      plot(-log10(allDiff$FDR), allDiff$logFC, xlab="-log10(FDR)",ylab="logFC",
           main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)
      diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC>foldChange,]
      points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="red",cex=0.4)
      diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC<(-foldChange),]
      points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="green",cex=0.4)
      abline(h=0,lty=2,lwd=3)
      dev.off()
      
      #heatmap
      hmExp=log10(heatmapData+0.001)
      library('gplots')
      hmMat=as.matrix(hmExp)
      pdf(file="heatmap.pdf",width=60,height=30)
      par(oma=c(10,3,3,7))
      heatmap.2(hmMat,col='greenred',trace="none")
      dev.off()

      上图是miRNA的差异结果图。

      下面,我们要进一步分析差异的mRNA和差异的miRNA有什么联系呢?

      我们把diffmiRNAExp.txt和 diffmRNAExp.txt 放到目录文件夹,用perl脚本将其融合到一个文件。

      miRNA = read.table("diffmiRNAExp.txt", row.names=1 ,header=T,sep="\t",check.names=F)
      RNA = read.table("diffmRNAExp.txt", row.names=1 ,header=T,sep="\t",check.names=F)
      colnames(miRNA)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3\\-\\4",colnames(miRNA))
      colnames(RNA)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3\\-\\4",colnames(RNA))
      sameSample=intersect(colnames(miRNA),colnames(RNA))
      merge=rbind(id=sameSample,miRNA[,sameSample],RNA[,sameSample])
      write.table(merge,file="merge.txt",sep="\t",quote=F,col.names=F)
      
      miRNALabel=cbind(rownames(miRNA),"miRNA")
      RNALabel=cbind(rownames(RNA),"gene")
      nodeLabel=rbind(c("ID","Classify"),miRNALabel,RNALabel)
      write.table(nodeLabel,file="nodeLabel.txt",sep="\t",quote=F,col.names=F,row.names=F)

      nodelabel就是指定哪些属于miRNA,哪些属于mRNA,给一个标签,后续用cytoscape画图时候用,保存好。

      然后下一步我们用WGCNA去分析共表达情况。

      待续。。。。。。

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

      • 分享:
      作者头像
      一览

      上一篇文章

      懒花花想要博客,于是我就…
      2020年4月8日

      下一篇文章

      这个人骑车骑出了一个向量,然后。。。
      2020年4月9日

      你可能也喜欢

      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年
      在线支付 激活码

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