1. 先準(zhǔn)備好外顯子
參考這篇
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("路徑/gencode.v22.annotation.gtf",format="gtf")
#通過exonsBy獲取每個(gè)gene上的所有外顯子的起始位點(diǎn)和終止位點(diǎn),然后用reduce去除掉重疊冗余的部分 #最后計(jì)算長(zhǎng)度
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)
#轉(zhuǎn)置一下
exons_gene_lens1 <- t(exons_gene_lens1)
#合并rawdata count和exon
counts = read.table("路徑/mut-1.txt", sep = "\t", header = T)
2. 因?yàn)槲业膔awdata基因名字是symbol涝桅，所以需要名稱轉(zhuǎn)換一下拜姿，用到了clusterProfiler里的bitr（但是由于bitr庫(kù)里可能不能及時(shí)更新基因list，所以會(huì)有遺漏的冯遂，這點(diǎn)我找不到解決辦法蕊肥，除非手動(dòng)轉(zhuǎn)換）
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. 這里發(fā)現(xiàn)exon里的ENSEMBL名里含有版本號(hào)，也就是有小數(shù)點(diǎn)
#去掉名字里的小數(shù)點(diǎn)
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矩陣?yán)锏那?46列的數(shù)值都除以kb
countdata <- count_with_exonlength$Mut_NPC_1
rpk <- countdata/kb
t(rpk)
#FPKM計(jì)算
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)
這樣就轉(zhuǎn)換好了壁却，我是一個(gè)一個(gè)樣品轉(zhuǎn)換的，也可以寫函數(shù)批量轉(zhuǎn)換裸准，之后下一步進(jìn)行相關(guān)性分析→GSEA分析
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
