根据samtool 统计的reads数计算tpm值

library(data.table)
library(dplyr)
#DNAtpm计算
SMP<-c(123,98,126,4,128,32,134,8,130,29,131,42,119,19,139,24,125,37,144,14,117,5,142,13,127,22,129,48,138,47,140,20,132,18,149,73,120,105,148,44,122,30,124,26,143,34,133,12,145,21,152,10,118,102,147,80,121,1,151,35,135,38,136,2,150,17,146,55,141,11,137,97)
a<-fread("123.DNA.abd",header=F)
a[,1:2]->a #获得表头
colnames(a)<-c("ORF","length")
for (i in SMP){
fread(paste(i,".DNA.abd",sep=""),header=F)->b
left_join(a,b[,c(1,3)],by=c("ORF"="V1"))->a
}
colnames(a)[3:(length(SMP)+2)]<-SMP
a[-which(a[,1]=="*"),]->a
write.table(a,"DNA.reads",sep="\t",quote=F,row.names=F)
t(t(1e6*a[,3:(length(SMP)+2)]/t(a[,2]))/colSums(a[,3:(length(SMP)+2)]/t(a[,2])))->c
data.frame(a[,1],c)->c
c[-which(rowSums(c[,-1])==0),]->d
write.table(d,file="DNAtpm.tsv",quote=F,sep="\t",row.names=F)

#RNAtpm
SMP<-c(123,98,126,107,128,32,134,8,130,29,131,42,119,108,139,24,125,37,144,83,117,5,142,74,127,22,129,81,138,99,140,20,132,82,149,73,120,105,148,44,122,101,124,100,143,76,133,12,145,79,152,77,118,102,147,80,121,106,151,104,135,38,136,2,150,78,146,55,141,11,137,97)
a<-fread("123.RNA.abd",header=F)
a[,1:2]->a
colnames(a)<-c("ORF","length")
for (i in SMP){
fread(paste(i,".RNA.abd",sep=""),header=F)->b
left_join(a,b[,c(1,3)],by=c("ORF"="V1"))->a
}
colnames(a)[3:(length(SMP)+2)]<-SMP
a[-which(a[,1]=="*"),]->a
write.table(a,"RNA.reads",sep="\t",quote=F,row.names=F)
t(t(1e6*a[,3:(length(SMP)+2)]/t(a[,2]))/colSums(a[,3:(length(SMP)+2)]/t(a[,2])))->c
data.frame(a[,1],c)->c
c[-which(rowSums(c[,-1])==0),]->d
write.table(d,file="RNAtpm.tsv",quote=F,sep="\t",row.names=F)