前面分享了FPKM數據該怎么分析？
后來看了一下別人分享的例子，應該把FPKM轉換為TPM后再進行分析。
為什么說FPKM和RPKM都錯了？
簡單小需求：如何將FPKM轉換成TPM

1、數據下載

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)#在調用as.data.frame的時，將stringsAsFactors設置為FALSE可以避免character類型自動轉化為factor類型
# 注意查看下載文件的大小，檢查數據 
  gset <- read.table("GSE115422_BM_EC_9GY_fpkms.txt.gz",
                     sep='\t',
                     quote = "",
                     fill = T, 
                     comment.char = "!",
                     header = T) # 提取表達矩陣       
  save(gset,file="GSE115422_gset.Rdata")   ## 保存到本地


head(gset)
View(gset)

#去重
# a <- gset$gene_name
# 
# gset1 <- gset[a,]
# gset1 <- gset[gset$gene_name,]
# 
# gset2  <- unique(gset1)

#不知這步處理之后，為啥會提示“Mar-01” Mar-02重復

#生信技能樹或者菜鳥團有相關的推文，找推文為啥基因或者探針會出現日期

#除去Mar-01，Mar-02
gset <- gset[!(gset$gene_name == 'Mar-01'),]

gset <- gset[!(gset$gene_name == 'Mar-02'),]

#把第一列添加為行名
rownames(gset) <- gset[,1]

#刪除第一列
gset <- gset[,-1]

2.將FPKM轉換為TPM

exprSet <- gset

fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
tpms <- apply(exprSet,2,fpkmToTpm)
tpms[1:3,]
colSums(tpms)
#輸出結果：
tpms[1:3,]

3、質控和檢測數據質量

View(tpms)

# group_list <- colnames(tpms)
# group_list <- as.factor(group_list)
# 
# group_list <- as.character(group_list)
# 
# library(stringr)
# 
# group_list <- str_split(as.character(group_list),'.', simplify = T)[,1]

group_list <- c(rep('BMC',3),rep('Lin',3),rep('LSK',3),rep('EC',3),rep('X9GY',3))

#表達矩陣數據校正
exprSet <- tpms
#exprSet <- exprSet[,10:15]
exprSet

boxplot(exprSet,las=2)


#數據里面有很多基因在某些樣本中的檢測值為零,除去某個這部分值；

#如果一個基因在三個及三個以上的樣本里面為零，則除去該基因

#apply(dat1,1,function(x){sum(floor(x)==0)>3})

exprSet <- exprSet[!apply(exprSet,1,function(x){sum(floor(x)==0)>3}),]

boxplot(exprSet,las=2)#不在同一水平線上

#limma歸一化
library(limma)
exprSet <- normalizeBetweenArrays(exprSet)
boxplot(exprSet,las=2)

#
exprSet <- log(exprSet+1)##表達值為成千上萬,需要取log值
boxplot(exprSet,las=2)

if(F){
  exprSet <- log(exprSet+1)#應該取log2,完全不在一個水平上面
  boxplot(exprSet,las=2)
  
  #歸一化
  exprSet <- normalizeBetweenArrays(exprSet)
  
  boxplot(exprSet,las=2)###歸一化后，表達值偏離
  
  #數據里面有很多基因在某些樣本中的檢測值為零,除去某個這部分值；
  
  #如果一個基因在三個及三個以上的樣本里面為零，則除去該基因
  
  #apply(dat1,1,function(x){sum(floor(x)==0)>3})
  
  exprSet <- exprSet[!apply(exprSet,1,function(x){sum(floor(x)==0)>3}),]
  
  boxplot(exprSet,outline=FALSE,notch=T,las=2)#不在同一水平線上
  #應該先做該步，后邊再做歸一化
  
}

#PCA看樣本分組情況
## PCA

library(ggfortify)
df=as.data.frame(t(exprSet))
df$group=group_list
png('pca.png',res=120)

pp <- autoplot(prcomp( df[,1:(ncol(df)-1)] ), 
               data=df,
               colour = 'group')+
  theme_bw()
pp
dev.off()



#boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(limma) 
exprSet <- normalizeBetweenArrays(exprSet)
boxplot(exprSet,las=2)##表達值很高，應該取log2,完全不在一個水平上面

#boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
#判斷數據是否需要轉換
exprSet <- log2(exprSet+1)#應該取log2,完全不在一個水平上面

#取EC，x9GY
exprSet <- exprSet[,10:15]

boxplot(exprSet,las=2)

exprSet <- log2(exprSet+1)

boxplot(exprSet,las=2)

exprSet <- normalizeBetweenArrays(exprSet)

boxplot(exprSet,las=2)


## 下面是畫PCA的必須操作，需要看說明書。

dat <- t(exprSet)#畫PCA圖時要求是行名時樣本名，列名時探針名，因此此時需要轉換
dat <- as.data.frame(dat)#將matrix轉換為data.frame
dat <- cbind(dat,group_list) #cbind橫向追加，即將分組信息追加到最后一列
#dat[1:5,1:5]
head(dat)

#BiocManager::install("FactoMineR")
#BiocManager::install("factoextra")
library("FactoMineR")#畫主成分分析圖需要加載這兩個包
library("factoextra") 
# The variable group_list (index = 54676) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#現在dat最后一列是group_list，需要重新賦值給一個dat.pca,這個矩陣是不含有分組信息的

fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = dat$group_list, # color by groups
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)
ggsave('all_samples_PCA.png')


#  ## hclust
colnames(exprSet) <- paste(group_list,1:ncol(exprSet),sep='_')
# Define nodePar
nodePar <- list(lab.cex = 0.4, pch = c(NA, 19),
                cex = 0.7, col = "blue")
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
png('hclust.tif',res=120)

plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
dev.off()
save(exprSet,group_list,file = 'step2-output.Rdata')
#樣品聚類比較好，本次因為本來就是來自不同類型的細胞分開的非常好，
#實際分析中可以刪除某些離群或者混淆的樣本

4、差異分析

#差異分析：

#單分組
if(F){dat <- exprSet
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
bp=function(g){
  library(ggpubr)
  df=data.frame(gene=g,stage=group_list)
  p <- ggboxplot(df, x = "stage", y = "gene",
                 color = "stage", palette = "jco",
                 add = "jitter")
  #  Add p-value
  p + stat_compare_means()
}
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg) 
}

#本例為多分組
group_list <- c(rep('BMC',3),rep('Lin',3),rep('LSK',3),rep('EC',3),rep('X9GY',3))

group <- factor(group_list)
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
design

#指定那類樣本比上那類樣本，特別注意有順序，橫杠前的樣本比上橫杠后的樣本
contrast.matrix <- makeContrasts(X9GY - EC, 
                                 Lin - BMC,
                                 LSK - BMC, 
                                 levels=design)

contrast.matrix


fit <- lmFit(exprSet, design)

fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2)

allDiff1=topTable(fit2,adjust='fdr',coef=1,number=Inf) 
allDiff2=topTable(fit2,adjust='fdr',coef=2,number=Inf) 
allDiff3=topTable(fit2,adjust='fdr',coef=3,number=Inf)

5、差異基因注釋分析

library(ggstatsplot);
library(cowplot);
library(clusterProfiler);
library(stringr);
library(dplyr);
library(tidyr);
library(ggplot2);
library(ggstatsplot);
## 不同的閾值，篩選到的差異基因數量就不一樣，后面的超幾何分布檢驗結果就大相徑庭。
deg <- allDiff1
if(T){
  logFC_t=1.5
  deg$g=ifelse(deg$P.Value>0.05,'stable',
               ifelse( deg$logFC > logFC_t,'UP',
                       ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
  )
  table(deg$g)
  head(deg)
  deg$symbol=rownames(deg)
  library(ggplot2)
  library(clusterProfiler)
  library(org.Mm.eg.db)
  df <- bitr(unique(deg$symbol), fromType = "SYMBOL",
             toType = c( "ENTREZID"),
             OrgDb = org.Mm.eg.db)
  head(df)
  DEG=deg
  head(DEG)
  
  DEG=merge(DEG,df,by.y='SYMBOL',by.x='symbol')
  head(DEG)
  
  save(DEG,file = 'anno_DEG.Rdata')
  gene_up= DEG[DEG$g == 'UP','ENTREZID'] 
  gene_down=DEG[DEG$g == 'DOWN','ENTREZID'] 
}

KEGG分析，超幾何分布檢驗

###這里就拿KEGG數據庫舉例吧，拿自己判定好的上調基因集進行超幾何分布檢驗，如下
if(T){
  gene_down
  gene_up
  enrichKK <- enrichKEGG(gene         =  gene_up,
                         organism     = 'mmu',
                         #universe     = gene_all,
                         pvalueCutoff = 0.05,
                         qvalueCutoff =0.05)
  head(enrichKK)[,1:6] 
  browseKEGG(enrichKK, 'hsa04512')
  dotplot(enrichKK)
  ggsave("enrichKK.png")
  enrichKK=DOSE::setReadable(enrichKK, OrgDb='org.Mm.eg.db',keyType='ENTREZID')
  enrichKK 
}
##最基礎的條形圖和點圖
#條帶圖
barplot(enrichKK,showCategory=20)
#氣泡圖
dotplot(enrichKK)

通路與基因之間的關系可視化

#通路與上調基因之間的關系可視化
###制作genlist三部曲：
## 1.獲取基因logFC
DEG_up <- DEG[DEG$g == 'UP',]
geneList <- DEG_up$logFC
## 2.命名
names(geneList) = DEG_up$ENTREZID
## 3.排序很重要
geneList = sort(geneList, decreasing = TRUE)
head(geneList)

cnetplot(enrichKK, categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(enrichKK, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
ggsave("enrichKK_cnetplot.png")

通路與通路之間的連接展示

#通路與通路之間的連接展示
emapplot(enrichKK)
ggsave("enrichKK_emapplot.png")

熱圖展現通路與基因之間的關系

#熱圖展現通路與基因之間的關系
heatplot(enrichKK)
ggsave("enrichKK_heatplot.png")

GO分析重點是ID轉換

library(clusterProfiler);
ego_bp_up<-enrichGO(gene       = DEG_up$ENTREZID,
                 OrgDb      = org.Mm.eg.db,
                 keyType    = 'ENTREZID',
                 ont        = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,#0.01
                 qvalueCutoff = 0.05)
goplot(ego_up)
ggsave("ego_bp_up_goplot.png")
head(ego)
library(stringr)
barplot(ego_bp_up,showCategory = 16,title="The GO_BP enrichment analysis of all DEGs ")+ 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(ego_bp) str_wrap(ego_bp,width = 25))
ggsave("ego_bp_up_barplot.png")

參考
RNAseq數據，下載GEO中的FPKM文件后該怎么下游分析