单细胞 GEO TCGA验证 实验验证是目前生信的新趋势。我们相信，随着单细胞测序的普及，生信分析和数据挖掘在科研中扮演越来越重要的角色。因此，R 语言已经成为每个科研人的必备技能。在掌握 R 处理GEO芯片数据后，我们以TCGA的数据为例，分享如何用 R 处理二代测序数据。

图片

对于生信论文，大趋势是临床样本测序 生信验证 实验验证，因为高分论文都是这么做的。二代测序已经非常普及！TCGA数据库中存放了各种肿瘤临床标本的多组学数据，样本信息也非常全面，是肿瘤研究极为重要的资源！仅仅是掌握在线工具已经无法满足个性化分析需求。因此，我们应该掌握数据下载、处理和可视化的方法。

图片

TCGAbiolinks是官方唯一推荐的专用下载及分析可视化一体的 R 包。TCGAbiolinks可以下载全部的TCGA数据（除了受限制的）。除了常规的转录组数据，还包括甲基化数据、SNP数据、突变数据、临床信息等数据；还能进行数据分析，例如差异分析、生存分析、聚类分析等。此外，TCGAbiolinks还具有强大的绘图功能，可绘制突变瀑布图等。

图片

首先环境准备，获取数据，清洗数据。##########--------TCGA实操教程---------############## 加载 R 包，准备环境## install.packages('BiocManager')## BiocManager::install(c('SummarizedExperiment','TCGAbiolinks'))library(SummarizedExperiment)library(TCGAbiolinks)#### 获取数据query.exp <- GDCquery( project = 'TCGA-LUAD', data.category = 'Transcriptome Profiling', data.type = 'Gene Expression Quantification', workflow.type = 'STAR - Counts', sample.type = c('Primary Tumor','Solid Tissue Normal'))GDCdownload( query = query.exp, files.per.chunk = 100)luad.exp <- GDCprepare( query = query.exp, save = TRUE, save.filename = 'luadExp.rda')# 获取亚型信息infomation.subtype <- TCGAquery_subtype(tumor = 'LUAD')# 获取临床信息information.clinical <- GDCquery_clinic(project = 'TCGA-LUAD', type = 'clinical') # 肿瘤样本samples.primary.tumour <- luad.exp$barcode[luad.exp$shortLetterCode == 'TP']# 对照样本samples.solid.tissue.normal <- luad.exp$barcode[luad.exp$shortLetterCode == 'NT']然后，差异分析，富集分析。

##########--------差异分析---------##########dataPrep <- TCGAanalyze_Preprocessing(  object = luad.exp,   cor.cut = 0.6)                      # BiocManager::install(c('EDAseq','DESeq2'))dataNorm <- TCGAanalyze_Normalization(  tabDF = dataPrep,  geneInfo = geneInfoHT,  method = 'gcContent')                dataFilt <- TCGAanalyze_Filtering(  tabDF = dataNorm,  method = 'quantile',   qnt.cut =  0.25)   # BiocManager::install('edgeR')dataDEGs <- TCGAanalyze_DEA(  mat1 = dataFilt[,samples.solid.tissue.normal],  mat2 = dataFilt[,samples.primary.tumour],  Cond1type = 'Normal',  Cond2type = 'Tumor',  fdr.cut = 0.01 ,  logFC.cut = 2,  method = 'glmLRT',  pipeline = 'edgeR')  

数据可视化，绘制火山图##########--------TCGA数据可视化--------########### 加载包library(tidyverse)library(ggplot2)matrix <- dataDEGs[dataDEGs$gene_type =='protein_coding',]#### 用aggregate包，取最大值（强烈推荐）matrix <- aggregate( . ~ gene_name,data = matrix, max) rownames(matrix) <- matrix$gene_namecolnames(matrix) <- c('Symbol', 'logFC','logCPM','LR','Pvalue','adj.P.Val','genetype')str(matrix) ## 字符型向量，无法进行运算allDiff <- as.data.frame(lapply(matrix[,c(2:6)],as.numeric))rownames(allDiff) <- rownames(matrix)### 火山图xMax <- max(-log10(allDiff$adj.P.Val)) yMax <- max(abs(allDiff$logFC))library(ggplot2) ##可视化的关键包allDiff$change <- ifelse(allDiff$adj.P.Val < 0.05 & abs(allDiff$logFC) > 1, ifelse(allDiff$logFC > 1,'UP','DOWN'), 'NOT')table(allDiff$change)p1 <- ggplot(data = allDiff, aes(x = -log10(adj.P.Val), y = logFC, color = change)) geom_point(alpha=0.8, size = 1) theme_bw(base_size = 15) theme(plot.title = element_text(hjust=0.5), # 标题居中 panel.grid.minor = element_blank(), panel.grid.major = element_blank()) # 网格线设置为空白 geom_hline(yintercept = 0 ,linetype = 2 ) scale_color_manual(name = '', values = c('red', 'blue', 'black'), limits = c('UP', 'DOWN', 'NOT')) xlim(0,xMax) ylim(-yMax,yMax) labs(title = 'Volcano', x = '-Log10(adj.P.Val)', y = 'LogFC')p1### 美化火山图##数据整理和条件设置data0 <- allDiffdata1 <- data0 %>% rownames_to_column('Genes') #行名转为Genes为列名的一列data2 <- data1 %>% mutate(regulate = case_when(logFC >= 1 & adj.P.Val <= 0.05 ~ 'up', logFC <= -1 & adj.P.Val <= 0.05 ~ 'down', TRUE ~ 'NS')) ## 使用ggplot2包绘制火山图# 基本绘图library(ggplot2)ggplot(data2,aes(logFC,-log10(adj.P.Val))) geom_point() labs(x=expression(Log[2]*' Fold Change'), y=expression(-Log[10]*' (p value)')) #修改坐标轴命名细节# 美化绘图ggplot(data2,aes(logFC,-log10(adj.P.Val), #分别给正负显著变化的基因在图中根据颜色、大小标注出来 color=factor(regulate), size=factor(regulate))) geom_point() labs(x=expression(Log[2]*' Fold Change'), y=expression(-Log[10]*' (p value)')) theme_grey(base_size = 15) scale_color_manual(values = c('blue','grey','red')) scale_size_manual(values = c(2,1,2)) theme(legend.title = element_blank(), #图例的设置参数 legend.position = 'right', #标签位置为right legend.background = element_rect(fill='transparent'))通过上述代码，我们就可以获得清洁数据和差异基因，进行后续分析了。 本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。