在做转录组分析,下载的数据有各种形式,比如count,FPKM,tpm。我今年2月份发表的一篇文章,小修审稿人提出FPKM过时了,建议我换成CPM。
我特地去pubmed查了查,发现CPM确实比FPKM和TPM要好。但是现在大多数的文章还是支持TPM格式,但是FPKM已经过时了。
我会陆续将TCGA网站上的数据转换成CPM和TPM格式然后上传到我的资源,这样大家就避免了数据前期处理的麻烦。
下面我们开始进入正题:
今天先讲FPKM转变为TPM,很多数据都是FPKM格式,例如TCGA里面的表达数据,除了count,就是FPKM的形式。
下面我们以TCGA的数据为例子:
选择一个胆管癌数据库:
选择FPKM格式的测序数据:
下载数据:
此处注意数据的格式是log2(fpkm+1)的形式:
我们要做fpkm转换为tpm,需要将表达数据转换成fpkm,即2^data-1。
我们先读取数据:
setwd("H:\\")
dir()
data <- read.csv("TCGA-CHOL.htseq_fpkm.tsv",header = T,sep = "\t")
rownames <- as.data.frame(data$Ensembl_ID)## 将第一列的基因名传给rownames
data <- data[,-1]## 减去第一列基因名,只剩下表达数据
data[1:5,1:5]
将log2(FPKM+1)转换成FPKM形式:
data <- 2^data-1
data[1:5,1:5]
> data[1:5,1:5]
TCGA.ZD.A8I3.01A TCGA.W5.AA2U.11A TCGA.W5.AA30.01A TCGA.W5.AA38.01A TCGA.W5.AA30.11A
1 0.06480644 0.0000000 0.000000 0.000000 0.0000000
2 0.00000000 0.0000000 0.000000 0.000000 0.0000000
3 2.77571417 1.1098491 3.140476 5.657167 1.6098847
4 0.00000000 0.0000000 0.000000 0.000000 0.0000000
5 4.32271632 0.5173048 2.536254 12.805596 0.5115492
FPKM转换成TPM:
fpkmToTpm <- function(fpkm)
{
exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
TPM <- apply(data,2,fpkmToTpm)
可以看到,每个样本名的总和都是1*1^6,说明已经完成tpm转换。
colSums(TPM)
> colSums(TPM)
TCGA.ZD.A8I3.01A TCGA.W5.AA2U.11A TCGA.W5.AA30.01A TCGA.W5.AA38.01A TCGA.W5.AA30.11A TCGA.YR.A95A.01A
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
TCGA.W5.AA36.01A TCGA.W5.AA2U.01A TCGA.ZH.A8Y2.01A TCGA.3X.AAVB.01A TCGA.W5.AA2Q.01A TCGA.W5.AA2Z.01A
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
TCGA.3X.AAVA.01A TCGA.W5.AA34.11A TCGA.ZH.A8Y6.01A TCGA.ZH.A8Y1.01A TCGA.ZU.A8S4.01A TCGA.W5.AA2R.01A
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
TCGA.W5.AA2H.01A TCGA.W5.AA2R.11A TCGA.W5.AA2I.11A TCGA.W5.AA39.01A TCGA.W5.AA2O.01A TCGA.3X.AAVC.01A
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
TCGA.3X.AAVE.01A TCGA.W5.AA2X.01A TCGA.ZH.A8Y5.01A TCGA.W5.AA33.01A TCGA.ZH.A8Y8.01A TCGA.WD.A7RX.01A
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
TCGA.W5.AA34.01A TCGA.W6.AA0S.01A TCGA.4G.AAZO.01A TCGA.W5.AA31.01A TCGA.W5.AA2W.01A TCGA.W5.AA2Q.11A
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
TCGA.W5.AA31.11A TCGA.ZH.A8Y4.01A TCGA.W5.AA2T.01A TCGA.W5.AA2X.11A TCGA.4G.AAZT.01A TCGA.3X.AAV9.01A
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
TCGA.ZU.A8S4.11A TCGA.W5.AA2I.01A TCGA.W5.AA2G.01A
1e+06 1e+06 1e+06
转换完以后,此时是TPM格式,记住将TPM转变成log2(TPM +1),然后再和之前保存的基因rownames合并。
TPM <- log2(TPM+1)
data <- cbind(rownames,TPM)
write.csv(data,"data.csv")## 读出数据