1. 先准备好外显子

参考这篇

library(GenomicFeatures) 
txdb <- makeTxDbFromGFF("路径/gencode.v22.annotation.gtf",format="gtf")

#通过exonsBy获取每个gene上的所有外显子的起始位点和终止位点,然后用reduce去除掉重叠冗余的部分 #最后计算长度 
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})

exons_gene_lens1 <- as.data.frame(exons_gene_lens)
#转置一下
exons_gene_lens1 <- t(exons_gene_lens1)

#合并rawdata count和exon
counts = read.table("路径/mut-1.txt", sep = "\t", header = T)

2. 因为我的rawdata基因名字是symbol，所以需要名称转换一下，用到了clusterProfiler里的bitr（但是由于bitr库里可能不能及时更新基因list，所以会有遗漏的，这点我找不到解决办法，除非手动转换）

library(DOSE)
library(org.Hs.eg.db)
library(clusterProfiler)
gene.df<-bitr(counts$SYMBOL, fromType = "SYMBOL" , 
              toType = c("ENSEMBL"),
              OrgDb = org.Hs.eg.db,
              drop = FALSE)
counts_change <- merge(gene.df, counts, all.x = F, all.y = F)
#去掉NA
all <- apply(counts_change, 1, function(x) all(x==0) )
counts_change <- counts_change[!all,]
counts_change <- counts_change[complete.cases(counts_change),]
sum(is.na(counts_change))

3. 这里发现exon里的ENSEMBL名里含有版本号，也就是有小数点

#去掉名字里的小数点
xc <- gsub("\\.(\\.?\\d*)","",rownames(exons_gene_lens1))
rownames(exons_gene_lens1) = xc
colnames(exons_gene_lens1) = "Length"
exons_gene_lens1 <- as.data.frame(exons_gene_lens1)

rownames <- as.data.frame(row.names(exons_gene_lens1))
exons_gene_lens1[,2] <- rownames
colnames(exons_gene_lens1) <- c("Length", "ENSEMBL")

#合并
count_with_exonlength <- merge(counts_change,exons_gene_lens1, by ="ENSEMBL")
write.csv(count_with_exonlength,file="C:/Users/wang/Desktop/wt-1.csv",quote=F)

#先把length除以1000，就是上面公式里说的单位要kb 
kb <- count_with_exonlength$Length/1000
kb
#把新count矩阵里的前546列的数值都除以kb 
countdata <- count_with_exonlength$Mut_NPC_1
rpk <- countdata/kb 
t(rpk) 
#FPKM计算 
fpkm <- t(t(rpk)/colSums(as.data.frame(countdata)) * 10^6) 
head(fpkm)
MUT_1_FPKM <- cbind(count_with_exonlength, fpkm)

write.csv(MUT_1_FPKM,file="C:/Users/wang/Desktop/MUT-1-fpkm.csv",quote=F)

这样就转换好了，我是一个一个样品转换的，也可以写函数批量转换，之后下一步进行相关性分析→GSEA分析