高通量測序數(shù)據一般公司都會提供兩種矩陣，一種是Row counts，前面說過的用于差異基因的篩選。第二種是FPKM值，可以理解為轉錄組基因的表達矩陣，可以用于做熱圖和基因表達變化的比較。但是數(shù)據挖掘中，我們只能下載到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數(shù)據庫，并選擇合適的物種，因為我們下載的是人的數(shù)據，所以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矩陣計算每個樣本的文庫大?。?/p>
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注釋。
首先加載人基因數(shù)據庫：
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")
這樣我們的問題就解決了！
接下來我們會著重于轉錄組數(shù)據的可視化和一些后續(xù)的分析！歡迎關注分享和交流討論！
