热点终结者？泛癌-目标基因集在肿瘤组织与正常组织的差异分析

热点终结者？泛癌-目标基因集在肿瘤组织与正常组织的差异分析

近年来，多维组学联合分析成为大趋势，经常出现在高分文章中，泛癌分析在现在还是非常火的，而“泛癌”pancancer研究尤其火热，在pubmed中搜索一个”pancer”,就可以发现今年来相关文献发表量在逐年增加，为了跟上时代潮流，小果今天为大家分享目标基因集在泛癌与正常组织的差异分析奥，那就和小果开始今天的学习之旅吧！

1. 如何做泛癌与正常组织差异分析？

在进行泛癌差异分析之前，小果想带着大家进行了解一下何为泛癌？泛癌就是在研究当下的癌种后，在扩大到其他癌种，可以是所有肿瘤的大泛癌，也可以是呼吸系统，消化系统的小泛癌奥！小果的解释还算简单易懂吧！

接下来小果为小伙伴讲讲如何做泛癌差异分析，主要包括泛癌表达矩阵和样本信息文件下载，数据处理获得目标基因集在不同肿瘤组织中表达矩阵，最后利用limma包进行差异分析，获得差异分析结果文件，这就是基本的分析流程啦。

1. 需要的R包

本次分析需要的R包最主要的是limma包进行差异分析，温馨提示，需要注意的是在加载R包时，千万不要忘记加载twoclasslimma.R脚本，否则会遇到报错奥！，可以通过以下命令进行安装，小果已经尝试过的奥！

#安装需要的R包

install.packages(“ggplot2”)

install.packages(“data.table”)

BiocManager::install(“impute”)

BiocManager::install(“limma”)

install.packages(“tidyverse”)

install.packages(“ggpubr”)

#导入需要的R包

library(ggplot2)

library(tidyverse)

library(impute)

library(limma)

library(ggpubr)

source(“twoclasslimma.R”)

1. 数据下载

在做泛癌分析之前，最重要的问题是泛癌基因表达矩阵和样本信息下载，小果所用的泛癌数据是来自该网站https://gdc.cancer.gov/about-data/publications/pancanatlas，只需要登录该网站就可以很方便的下载到自己想要的文件奥，由于泛癌基因表达矩阵文件较大，下载会有点慢，有需要的可以call 小果哈。。。。

#泛癌样本信息文件,主要包括样本名和肿瘤组织类型，其他列可以忽略奥！

merged_sample_quality_annotations.tsv

#泛癌基因表达矩阵文件，行名为基因名，列名为样本名，认真的小伙伴会发现基因名不是我们需要的奥，下面的分析中会进行转化操作奥，可以重点关注奥！

EBPlusPlusAdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.tsv

1. 目标基因集在不同肿瘤组织中差异分析

在成功的下载到数据后，就要进入到最主要的环节了，进行目标基因集在不同肿瘤组织中差异分析的工作了，首先准备感兴趣的肿瘤名称和感兴趣的基因名称，然后提取感兴趣基因集在不同肿瘤组织中的表达矩阵文件，然后利用limma软件进行不同肿瘤组织与正常组织差异分析，获得差异分析结果文件，该分析中包含大量的数据处理技巧，小果强烈推介认真在认真阅读奥。

# 感兴趣的肿瘤名称

tumors <- c(“BLCA”,”BRCA”,”CESC”,”CHOL”,”COAD”,

“ESCA”,”GBM”,”HNSC”,”KICH”,”KIRC”,

“KIRP”,”LIHC”,”LUAD”,”LUSC”,”PAAD”,

“PRAD”,”READ”,”STAD”,”THCA”,”UCEC”)

#感兴趣的基因名称

Gene<- c(“CDKN1A”,”HSPA5″,”TTC35″,”SLC7A11″,”NFE2L2″,”MT1G”,”HSPB1″,”GPX4″,”FANCD2″,”CISD1″,”FDFT1″,”SLC1A5″,”SAT1″,”TFRC”,”RPL8″,”NCOA4″,”LPCAT3″,”GLS2″,”DPP4″,”CS”,”CARS”,”ATP5G3″,”ALOX15″,”ACSL4″,”EMC2″)

#处理和简化样本信息文件(此处开始收费)

anno<-read.delim(“merged_sample_quality_annotations.tsv”,sep=”\t”,row.names=NULL,check.names = F,stringsAsFactors = F,header = T)

#提取aliquot_barcode的1-15个字符串

anno$simple_barcode <- substr(anno$aliquot_barcode,1,15)

#提取cancer_type和simple_barcode两列

samAnno <- anno[!duplicated(anno$simple_barcode),c(“cancer type”, “simple_barcode”)]

#保留有cancer type不为空的样本

samAnno <- samAnno[which(samAnno$`cancer type` != “”),]

#保存简化后的样本信息文件

write.table(samAnno,”simple_sample_annotation.txt”,sep = “\t”,row.names = F,col.names = T,quote = F)

#表达矩阵文件处理

expr<-fread(“EBPlusPlusAdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.tsv”,sep= “\t”,stringsAsFactors = F,check.names = F,header = T)

#将数据转化为数据框格式

expr <- as.data.frame(expr)

#第一列作为行名

rownames(expr) <- expr[,1]

#删除第一列

expr <- expr[,-1]

#对行名通过“|”来分割并保存第一个分割变量

gene <- sapply(strsplit(rownames(expr),”|”,fixed = T), “[“,1)

expr$gene <- gene

# 移除重复样本

expr <- expr[!duplicated(expr$gene),] # 移除重复样本

rownames(expr) <- expr$gene;

#删除最后一列

expr <- expr[,-ncol(expr)]

#提取感兴趣基因集的表达矩阵文件

comgene <- intersect(rownames(expr),Gene)

expr_sub <- expr[comgene,]

colnames(expr_sub) <- substr(colnames(expr_sub),1,15)

expr_sub <- expr_sub[,!duplicated(colnames(expr_sub))]

#通过循环的方式进行每种肿瘤的差异分析

exprTab <- ndegs <- NULL

# 设置差异表达的阈值

log2fc.cutoff <- log2(1.5)

fdr.cutoff <- 0.05

for (i in tumors) {

message(“–“,i,”…”)

#获取相应肿瘤的样本信息文件

sam <- samAnno[which(samAnno$`cancer type` == i),”simple_barcode”]

comsam <- intersect(colnames(expr_sub), sam)

# 获得肿瘤样本

tumsam <- comsam[substr(comsam,14,14) == “0”]

#获得正常样本

norsam <- comsam[substr(comsam,14,14) == “1”]

#获得对应的肿瘤和正常组织表达矩阵文件

expr_subset <- expr_sub[,c(tumsam,norsam)]

expr_subset[expr_subset < 0] <- 0

expr_subset <- as.data.frame(impute.knn(as.matrix(expr_subset))$data)

write.table(expr_subset, paste0(“TCGA_”,i,”_expr_subset.txt”),sep = “\t”,row.names = T,col.names = NA,quote = F)

subt <- data.frame(condition = rep(c(“tumor”,”normal”),c(length(tumsam),length(norsam))),

row.names = colnames(expr_subset),

stringsAsFactors = F)

twoclasslimma(subtype = subt, # 亚型分组信息

featmat = expr_subset, # 表达矩阵文件

treatVar = “tumor”, # “治疗组”的名（比较的组）

ctrlVar = “normal”, # “对照组”的名（被比较的组）

prefix = paste0(“TCGA_”,i), # 差异表达的文件的前缀

overwt = T, # 是否覆盖已经存在的差异表达文件

sort.p = F, # 是否排序p值

verbose = TRUE, # 是否简化输出

res.path = “.”) # 输出结果路径

# 加载差异表达文件

res<-read.table(paste0(“TCGA_”,i,”_limma_test_result.tumor_vs_normal.txt”),sep = “\t”,row.names = 1,check.names = F,stringsAsFactors = F,header = T)

#判断上调和下调基因

upgene <- res[which(res$log2fc > log2fc.cutoff & res$padj < fdr.cutoff),]

dngene <- res[which(res$log2fc < -log2fc.cutoff & res$padj < fdr.cutoff),]

# 基因差异表达的数目

if(nrow(upgene) > 0) {

nup <- nrow(upgene)

} else {nup <- 0}

if(nrow(dngene) > 0) {

ndn <- nrow(dngene)

} else {ndn <- 0}

#合并每种肿瘤的差异分析结果，并转化为数据框

exprTab <- rbind.data.frame(exprTab,

data.frame(gene = rownames(res),

log2fc = res$log2fc,

FDR = res$padj,

tumor = i,

stringsAsFactors = F),

stringsAsFactors = F)

#合并每种肿瘤的差异基因上下调数目结果，并转化为数据框

ndegs <- rbind.data.frame(ndegs,

data.frame(tumor = i,

Group = c(“UP”,”DOWN”),

Number = c(nup,ndn),

stringsAsFactors = F),

stringsAsFactors = F)

}

1. 绘制差异基因堆叠柱状图和基因表达量气泡图

#绘制差异基因数目堆叠柱状图

top <- ggplot(data = ndegs) +

geom_bar(mapping = aes(x = tumor, y = Number, fill = Group),

stat = ‘identity’,position = ‘stack’) +

scale_fill_manual(values = c(orange,green)) +

theme_classic() +

theme(axis.text.x = element_blank(),

axis.title.x = element_blank(),

plot.margin = unit(c(1,0,0,1), “lines”))

# 绘制泡泡图

#确定绘图基因的因子顺序

exprTab$gene<-factor(exprTab$gene, levels=rev(c(“CDKN1A”,”HSPA5″,”TTC35″,”SLC7A11″,”NFE2L2″,”MT1G”,”HSPB1″,”GPX4″,”FANCD2″,”CISD1″, “FDFT1″,”SLC1A5″,”SAT1″,”TFRC”,”RPL8″,”NCOA4″,”LPCAT3″,”GLS2″,”DPP4″,”CS”,”CARS”,”ATP5G3″,”ALOX15″,”ACSL4″)))

#设置渐变色

my_palette <- colorRampPalette(c(green,”white”,orange), alpha=TRUE)(n=128)

center <- ggplot(exprTab, aes(x=tumor,y=gene)) +

geom_point(aes(size=-log10(FDR),color=log2fc)) +

scale_color_gradientn(‘log2(FC)’,

colors=my_palette) +

theme_bw() +

#修改主题

theme(axis.text.x = element_text(angle = 45, size = 12, hjust = 0.3, vjust = 0.5, color = “black”),

axis.text.y = element_text(size = 12, color = rep(c(red,blue),c(14,10))),

axis.title = element_blank(),

panel.border = element_rect(size = 0.7, linetype = “solid”, colour = “black”),

plot.margin = unit(c(0,0,1,1), “lines”))

#利用ggpubr中的ggrange函数进行拼图

ggarrange(top,

center,

nrow = 2, ncol = 1,#排列方式两行三列

align = “v”, #排列方向

heights = c(2,6),#两张图片的大小比例

common.legend = F)#对图例不进行合并

#保存图片

ggsave(“differention_expression_pancancer.pdf”,heigth=8,weight=8)

1. 结果文件解读
2. differential_expression_of_interested_genes_in_pancancer.pdf

该图为不同肿瘤组织中差异基因堆叠柱状图和不同肿瘤组织中目标基因表达量气泡图组合图

1. TCGA_BLCA_limma_test_result.tumor_vs_normal.txt

该文件为目标基因集在相应的肿瘤组织中的limma差异分析结果

1. TCGA_BLCA_expr_subset.txt

该文件为目标基因集在肿瘤组织中的表达矩阵文件

1. simple_sample_annotation.txt

该文件为简化的样本信息文件，第一列为肿瘤类型，第二列为对应的样本名

最终小果从数据下载，差异分析和绘图，完成了铁死亡相关基因在泛癌和正常组织的差异分析；泛癌相关的分析都可以用本公司新开发的云平台生物信息分析小工具实现奥，零代码完成分析，欢迎小伙伴来尝试奥，云平台网址：http://www.biocloudservice.com/home.html，例如基因在肿瘤与正常组织中的甲基化分析(http://www.biocloudservice.com/771/771.php)，泛癌中特定通路基因的突变分析(http://www.biocloudservice.com/740/740.php)，泛癌中拷贝数变异分析(http://www.biocloudservice.com/700/700.php)等小工具；欢迎大家和小果一起讨论学习奥，下期再见奥。