需要‘All_hg19gene_len.txt’文件者私信留言
#counts值轉(zhuǎn)化為TPM，F(xiàn)PKM
#讀入length
len=read.table('All_hg19gene_len.txt',header = T,sep = '\t')
colnames(len)[1] = 'SYMBOL'
# read expression Matter
exp=read.csv('exp-ID-counts.csv')
head(exp)
#鏈接
library(dplyr)
merge<-left_join(exp,len,by="SYMBOL")#根據(jù)基因那列進(jìn)行合并
merge <- na.omit(merge)#刪除錯(cuò)誤值行
write.csv(merge,file = "merge.csv",sep = "\t")#讀出文件琴昆，直接往下運(yùn)行也許
#計(jì)算TPM
head(merge)
merge1 <- merge[!duplicated(merge$SYMBOL),]
rownames(merge1)<-merge1$SYMBOL
merge1<-merge1[,-1]
head(merge2)#最后一列Length是基因長(zhǎng)度
merge2=merge1[,-c(95,96)]
kb <- merge2$Length/1000
kb
countdata <- merge2[,1:94]
rpk <- countdata / kb
rpk
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
head(tpm)
write.table(tpm,file="tpm.xls",sep="\t",quote=F)
##計(jì)算FPKM
fpkm <- t(t(rpk)/colSums(countdata) * 10^6)
head(fpkm)
write.table(fpkm,file="fpkm.xls",sep="\t",quote=F)
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
