TCGA差异表达分析|2022.5.1更新

news/2024/11/28 19:34:05/

作者：Squirrelity

2022-07-18 补充说明
最近R更新了，很多包都用不了，如果遇到报错或者是运行不了有可能是因为版本问题。

一、加载对应的R包
这里用到十三个包（距离上次更新之后又新增了不少方法/包）：

library("TCGAbiolinks")
library("plyr")
library("limma")
library("biomaRt")
library("SummarizedExperiment")
library("stringr")
library("ggplotify")
library("patchwork")
library("cowplot")
library('DESeq2')
library('edgeR')
library("dplyr")
library("rtracklayer")

光下载都费了不少功夫www，下面把install代码放出来(｡･∀･)ﾉﾞ直接install.packages()日常失败我就不说了
大部分包都可以在bioconductor中找到，有遗漏的可以去官网找下载代码
https://bioconductor.org/

#plyr下载 唯一可以直接install的包啊哈哈
install.packages("plyr")
#if(!require(stringr))
#install.packages（）下载法
if(!require(stringr))install.packages('stringr')
if(!require(ggplotify))install.packages("ggplotify")
if(!require(patchwork))install.packages("patchwork")
if(!require(cowplot))install.packages("cowplot")#更新BiocManager3.15，R版本为4.2.0
if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install(version = "3.15")#tccgbiolinks包稳定版本安装
if (!requireNamespace("BiocManager", quietly=TRUE))install.packages("BiocManager")
BiocManager::install("TCGAbiolinks")#limma包下载
if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install("limma")#biomaRt包下载
if (!require("BiocManager", quietly = TRUE))install.packages("biomaRt")
BiocManager::install("biomaRt")#最麻烦的SummarizedExperiment包：force = TRUE是根据warning后来加的if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install("SummarizedExperiment",force = TRUE)
#edgeR包下载
if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install("edgeR")
#DESeq2包下载
if (!requireNamespace("BiocManager", quietly=TRUE))install.packages("BiocManager")
BiocManager::install("DESeq2")

在使用之前需要加载BIocManager，代码参考这个：https://bioconductor.org/install/

if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install(version = "3.15")

若提示warning：A version of this package for your version of R might be available elsewhere,see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

可以的解决方法有：

1.乖乖在官网下载cran（最不推荐）
2.在Rstudio上方的tools-global options处找到packages，修改默认的global cran，选择apply-ok即可
3.找到包的官网，用官网提供的代码下载（首推♥♥♥）

二、下载数据
在下载之前设置工作路径：
dir.create（）创建目录，getwd（）获取工作路径，setwd（）设置工作路径，由于TCGA下载下来的数据包都挺大的，建议还是选一个比较富裕的盘来作为工作路径。

dir.create("D:\\BioInfoCloud\\TCGABiolinks\\case_study")
setwd("D:\\BioInfoCloud\\TCGABiolinks\\case_study")

这里用到的是R包TCGAbiolinks：
可以参照R包TCGAbiolinks官网使用http://www.bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/casestudy.html#Case_study_n_1:_Pan_Cancer_downstream_analysis_BRCA
示例：

query <- GDCquery(project = "TCGA-BRCA",data.category = "Transcriptome Profiling",data.type = "Gene Expression Quantification", workflow.type = "STAR - Counts")

project选出的是肿瘤项目，而里面用到的都是缩写，详见https://blog.csdn.net/Squirrelity/article/details/124259330?spm=1001.2014.3001.5501

建议去TCGA官网repository一边对照着选所需要的样本
https://portal.gdc.cancer.gov/repository?facetTab=cases

#下载到本地
GDCdownload(query = query, method = "api")
#查看下载的数据
View(query)
BRCA.Rnaseq.SE <- GDCprepare(    query = query, save = TRUE, save.filename = "brca.rda")
BRCAMatrix <- assay(BRCA.Rnaseq.SE,"unstranded")
#记住这个文件BRCA.RNAseq_CorOutliers
BRCA.RNAseq_CorOutliers <- TCGAanalyze_Preprocessing(BRCA.Rnaseq.SE)

三、ID转换
下载下来的BRCA.RNAseq_CorOutliers为ENTREZID，而我们肯定是需要图片显示基因名而不是ENTREZID，因此进行数据转换，这里用到包括但不限于的dplyr包和rtracklayer包

library("dplyr")
library("rtracklayer")

ID转换分为四步：
1.导入数据：BRCA.RNAseq_CorOutliers和人类基因组注释文件；

data=read.table(BRCA.RNAseq_CorOutliers,header=T,sep='\t')
#把行名改为gene_id，与gtf保持一致
colnames(data)[1] <- "gene_id"

对照人类基因组注释文件，对BRCA.RNAseq_CorOutliers进行ID转换
其中，人类基因组注释文件参考http://www.360doc.com/content/21/1028/10/77506210_1001626502.shtml

#处理人类基因组注释文件的数据
gtf <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf.gz')
gtf <- as.data.frame(gtf)
save(gtf,file="人类基因组注释文件.Rda")
gtf <-load(file="人类基因组注释文件.Rda")
#根据条件筛选基因（大筛选）
a = dplyr::filter(gtf,type=="gene")
dim（a）
#只要gene_name,gene_id,gene_biotype这三行
b = dplyr::select(a,c(gene_name,gene_id,gene_biotype))

2.数据预处理

#ENTREZID带有，这里去除小数点及后边的数字（我用excel处理的，ctrl+F无字符替换.*）
data1 <- separate(data,gene_id,into = c("gene_id"),sep="\\.")

3.数据处理

#根据gene id 合并文件
c = dplyr::inner_join(b,data,by="gene_id")
#去掉2,3列，基因名再去重
d=select(c,-gene_id,-gene_biotype)
data1=distinct(d,gene_name,.keep_all = T)
#把行名由数字换成基因
rownames(data1)<- data1[,1]
data1<-data1[,-1]

4.数据后处理

#如下载的数据取了log2（count-1）这里再返回count
data2 <- 2^data1 -1
write.csv(data2,file="data2.csv")
data2 <- read.csv("data2.csv")
#重新用read打开整行的-会变成.因此需要恢复原来的行名
colnames(data2) <- colnames(BRCA.RNAseq_CorOutliers)

四、差异表达分析
1.参考网址
代码示例参考这个包的文档http://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html
TCGA可视化教程
https://www.jianshu.com/p/d3e481f0187a
https://cloud.tencent.com/developer/article/1778874
http://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/casestudy.html

2.数据预处理

#处理后可得热图判断样本相似性
dataPrep <- TCGAanalyze_Preprocessing(object = BRCA.RNAseq_CorOutliers, cor.cut = 0.6
)

*上面代码生成的图如下所示。组内的样本相似性都很高，符合预期。在这里插入图片描述

3.对数据进行标准化处理+质控+差异化分析
**TCGAanalyze_LevelTab（）**将差异表达基因在正常和肿瘤组织中的表达量数据添加到差异表达分析结果中的主要用法：

TCGAanalyze_LevelTab(FC_FDR_table_mRNA, typeCond1, typeCond2, TableCond1,
TableCond2, typeOrder = TRUE)

在这里插入图片描述

            
#数据标准化
dataNorm <- TCGAanalyze_Normalization(tabDF = data2,geneInfo = geneInfoHT,method = "gcContent"
)                
#数据质控
dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,method = "quantile", qnt.cut =  0.25
)   #分组：NT正常组织组 TP癌症组织组*# selection of normal samples "NT"
samplesNT <- TCGAquery_SampleTypes(barcode = colnames(dataFilt),typesample = c("NT")
))
# selection of tumor samples "TP"
samplesTP <- TCGAquery_SampleTypes(barcode = colnames(dataFilt), typesample = c("TP")
)# 差异表达分析
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"
)  #在正常和肿瘤样本中差异基因的表达值
dataDEGsFiltLevel <- TCGAanalyze_LevelTab(FC_FDR_table_mRNA = dataDEGs,typeCond1 = "Tumor",typeCond2 = "Normal",TableCond1 = dataFilt[,samplesTP1],TableCond2 = dataFilt[,samplesNT1]
)

得到的dataDEGsFiltLevel文件按logFC绝对值排序可得最显著的top差异表达基因（excel表处理）
在这里插入图片描述

五、可视化
ps：得到的图片有的可以直接看，有的保存在工作路径上了。
1.PCA主成分分析

TCGAvisualize_PCA()实现主成分分析的主要用法:

TCGAvisualize_PCA(dataFilt, dataDEGsFiltLevel, ntopgenes, group1, group2)

在这里插入图片描述

#标准化
dataNorm <- TCGAbiolinks::TCGAanalyze_Normalization(data2, geneInfo)#质量控制
dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,method = "quantile", qnt.cut =  0.25)#选择正常样本
group1 <- TCGAquery_SampleTypes(colnames(dataFilt), typesample = c("NT"))
#选择癌症样本
group2 <- TCGAquery_SampleTypes(colnames(dataFilt), typesample = c("TP"))#Principal Component Analysis plot for ntop selected DEGs
pca.top200 <- TCGAvisualize_PCA(dataFilt,dataDEGsFiltLevel, ntopgenes = 200,group1, group2)

上面代码生成的图如下所示。
在这里插入图片描述

2.火山图

#为了做图的需要，突出显示logFC≥8的gene名称
DEG.BRCA.filt<-dataDEGs[which(abs(dataDEGs$logFC) >= 8), ]
str(DEG.BRCA.filt)
#'data.frame':    29 obs. of  5 variables:
#说明共有29个基因满足|logFC|≥8TCGAVisualize_volcano(dataDEGs$logFC, dataDEGs$FDR,filename = "TumorvsNormal_FC8.edgeR.pdf", xlab = "logFC",names = rownames(dataDEGs), show.names = "highlighted",x.cut = 1, y.cut = 0.01, highlight = rownames(dataDEGs)[which(abs(DEG.LIHC.edgeR$logFC) >= 8)],highlight.color = "orange",title = "volcano plot by edgeR")

上面代码生成的图如下所示。突出显示了logFC≥8的gene名称
在这里插入图片描述

3.GO功能分析条形图
TCGAbiolinks 输出条形图，其中包含三个本体的主要类别（分别为GO：生物过程、GO：细胞成分和GO：分子功能）的基因数量。

ansEA <- TCGAanalyze_EAcomplete(TFname = "DEA genes Normal Vs Tumor",RegulonList = dataDEGs$gene_name
)  TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP),GOBPTab = ansEA$ResBP,GOCCTab = ansEA$ResCC,GOMFTab = ansEA$ResMF,PathTab = ansEA$ResPat,nRGTab = dataDEGs$gene_name,nBar = 10
)

上面代码生成的图如下所示。
上面代码生成的图如下所示。


http://www.ppmy.cn/news/966965.html

相关文章

TCGA/GTEx泛癌数据任意基因表达量展示

TCGA/GTEx泛癌数据任意基因表达量展示

有了泛癌的数据之后就可以进行各种分析了&#xff0c;当然这些都是在R语言的基础上进行的。如果你不会R语言&#xff0c;也可以通过各种各样的网页工具实现。 我们今天就简单展示下任意基因在泛癌图谱中的表达量情况。 TCGA,GTEx,TCGAGTEx的泛癌数据都整理好了&#xff0c;大…
阅读更多...
GATK变异检测

GATK变异检测

变异检测 GATK 变异检测 GATK是Genome Analysis Toolkit的缩写&#xff0c;是用来处理高通量测序数据的一套软件。最初&#xff0c;GATK被设计用来分析人类基因组和外显子&#xff0c;主要用来寻找SNP和indel。后开&#xff0c;GATK的功能越来越丰富&#xff0c;增加了short v…
阅读更多...
四、肿瘤全基因组学体细胞点突变特征（The repertoire of mutational signatures in human cancer）

四、肿瘤全基因组学体细胞点突变特征（The repertoire of mutational signatures in human cancer）

全文链接 一、肿瘤突变特征&#xff1a;碱基置换及插入、缺失突变 单碱基置换&#xff08;49种特征类型&#xff0c;single-base-substitution&#xff0c;SBS&#xff09; 双碱基置换&#xff08;11种特征类型&#xff0c;doublet-base-substitution&#xff0c;DBS&#xf…
阅读更多...
TCGA_联合GTEx分析2_查看批次效应

TCGA_联合GTEx分析2_查看批次效应

在 TCGA_联合GTEx分析1_得到表达矩阵.tpm_老实人谢耳朵的博客-CSDN博客 中&#xff0c;获取了TCGA和GTEx中样本的表达矩阵数据&#xff0c;数据格式均为tpm。本文对二者进行合并后&#xff0c;通过PCA分析、绘制内参箱线图等方法&#xff0c;查看是否存在批次效应。 关于批次效…
阅读更多...
TCGA 亚型突变负荷代码

TCGA 亚型突变负荷代码

#1、准备文件/数据并加载相应的包 #1.1下载并加载相应的包&#xff0c;有就直接加载&#xff0c;没有就下载后再加载。 install.packages("pacman") library(pacman) p_load(TCGAbiolinks,DT,tidyverse) BiocManager::install("TCGAbiolinks") library(t…
阅读更多...
基因组变异检测SNPcalling（GATK）

基因组变异检测SNPcalling（GATK）

基因组变异检测SNPcalling&#xff08;GATK&#xff09; 第一步&#xff0c;将数据软连接到自己的工作路径下第二步&#xff0c;用BWA index功能为ref文件做index&#xff0c;用Samtools软件为ref做index第三步&#xff0c;用BWA软件做序列比对&#xff0c;得到序列比对的sam文…
阅读更多...
提取TCGA 中体细胞突变数据的表达矩阵

提取TCGA 中体细胞突变数据的表达矩阵

#因为之前的命令调用GDCquery_Maf 发现用不了 #故找到了一些其他的方法&#xff0c;并且自己试着将其弄成了一个表达矩阵。 #代码如下 #1、下载加载相应的包 install.packages("pacman") library(pacman) p_load(TCGAbiolinks,DT,tidyverse) BiocManager::insta…
阅读更多...
maftools|TCGA肿瘤突变数据的汇总，分析和可视化

maftools|TCGA肿瘤突变数据的汇总，分析和可视化

之前介绍了使用maftools | 从头开始绘制发表级oncoplot&#xff08;瀑布图&#xff09; R-maftools包绘制组学突变结果&#xff08;MAF&#xff09;的oncoplot或者叫“瀑布图”&#xff0c;以及一些细节的更改和注释。 本文继续介绍maftools对于MAF文件的其他应用&#xff0c;为…
阅读更多...
最新文章

Copyright @ 2022~2023