序列分析
TCGA及lncRNA矩阵整理
一、目的
二、代码
#这个函数有三个参数，
#@metadata是从TCGA数据下载的sample sheet
#@path是保存样本表达文件的路径
#@data.type是要合并的数据的类型，这里支持RNAseq和miRNAs两种
#@mRNA_expr_type是RNA表达谱的类型，这里支持TPMcounts，TPM，FPKM和FPKM-UQ
#@symbol,为T提取基因名字，为F，不提取
#@type,为T提取RNA类型，为F，不提取
library(limma)
setwd("E:/微生物-赵建组/课题进展/生信挖掘相关基因/03.TCGA胃癌数据/01.TCGA下载的胃癌的数据") #工作目录（需修改）
merge_TCGA <- function(metadata, path, data.type, mRNA_expr_type="TPM", symbol=T, RNA_type=T){
#通过合并path,还有sample sheet前两列得到每一个文件的完整路径
filenames <- file.path(path, metadata$file_id, metadata$file_name,
fsep = .Platform$file.sep)
#判断需要合并的是什么数类型，如果是RNAseq执行下面的代码
if (data.type=='RNAseq') {
message ('############### 正在进行提取，请稍后 ################\n',
' ### 正在进行提取 ###\n',
' ### 正在进行提取 ###\n',
' ### 正在进行提取 ###\n')
#根据需要的RNAseq表达谱数据类型，提取相应列的数据
if(mRNA_expr_type=="TPM"){
column=4
}else if(mRNA_expr_type=="TPM"){
column=7
}else if(mRNA_expr_type=="FPKM"){
column=8
}else if(mRNA_expr_type=="FPKM_UQ"){
column=9
}
#通过lapply循环去读每一个样本的表达，然后通过cbind合并成矩阵
rnaMatrix <- do.call("cbind", lapply(filenames, function(fl)
read.table(fl,skip=6,sep="\t")[,column]))
#获取第一个文件的第一列作为矩阵的行名
ensembl <- read.table(filenames[1],skip=6,sep="\t",stringsAsFactors = F)$V1
gene_symbol <- read.table(filenames[1],skip=6,sep="\t",stringsAsFactors = F)$V2
type <- read.table(filenames[1],skip=6,sep="\t",stringsAsFactors = F)$V3
#排除掉_PAR_Y为后缀的转录本
index=grepl("^\\d+$",sapply(strsplit(ensembl, '.', fixed=TRUE), '[',2))
rnaMatrix=rnaMatrix[index,]
#去掉Ensembl ID后面的.和数字，eg.ENSG00000000003.13
rownames(rnaMatrix) <- sapply(strsplit(ensembl[index], '.', fixed=TRUE), '[',1)
gene_symbol=gene_symbol[index]
type=type[index]
#将sample sheet的sample id这一列作为表达矩阵的列名
colnames(rnaMatrix) <- metadata$sample
#统计样本数和基因数
nSamples = ncol(rnaMatrix)
nGenes = nrow(rnaMatrix)
if(RNA_type){
rnaMatrix=data.frame(type,rnaMatrix,stringsAsFactors = F,check.names = F)
}
if(symbol){
rnaMatrix=data.frame(gene_symbol,rnaMatrix,stringsAsFactors = F,check.names = F)
}
#输出样本数和基因数
message (paste('Number of samples: ', nSamples, '\n', sep=''),
paste('Number of genes: ', nGenes, '\n', sep=''))
#返回最后的基因表达矩阵
return (rnaMatrix)
}else if (data.type=='miRNAs') { #如果需要合并的是miRNA的数据，执行下面代码
message ('############### Merging miRNAs data ###############\n')
#利用lapply来读取每个样本miRNA的表达数据，这里需要用到filtermir这个函数
#filtermir主要提取mature mir的counts数
mirMatrix <- lapply(filenames, function(fl) filtermir(fl))
#mirbase是目前人的说有miRNA成熟提的ID号，eg.MIMAT0027618
mirs <- mirbase$V1
#cbind合并成矩阵合并成矩阵，这里注意并不是每个样本中都表达所有的miRNA
#如果不存在某个miRNA,在这一步表达量会用NA表示
mirMatrix <- do.call('cbind', lapply(mirMatrix,
function(expr) expr[mirs]))
#设置表达矩阵的行名为miRNA的名字，mirbase的第二列就是所有人的miRNA的名字
rownames(mirMatrix) <- mirbase$V2
#设置表达矩阵的列名为sample sheet的sample id这一列
colnames(mirMatrix) <- metadata$sample
#将表达量为NA的地方转换成0
mirMatrix[is.na(mirMatrix)] <- 0
#统计样本数和miRNA数目
nSamples = ncol(mirMatrix)
nGenes = nrow(mirMatrix)
#输出样本数和miRNA的数目
message (paste('Number of samples: ', nSamples, '\n', sep=''),
paste('Number of miRNAs: ', nGenes, '\n', sep=''))
#返回miRNA的表达矩阵
return (mirMatrix)
}else{ #如果data.type不是上面提到的两种，就报错，停止执行
stop('data type error!')
}
}
#定义filtermir函数
#@fl参数为所有样本miRNA表达文件的链接
filtermir <- function(fl) {
#readtable读取文件内容
expr <- read.table(fl, header=TRUE, stringsAsFactors = FALSE)
#找到miRNA_region这一列以mature开头的行，eg.mature,MIMAT0000062
expr <- expr[TPMtsWith(expr$miRNA_region, "mature"),]
#将同一个成熟体的所有counts累加起来，作为这个成熟体counts
expr <- aggregate(expr$read_count, list(expr$miRNA_region), sum)
#将miRNA_region这一列（mature,MIMAT0000062）按头号分开，取第二列，得到成熟体ID号
mirs <- sapply(strsplit(expr$Group.1, ',', fixed=TRUE),'[',2)
#去掉第一列，只保留miRNA成熟体的counts数
expr <- expr[,-1]
#加上对应的成熟体ID作为counts的名字
names(expr) <- mirs
#返回过滤之后的miRNA成熟体的counts
return(expr)
}
#定义去除重复样本的函数FilterDuplicate
FilterDuplicate <- function(metadata) {
filter <- which(duplicated(metadata[,'sample']))
if (length(filter) != 0) {
metadata <- metadata[-filter,]
}
message (paste('Removed', length(filter), 'samples', sep=' '))
return (metadata)
}
#定义过滤样本类型的函数FilterSampleType，只保留PrimaryTumor和SolidTissueNormal这两种样本类型
FilterSampleType <- function(metadata) {
filter <- which(! metadata$sample_type %in%
c('PrimaryTumor', 'SolidTissueNormal'))
if (length(filter) != 0) {
metadata <- metadata[-filter,]
}
message (paste('Removed', length(filter), 'samples', sep=' '))
return (metadata)
}
#读入RNAseq的sample sheet文件
metaMatrix.RNA=read.table("RNAseq_sample_sheet.tsv",sep="\t",header=T)
#替换.为下划线，转换成小写，sample_id替换成sample
names(metaMatrix.RNA)=gsub("sample_id","sample",gsub("\\.","_",tolower(names(metaMatrix.RNA))))
#删掉最后一列sample_type中的空格
metaMatrix.RNA$sample_type=gsub(" ","",metaMatrix.RNA$sample_type)
#删掉重复的样本
metaMatrix.RNA <- FilterDuplicate(metaMatrix.RNA)
#删掉非PrimaryTumor和SolidTissueNormal的样本
#会删掉像Recurrent Tumor这样的样本
metaMatrix.RNA <- FilterSampleType(metaMatrix.RNA)
#调用merge_TCGA函数合并RNAseq的表达矩阵
#TPM
RNA_TPM=merge_TCGA(metadata=metaMatrix.RNA,
path="RNAseq",
data.type="RNAseq",
mRNA_expr_type="TPM",
symbol = T,
RNA_type=T
)
RNA_TPM=as.matrix(RNA_TPM)
rownames(RNA_TPM)=RNA_TPM[,1]
RNA_TPM_exp=RNA_TPM[,3:ncol(RNA_TPM)]
RNA_TPM_dimnames=list(rownames(RNA_TPM_exp),colnames(RNA_TPM_exp))
RNA_TPM_data=matrix(as.numeric(as.matrix(RNA_TPM_exp)),nrow=nrow(RNA_TPM_exp),dimnames=RNA_TPM_dimnames)
RNA_TPM_data=avereps(RNA_TPM_data)
write.table(file="combined_RNAseq_TPM.txt",RNA_TPM_data,sep="\t",quote=F)
#提取lncRNA
RNA_TPM=as.data.frame(RNA_TPM)
RNA_TPM_lnc=RNA_TPM[RNA_TPM$type=="lncRNA",]
RNA_TPM_lnc=as.matrix(RNA_TPM_lnc)
rownames(RNA_TPM_lnc)=RNA_TPM_lnc[,1]
RNA_TPM_lnc_exp=RNA_TPM_lnc[,3:ncol(RNA_TPM_lnc)]
RNA_TPM_lnc_dimnames=list(rownames(RNA_TPM_lnc_exp),colnames(RNA_TPM_lnc_exp))
RNA_TPM_lnc_data=matrix(as.numeric(as.matrix(RNA_TPM_lnc_exp)),nrow=nrow(RNA_TPM_lnc_exp),dimnames=RNA_TPM_lnc_dimnames)
RNA_TPM_lnc_data=avereps(RNA_TPM_lnc_data)
write.table(file="combined_RNAseq_TPM_lncRNA.txt",RNA_TPM_lnc_data,sep="\t",quote=F)