专栏首页小明的数据分析笔记本使用R语言下载TCGA数据库癌症基因表达数据小例子

使用R语言下载TCGA数据库癌症基因表达数据小例子

参考资料
  • 生信技能树 公众号文章 TCGA数据下载—TCGAbiolinks包参数详解
  • 生信技能树 公众号文章 批量COX回归生存分析图,指定挑选lncRNA基因,森林图,ROC曲线打包给你
  • 生信星球 公众号文章 重磅!TCGA数据分析流程梳理总结
  • 生信星球 公众号文章 TCGA3.R包TCGAbiolinks下载数据
  • 生信星球 公众号文章 TCGA的样本id里藏着分组信息
  • 简书文章 TCGA癌症缩写、癌症中英文对照
  • CDSN文章 通过整理TCGA数据,探索某癌症的癌组织和正常组织的差异基因
  • TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages
  • TCGAbiolinks包下载TCGA数据进行表达差异分析-乳腺癌案例
代码

数据下载

BiocManager::install("TCGAbiolinks")
library(TCGAbiolinks)
?GDCquery
x<-getGDCprojects()$project_id
x == "TCGA-COAD"
x == "TCGA-READ"
TCGAbiolinks:::getProjectSummary("TCGA-COAD")
query<-GDCquery(project = "TCGA-COAD",
                legacy=F,
                experimental.strategy ="RNA-Seq",
                data.category = "Transcriptome Profiling",
                data.type="Gene Expression Quantification",
                workflow.type="HTSeq - Counts")
GDCdownload(query)

这里遇到的问题是:所有数据都下载下来了,如何将数据整合成表达量矩阵的形式呢?

更新

上面遇到的问题今天找到了解决办法:可以使用SummarizedExperiment包中的函数assay()函数将表达矩阵提取出来;colData()函数好像是获得一些样本信息

查看TCGAbiolinks包的帮助文档 browseVignettes("TCGAbiolinks") 尝试重复帮助文档中的第3项内容 Downloading and preparing files for analysis 这篇帮助文档中用到的例子是TCGA-GBM——Glioblastoma multiforme——多形成性胶质细胞瘤

获取基因表达矩阵
query_3 <- GDCquery(project = "TCGA-GBM",
                    data.category = "Gene expression",
                    data.type = "Gene expression quantification",
                    platform = "Illumina HiSeq", 
                    file.type  = "results",
                    experimental.strategy = "RNA-Seq",
                    barcode = barcode,
                    legacy = TRUE)
GDCdownload(query_3)
dataPrep1<-GDCprepare(query=query_3)
library(SummarizedExperiment)
expdat<-assay(dataPrep1)
dim(expdat)
head(expdat)

备注:file.type 的另外一个参数是 normalized_results

DESeq2基因差异表达分析
library(stringr)
group_list<-ifelse(str_sub(colnames(expdat),14,15) == "01","tumor","normal")
group_list
coldata<-data.frame(row.names = colnames(expdat),
                    condition=group_list)
coldata
library(DESeq2)
expdat<-round(expdat,0)
dds<-DESeqDataSetFromMatrix(countData = expdat,
                            colData = coldata,
                            design = ~ condition)
keep<-rowSums(counts(dds))>=10
dds<-dds[keep,]
dds<-DESeq(dds)
res<-results(dds,contrast = c("condition","tumor","normal"))
resOrdered <- res[order(res$pvalue),]
DEG <- as.data.frame(resOrdered)
DEG <- na.omit(DEG)
logFC_cutoff<-with(DEG,mean(abs(log2FoldChange))+2*sd(abs(log2FoldChange)))
logFC_cutoff
DEG$change<-as.factor(ifelse(DEG$pvalue<0.05&abs(DEG$log2FoldChange)>logFC_cutoff,
                             ifelse(DEG$log2FoldChange>logFC_cutoff,"UP","DOWN"),
                             "NOT"))
library(ggplot2)
this_title <- paste0('Cutoff for logFC 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',]))
ggplot(data=DEG,aes(x=log2FoldChange,y=-log10(pvalue),color=change))+
  geom_point(alpha=0.4,size=1.75)+
  labs(x="log2 fold change")+ ylab("-log10 pvalue")+
  ggtitle(this_title)+theme_bw(base_size = 10)+
  theme(plot.title = element_text(size=15,hjust=0.5))+
  scale_color_manual(values=c('blue','black','red'))

Rplot.png

使用R语言包 clusterProfiler 差异表达基因的GO富集分析

先看一下这个包的帮助文档

browseVignettes("clusterProfiler")
help(package="clusterProfiler")
library(clusterProfiler)
data(geneList,package="DOSE")
names(geneList)[1:100]

准备自己的输入数据

df<-DEG[which(DEG$change=="UP"),]
dim(df)
rownames(df)
gene_list<-str_split(rownames(df),"\\|")
gene_list
class(gene_list)
gene_list<-as.data.frame(matrix(unlist(a),ncol=2,byrow = T))
head(gene_list)
yy<-enrichGO(gene=gene_list$V2,
             OrgDb='org.Hs.eg.db',
             ont="BP",
             pvalueCutoff = 0.01)
dotplot(yy)

Rplot01.png

有时间试着用今天分析的数据准备一下GOplot包的输入数据。

本文分享自微信公众号 - 小明的数据分析笔记本(gh_0c8895f349d3),作者:Punicagranatum

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

原始发表时间:2020-02-18

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R语言画地图的几篇文章记录

    ggplot2添加箭头https://ggplot2.tidyverse.org/reference/geom_segment.html

    用户7010445
  • ggplot2:在一幅图中插入另外一幅图

    用户7010445
  • 生存分析中给基因表达量(连续变量)设置阈值

    在论文 Construction of a Competitive endogenous RNA network and identification of p...

    用户7010445
  • linux系统环境变量一文就够

    Linux是一个多用户的操作系统。每个用户登录系统后,都会有一个专用的运行环境。 通常每个用户默认的环境都是相同的,这个默认环境实际上就是一组环境变量的定义。 ...

    生信技能树
  • [译]RabbitMQ教程C#版 - 路由

    在本教程中,我们会日志系统其再添加一个特性,使其可以只订阅消息的一个子集。例如,将所有日志消息打印到 控制台的同时,只会将严重错误消息写入日志文件(保存到磁盘...

    Esofar
  • Python 中gzip模块完成对文件的

    gzip块主要支持打开对应格式的压缩文件,并可以完成对压缩文件的读出和写入操作。压缩文件被打开后,可以使用文件对象一样的方法,如read、readline、re...

    py3study
  • ASP.NET Core中的ActionFilter与DI

      前几篇文章都是讲ASP.NET Core MVC中的依赖注入(DI)与扩展点的,也许大家都发现在ASP.NET CORE中所有的组件都是通过依赖注入来扩展的...

    yoyofx
  • unicode中的‘\xa0’字符在转换成gbk编码时会出现问题,gbk无法转换'\xa0'字符。

    unicode中的‘\xa0’字符在转换成gbk编码时会出现问题,gbk无法转换’\xa0’字符。 所以,在转换的时候必需进行一些前置动作:

    学到老
  • notification通知栏

    https://blog.csdn.net/fzkf9225/article/details/81119780

    平凡的学生族
  • Python 协程检测Kubernetes服务端口

    https://www.cnblogs.com/xiao987334176/p/10237551.html

    py3study

扫码关注云+社区

领取腾讯云代金券