因为是癌症方面，自己不研究这一方面，所以不常用，但是GEO的转录组数据，是根据这个文件改写的

0.安装包

options("repos" = c(CRAN="http://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")cran_packages <- c('tidyr','tibble','dplyr','stringr','ggplot2','ggpubr','factoextra','FactoMineR','pheatmap',"survival","survminer","patchwork","ggstatsplot","ggplotify","cowplot","glmnet","ROCR","caret","randomForest","survminer","Hmisc","e1071","deconstructSigs","AnnoProbe","timeROC","circlize","VennDiagram","tinyarray"
)
Biocductor_packages <- c("limma","clusterProfiler","org.Hs.eg.db","SummarizedExperiment","DESeq2","edgeR","ggpubr","rtracklayer","genefilter","maftools","ComplexHeatmap","GDCRNATools"
)for (pkg in cran_packages){if (! require(pkg,character.only=T) ) {install.packages(pkg,ask = F,update = F)require(pkg,character.only=T) }
}# use BiocManager to install
for (pkg in Biocductor_packages){if (! require(pkg,character.only=T) ) {BiocManager::install(pkg,ask = F,update = F)require(pkg,character.only=T) }
}#前面的任何信息都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){require(pkg,character.only=T)
}#没有任何提示就是成功了，如果有warningxx包不存在，用library检查一下。#报错就回去重新安装。如果你没有安装xx包，却提示你xx包不存在，这也正常，是因为依赖关系，缺啥补啥。

1. 数据下载

1.数据下载

最方便的是xena。可以网页下载，也可以用代码下载。

proj = "TCGA-CHOL"
if(F){download.file(url = paste0("https://gdc.xenahubs.net/download/",proj, ".htseq_counts.tsv.gz"),destfile = paste0(proj,".htseq_counts.tsv.gz"))download.file(url = paste0("https://gdc.xenahubs.net/download/",proj, ".GDC_phenotype.tsv.gz"),destfile = paste0(proj,".GDC_phenotype.tsv.gz"))download.file(url = paste0("https://gdc.xenahubs.net/download/",proj, ".survival.tsv"),destfile = paste0(proj,".survival.tsv"))
}

2.生存信息与临床信息

这里仅仅是查看一下，到生存信息部分再整理。

clinical = read.delim(paste0(proj,".GDC_phenotype.tsv.gz"),fill = T,header = T,sep = "\t")
surv = read.delim(paste0(proj,".survival.tsv"),header = T)
clinical[1:4,1:4]
head(surv)

3.表达矩阵行名ID转换

dat = read.table(paste0(proj,".htseq_counts.tsv.gz"),check.names = F,row.names = 1,header = T)
#逆转log
dat = as.matrix(2^dat - 1)
dat[1:4,1:4]
# 深坑一个
dat[97,9]
as.character(dat[97,9]) #眼见不一定为实吧。# 用apply转换为整数矩阵
exp = apply(dat, 2, as.integer)
exp[1:4,1:4] #行名消失
rownames(exp) = rownames(dat)
exp[1:4,1:4]

gdc下载的数据从此处开始衔接

library(stringr)
# 行名ID转换：方法1(推荐)
head(rownames(exp))
library(AnnoProbe)
annoGene(rownames(exp),ID_type = "ENSEMBL") # 转换不了
rownames(exp) = str_split(rownames(exp),"\\.",simplify = T)[,1];head(rownames(exp))
re = annoGene(rownames(exp),ID_type = "ENSEMBL");head(re)
library(tinyarray)
exp = trans_array(exp,ids = re,from = "ENSEMBL",to = "SYMBOL")
exp[1:4,1:4]
# 方法2： zz_gtf_ID转换.Rmd

4.基因过滤

需要过滤一下那些在很多样本里表达量都为0或者表达量很低的基因。过滤标准不唯一。

过滤之前基因数量：

nrow(exp)

常用过滤标准1：

仅去除在所有样本里表达量都为零的基因

exp1 = exp[rowSums(exp)>0,]
nrow(exp1)

常用过滤标准2(推荐)：

仅保留在一半以上样本里表达的基因

exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]
nrow(exp)

5.分组信息获取

根据样本ID的第14-15位，给样本分组（tumor和normal）

table(str_sub(colnames(exp),14,15))
Group = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
Group = factor(Group,levels = c("normal","tumor"))
table(Group)

6.保存数据

TCGA以外的数据没有clinical,surv，从下面代码里去掉。

save(exp,Group,proj,clinical,surv,file = paste0(proj,".Rdata"))

其他数据下载方法

gdc-client

是从官方网站下载，代码在gdc文件夹，难度较大，结果与xena整理的一致，作为补充学习浏览一下即可。

GDCRNAtools

http://bioconductor.org/packages/devel/bioc/vignettes/GDCRNATools/inst/doc/GDCRNATools.html

1.三大R包差异分析

rm(list = ls())
load("TCGA-CHOL.Rdata")
table(Group)
#deseq2----
library(DESeq2)
colData <- data.frame(row.names =colnames(exp), condition=Group)
if(!file.exists(paste0(proj,"_dd.Rdata"))){dds <- DESeqDataSetFromMatrix(countData = exp,colData = colData,design = ~ condition)dds <- DESeq(dds)save(dds,file = paste0(proj,"_dd.Rdata"))
}
load(file = paste0(proj,"_dd.Rdata"))
class(dds)
res <- results(dds, contrast = c("condition",rev(levels(Group))))
#constrast
c("condition",rev(levels(Group)))
class(res)
DEG1 <- as.data.frame(res)
DEG1 <- DEG1[order(DEG1$pvalue),]
DEG1 = na.omit(DEG1)
head(DEG1)#添加change列标记基因上调下调
logFC_t = 2
pvalue_t = 0.05k1 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange < -logFC_t);table(k1)
k2 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange > logFC_t);table(k2)
DEG1$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG1$change)
head(DEG1)#edgeR----
library(edgeR)dge <- DGEList(counts=exp,group=Group)
dge$samples$lib.size <- colSums(dge$counts)
dge <- calcNormFactors(dge) design <- model.matrix(~Group)dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)fit <- glmFit(dge, design)
fit <- glmLRT(fit) DEG2=topTags(fit, n=Inf)
class(DEG2)
DEG2=as.data.frame(DEG2)
head(DEG2)k1 = (DEG2$PValue < pvalue_t)&(DEG2$logFC < -logFC_t);table(k1)
k2 = (DEG2$PValue < pvalue_t)&(DEG2$logFC > logFC_t);table(k2)
DEG2$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))head(DEG2)
table(DEG2$change)
###limma----
library(limma)
dge <- edgeR::DGEList(counts=exp)
dge <- edgeR::calcNormFactors(dge)
design <- model.matrix(~Group)
v <- voom(dge,design, normalize="quantile")design <- model.matrix(~Group)
fit <- lmFit(v, design)
fit= eBayes(fit)DEG3 = topTable(fit, coef=2, n=Inf)
DEG3 = na.omit(DEG3)k1 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC < -logFC_t);table(k1)
k2 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC > logFC_t);table(k2)
DEG3$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG3$change)
head(DEG3)tj = data.frame(deseq2 = as.integer(table(DEG1$change)),edgeR = as.integer(table(DEG2$change)),limma_voom = as.integer(table(DEG3$change)),row.names = c("down","not","up"));tj
save(DEG1,DEG2,DEG3,Group,tj,file = paste0(proj,"_DEG.Rdata"))

2.可视化

library(ggplot2)
library(tinyarray)
exp[1:4,1:4]
# cpm 去除文库大小的影响
dat = log2(cpm(exp)+1)
pca.plot = draw_pca(dat,Group);pca.plot
save(pca.plot,file = paste0(proj,"_pcaplot.Rdata"))cg1 = rownames(DEG1)[DEG1$change !="NOT"]
cg2 = rownames(DEG2)[DEG2$change !="NOT"]
cg3 = rownames(DEG3)[DEG3$change !="NOT"]h1 = draw_heatmap(dat[cg1,],Group,n_cutoff = 2)
h2 = draw_heatmap(dat[cg2,],Group,n_cutoff = 2)
h3 = draw_heatmap(dat[cg3,],Group,n_cutoff = 2)v1 = draw_volcano(DEG1,pkg = 1,logFC_cutoff = logFC_t)
v2 = draw_volcano(DEG2,pkg = 2,logFC_cutoff = logFC_t)
v3 = draw_volcano(DEG3,pkg = 3,logFC_cutoff = logFC_t)library(patchwork)
(h1 + h2 + h3) / (v1 + v2 + v3) +plot_layout(guides = 'collect') &theme(legend.position = "none")ggsave(paste0(proj,"_heat_vo.png"),width = 15,height = 10)

3.三大R包差异基因对比

UP=function(df){rownames(df)[df$change=="UP"]
}
DOWN=function(df){rownames(df)[df$change=="DOWN"]
}up = intersect(intersect(UP(DEG1),UP(DEG2)),UP(DEG3))
down = intersect(intersect(DOWN(DEG1),DOWN(DEG2)),DOWN(DEG3))
dat = log2(cpm(exp)+1)
hp = draw_heatmap(dat[c(up,down),],Group,n_cutoff = 2)#上调、下调基因分别画维恩图
up_genes = list(Deseq2 = UP(DEG1),edgeR = UP(DEG2),limma = UP(DEG3))down_genes = list(Deseq2 = DOWN(DEG1),edgeR = DOWN(DEG2),limma = DOWN(DEG3))up.plot <- draw_venn(up_genes,"UPgene")
down.plot <- draw_venn(down_genes,"DOWNgene")#维恩图拼图，终于搞定library(patchwork)
#up.plot + down.plot
# 拼图
pca.plot + hp+up.plot +down.plot+ plot_layout(guides = "collect")
ggsave(paste0(proj,"_heat_ve_pca.png"),width = 15,height = 10)

6. 生信技能树——TCGA癌症数据1相关推荐

CancerSubtypes包的介绍（根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理）
CancerSubtypes包的介绍(根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理,感谢Jimmy老师!) 1. 引言 2. 数据处理 2.1 基本处理 2.1.1 通过检查数据分布来分 ...
ProTICS包的介绍（根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理）
ProTICS包的介绍(根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理,感谢Jimmy老师!) 1.设置环境 2.Part1的结果 3.Part2的结果 4.Part3的结果 5.相关函数 ...
R：生信技能树学习笔记一
生信技能树小破站:R应该这样学1-4 1.查看已经安装的包的地址 .libPaths() 2.怎么查看函数用法 #在RStudio的右下角窗口的help可以看到 ?函数名 3.三个有用的函数 1.he ...
R：生信技能树学习笔记二
生信技能树小破站:R应该这样学5-7 1.热图 rm(list=ls()) library(pheatmap) a1=rnorm(100) dim(a1)=c(5,20) #设置维度 pheatmap ...
这是入门生信，学习生信分析思路和数据可视化的首选？
封面来源:https://www.zhihu.com/question/304747766 常规转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,是大部分CNS必备的技术,以后 ...
.md是什么文件_生信中常见的数据文件格式
TCGA | GEO | 文献阅读 | 数据库 | 理论知识 R语言 | Bioconductor | 服务器与Linux 前面我们介绍了各种测序技术的原理:illumina.Sanger.第三代和第 ...
生信技能树【代码大全搜录】
生信技能术代码大全: rm(list = ls()) options()$repos options()$BioC_mirror #options(BioC_mirror="https:// ...
TCGA学习笔记一(生信技能树概述版)
1.背景介绍重要数据外显子数据表达数据小RNA测序数据拷贝数芯片甲基化数据蛋白质组学数据临床信息癌症背景知识网页工具大全 GDC cbioportal:按照paper来分类的 UC ...
生信工具 | TCGA数据分析工具GEPIA最新更新，用于免疫细胞浸润分析
GEPIA(http://gepia.cancer-pku.cn/index.html)这个工具可以说是分析TCGA数据库数据分析工具中比较简单好用的工具了,包括生存分析,表达差异分析,相关性分析等, ...
生信技能树课程记录笔记（七）20220531
一.数据框排序法一:sort函数默认升序. 例:sort(test$Sepal.Length) 法二:order函数默认升序,返回数据下标组成的数组. 可以给向量排序,也可以给数据框排序例:t ...

6. 生信技能树——TCGA癌症数据1