一般来说,TCGA差异分析或者其他分析时,使用的是包含所有肿瘤样本和正常样本。但是,一般情况下,多数类型的TCGA肿瘤类型是肿瘤远远多于正常样本数量。这是因为有的肿瘤样本配对的有癌旁正常样本,但是大多数是没有配对的正常样本的。因此,这里想要通过代码,十分简便的取出仅包含配对样本们的表达矩阵。 代码比较简单,其中值得注意的是,有时可能会出现,仅有正常样本却没有配对的肿瘤样本的情况,这个时候也要将多余的正常样本删除。
#######################################################
##
#######################################################
rm(list=ls())
load( "mRNA_exprSet.Rda")
metadata <- data.frame(colnames(mRNA_exprSet ))
colnames(metadata) <- 'barcode'
for (i in 1:length(metadata[,1])) {
num <- as.numeric(as.character(substring(metadata[i,'barcode'],14,15)))
if (num == 1 ) {metadata[i,2] <- "T"}
if (num != 1) {metadata[i,2] <- "N"}
}
names(metadata)[2] <- 'Barcode'
table(metadata$Barcode)
N <- subset(metadata, metadata$Barcode == 'N')
N$barcode <- substr(x=N$barcode, start = 1, stop = 12)
#######################################################
##
#######################################################
dt <- as.data.frame(t(mRNA_exprSet))
dt[1:3,1:3]
dt$Barcode <- rownames(dt)
dt$sample <- substr(x=dt$Barcode, start = 1, stop = 12)
dt <- subset(dt, select=c('Barcode', 'sample'))
dt <- dt[which(dt$sample %in% N$barcode),]
table(dt$sample)
df <- as.data.frame(table(dt$sample))
df <- subset(df, df$Freq > 1)
dt <- dt[which(dt$sample %in% df$Var1),]
mRNA_exprSet <- mRNA_exprSet[,which(colnames(mRNA_exprSet) %in% dt$Barcode)]
#######################################################
##
#######################################################
metadata <- data.frame(colnames(mRNA_exprSet ))
colnames(metadata) <- 'barcode'
for (i in 1:length(metadata[,1])) {
num <- as.numeric(as.character(substring(metadata[i,'barcode'],14,15)))
if (num == 1 ) {metadata[i,2] <- "T"}
if (num != 1) {metadata[i,2] <- "N"}
}
names(metadata)[2] <- 'Barcode'
table(metadata$Barcode)
select_colname <- function(expr,name){
expr <- expr[,which(colnames(expr) %in% name[,1])]
expr <- expr[,order(names(expr))]
name <- name[which(name[,1] %in% colnames(expr)),]
name <- name[order(name[,1]),]
expr_name <- list( expr,name)
return( expr_name)
}
datexpr <- select_colname(mRNA_exprSet ,metadata)
data <- datexpr[[1]]
metadata <- datexpr[[2]]
group <- metadata$Barcode
exprSet <- data
save(data, file = 'exprSet.Rdata')
save(group , file = 'group.Rdata')
##########################################################################################
##
###########################################################################################