2019-07-13 - TCGA胃癌的miRNA成熟体在肿瘤不
思路
首先要依据临床特征建立分组，无论是复发vs原发还是TNM分型中的T/N/M。将肿瘤每个样本的表型特征，与表达数据矩阵一一对应到同一张表格里来，随后再用作图工具出图。因为比较懒，就直接懒人方便包，从Xena数据库把做好的TCGA胃癌miRNA成熟体以及表型(phenotype)以及生存相关(survival data)的文件全部下载来，同时合并分析处理。
流程如下
1. Metadata临床数据的处理
phenotype <- data.table::fread("TCGA-BLCA.GDC_phenotype.tsv.gz")
patho_stage <- phenotype$tumor_stage.diagnoses
barcode_id <- phenotype$submitter_id.samples
抽取了胃癌表型特征，并把肿瘤级别(stage)和对应的TCGA_ID数据抽提出来。看一眼。
最初始的metadata
metadata <- metadata[substring(metadata$barcode_id,14,15)!="11",]
metadata$barcode_id <- substring(metadata$barcode_id,1,15)
colnames(metadata)[1] <- "TCGA_ID"
因为我们只关心肿瘤样本，但数据还包含有正常对照(11A)的数据，所以我们用过滤工具将正常样本清理掉，并把TCGA的ID改名。
2. 表达矩阵的获得和合并处理
将表达矩阵读取，并试图利用TCGA_ID进行合并，因为原始矩阵中TCGA_ID是列名(colnames)，所以要将原始矩阵的第1列读进列名，然后删除第1列，并进行转置。在合并过程中发现TCGA_ID中的"-"被读成了点“."，所以用gsub替换命令换了回来。
miR<- read.table("miRNA_HiSeq_gene.gz",stringsAsFactors = F,header = T)
rownames(miR) <- miR[,1]
miR <- miR[,-1]
miR <- data.frame(t(miR))
miR$TCGA_ID <- gsub("\\.","\\-",rownames(miR))
colnames(metadata)[1] <- "TCGA_ID"
miRNA数据的最初矩阵
library(dplyr)
metadata <- metadata %>%
inner_join(miR,by="TCGA_ID") %>%
select(-TCGA_ID)
再利用dplyr包的inner_join命令合并，随后删除列名中的TCGA_ID。
得到一个与stage对应的表达矩阵
3. Stage数据的预处理
因为Stage信息中包含有stage有细分成stage iia, stage iib之类的，为了简化，我们将metadata中的stage i / ii /iii都进行了合并同类项。
metadata[metadata$patho_stage %in% c("stage iiia", "stage iiib","stage iiic"),]$patho_stage <- "stage iii"
metadata[metadata$patho_stage %in% c("stage iia", "stage iib"),]$patho_stage <- "stage ii"
metadata[metadata$patho_stage %in% c("stage ia", "stage ib","stage i"),]$patho_stage <- "stage i"
4. 利用function函数出图
这里还有一个问题，我们miR矩阵是以MIRBASE ID命名的(如MITAT00198668)，所以我们从miRBASE获取了miRNA和对应的MIRBASE ID表，以希望输入miRNA名称直接出图。
stage_plot <- function(miRnames){
sample <- as.character(metadata$patho_stage)
mirnames <- as.character(miRBASE[grep(paste0("^",miRnames,"$"),miRBASE$miRNA),]$MIRBASE_ID)
dat <- as.numeric(metadata[,grep(mirnames,colnames(metadata))])
test <- data.frame(cbind(sample,dat))
test <- na.omit(test)
test$sample <- as.factor(test$sample)
test$dat <- as.numeric(test$dat)
test <- test[!test$sample %in% c("","not reported"),]
my_comparisons <- list(c("stage i","stage ii"),c("stage i","stage iii"), c("stage i","stage iv"),c("stage ii","stage iii"),c("stage iii","stage iv"))
test[is.na(test$dat)] <- 0
ggplot(test,aes(x=sample,y=dat,fill=sample),color =sample,line.color = "gray",
line.size = 0.4, palette = "jco")+
geom_boxplot()+
geom_jitter()+
theme_bw()+
stat_compare_means(comparisons=my_comparisons, label.y = c(max(test$dat)+50, max(test$dat)+100,max(test$dat)+150,max(test$dat)+200,max(test$dat)+250))+ # Add pairwise comparisons p-value
stat_compare_means(label.y = max(test$dat)+350) # Add global p-value
ggsave(paste0(miRnames,"_expression_patho_stage.png"),width=6,heigh=6)
}
stage_plot("hsa-miR-123b-5p")
出图结果
根据组织病理级别(Stage)进行的表达差异作图
