高通量测序数据一般公司都会提供两种矩阵,一种是Row counts,前面说过的用于差异基因的筛选。第二种是FPKM值,可以理解为转录组基因的表达矩阵,可以用于做热图和基因表达变化的比较。但是数据挖掘中,我们只能下载到counts矩阵,所以要得到基因的表达量还需要通过计算。 这里我们使用biomaRt包举个例子,由counts值计算FPKM。 一、FPKM的计算 加载counts矩阵: setwd("E:/生物信息学") ma <- read.csv("GSE169758_markdup.featco.2.counts.csv",header = T,row.names = 1) 下载安装R包: BiocManager::install("biomaRt") library(biomaRt) 链接到ensembl数据库,并选择合适的物种,因为我们下载的是人的数据,所以dataset = "hsapiens_gene_ensembl",小鼠的则选择dataset = “mmusculus_gene_ensembl”。计算基因长度。 ensembl <- useMart("ensembl") ensembl = useDataset("hsapiens_gene_ensembl",mart=ensembl) test <- getBM(attributes=c('ensembl_gene_id', 'start_position', 'end_position','ensembl_transcript_id', 'transcript_length'),mart = ensembl) head(test) # ensembl_gene_id start_position end_position ensembl_transcript_id transcript_length #1 ENSG00000210049 577 647 ENST00000387314 71 #2 ENSG00000211459 648 1601 ENST00000389680 954 #3 ENSG00000210077 1602 1670 ENST00000387342 69 #4 ENSG00000210082 1671 3229 ENST00000387347 1559 #5 ENSG00000209082 3230 3304 ENST00000386347 75 #6 ENSG00000198888 3307 4262 ENST00000361390 956 将得到的test文件排序,并去除重复的ID: test <- test[order(test$ensembl_gene_id,test$transcript_length,decreasing = T),] g <- test[!duplicated(test$ensembl_gene_id),] g <- g[,c(1,5)] head(g) dim(g) summary(g) 取基因长度文件和count文件交集的基因: ng=intersect(rownames(ma),g$ensembl_gene_id) length(ng) lengths=g[match(ng,g$ensembl_gene_id),2] names(lengths) <- g[match(ng,g$ensembl_gene_id),1] head(lengths) 根据最终选取的counts矩阵计算每个样本的文库大小: ma <- ma[names(lengths),] total_count <- colSums(ma) head(total_count) #Mcc1 Mcc2 Mcc3 Mcc4 Mcc5 Mcc6 #39315944 15971423 16166354 25798072 33085950 1772396 计算FPKM,并保存文件: ma_fpkm <- t(do.call( rbind, lapply(1:length(total_count), function(i){ 10^9*ma[,i]/lengths/total_count[i] #lengths向量自动遍历 }) )) ma_fpkm[1:4,1:4] #[,1] [,2] [,3] [,4] #ENSG00000000003 0.01340093 0.000000 0.000000 0.01021143 #ENSG00000000005 0.00000000 0.000000 0.000000 0.00000000 #ENSG00000000419 28.87258994 24.594527 23.186029 27.78157422 #ENSG00000000457 2.23382616 1.677457 1.696455 2.31666061 write.csv(ma_fpkm, file = "ma_fpkm.csv") 二、基因ID转gene symbol 从前面的结果我们很清楚的看到,每个基因都是用ID表示的,类似于ENSG00000000003,但看这个名称,根本不知道它是什么,我们通常也不习惯这样表示基因,还是习惯于Gene symbol。接下来我们将这个FPKM文件进行基因ID注释。 首先加载人基因数据库: library(stringr) BiocManager::install("org.Hs.eg.db") library(org.Hs.eg.db) 准备三个文件,一个是gene ID,一个是gene symbol,还有一个是我们需要注释的文件,准备一列ensembl_id: gs <- toTable(org.Hs.egSYMBOL) head(gs) #gene_id symbol #1 1 A1BG #2 2 A2M #3 3 A2MP1 #4 9 NAT1 #5 10 NAT2 #6 11 NATP ge <- toTable(org.Hs.egENSEMBL) head(ge) #gene_id ensembl_id #1 1 ENSG00000121410 #2 2 ENSG00000175899 #3 3 ENSG00000256069 #4 9 ENSG00000171428 #5 10 ENSG00000156006 #6 12 ENSG00000196136 colnames(ma_fpkm) <- colnames(ma) ma_fpkm <- as.data.frame(ma_fpkm) ma_fpkm$ensembl_id <- row.names(ma_fpkm) 接着就是将这几个文件一一比对merge,就可以得到注释的gene symbol: b <- merge(ma_fpkm,ge,by="ensembl_id",all.x=T) c <- merge(b,gs,by="gene_id",all.x=T) c=c[order(c$ensembl_id),] c=c[!duplicated(c$ensembl_id),] c=c[match(ma_fpkm$ensembl_id,c$ensembl_id),] head(c) # Pan1 Pan2 Pan3 Pan4 Pan5 Pan6 symbol #29657 4.55462494 2.612894 3.31715074 3.80483778 4.487819 2.49863368 TSPAN6 #27673 0.00000000 0.000000 0.00000000 0.02566601 0.000000 0.00000000 TNMD #33278 27.28248482 18.442794 22.18114506 23.02063714 25.941927 20.66173442 DPM1 #26324 1.56335564 2.044088 2.33643810 2.39752223 1.953669 2.30114296 SCYL3 #25673 4.02781699 4.294728 5.45416997 4.60184848 3.717569 5.34177442 C1orf112 #15388 0.01591371 0.000000 0.01808745 0.01172831 0.000000 0.01215142 FGR write.csv(c,file = "genesymbol_fpkm.csv") 这样我们的问题就解决了! |
|