https://cloud.tencent.com/developer/article/1906294
####################定义基因长度为 非冗余exon长度之和
library("TxDb.Hsapiens.UCSC.hg38.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
if(!file.exists('gene_length.Rdata')){
exon_txdb=exons(txdb)
genes_txdb=genes(txdb)
o = findOverlaps(exon_txdb,genes_txdb)
o
t1=exon_txdb[queryHits(o)]
t2=genes_txdb[subjectHits(o)]
t1=as.data.frame(t1)
t1$geneid=mcols(t2)[,1]
g_l = lapply(split(t1,t1$geneid),function(x){
# x=split(t1,t1$geneid)[[1]]
head(x)
tmp=apply(x,1,function(y){
y[2]:y[3]
})
length(unique(unlist(tmp)))
})
head(g_l)
g_l=data.frame(gene_id=names(g_l),length=as.numeric(g_l))
save(g_l,file = 'gene_length.Rdata')
}
#检查表达量矩阵额基因长度信息
load(file = 'gene_length.Rdata')
head(g_l)
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
genelength <- merge(g_l,g2s,by='gene_id')
head(genelength)
genelength=genelength[,c(3,2)]
colnames(genelength) <- c("gene","length")
#排序
counts <- exp[rownames(exp) %in% genelength$gene,]
#转换成矩阵
count<- as.matrix(counts)
#存贮基因长度信息
genelength[match(rownames(counts),
genelength$gene),"length"]
# eff_length<- genelength$length
eff_length <- genelength[match(rownames(counts),genelength$gene),"length"]/1000
class(count[1,1])
#count转为数值
data=apply(count,2,as.numeric)
row.names(data)<-row.names(count)
count<-data
x <- count/eff_length
mat_tpm <- t( t(x) / colSums(x) ) * 1e6
rownames(mat_tpm)<- rownames(count)
mat_tpm[1:4,1:4]
#mat_tpm<-mat_tpm[-1,]
write.table(mat_tpm, file="paca.tidy.tpm.txt", sep="\t", col.names=T, quote=F,row.names = T)
