RNAseq下遊分析(一)--標準化及簡單可視化

前言

做轉錄組一般拿到基因表達矩陣之後工作即可開始做差異分析，在此之前還有一步就是對矩陣做標準化，常見的幾種RPKM、FPKM、TMM等，雖然RPKM、FPKM方法被吐槽的尤為厲害，但是大多數測序公司給出的結果依然還是很多在使用這種方法,同樣網上對各種標準化的方法也各有見解，如何選擇建議大家多去google了解每種方法的優劣，適用於何種情形，根據自己的數據選擇最合適的標準化方法，比如TPM這種方法會對各個樣品的總基因表達量均一化處理使其一致，會更適合觀察哪個基因在哪個組織裡高表達這種類型的實驗~

RPKM和FPKM的公式一致為read counts / (mapped reads (Millions) * exon length(KB)),可以看到其實標準化就是倆點 ，length指的就是基因的長度，mapped reads指的是所有基因上計數後的reads總數,TPM的方法和這兩種方法的最大區別就是長度和reads總數彼標準化順序先後的問題，之前看到了一個帖子(參考連結1)，解釋的很清晰可以，對基因長度和mapped reads這倆點的一些常見誤區做了介紹，大家可去原貼瀏覽，且文中推薦了一種新的方法計算這個基因長度

計算基因長度需要用到gtf文件,代碼如下:

# 根據gtf文件求基因長度，載入GTF文件，其實和feature結果一致
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("hg38.gtf",format="gtf")
#通過exonsBy獲取每個gene上的所有外顯子的起始位點和終止位點，然後用reduce去除掉重疊冗餘的部分，最後計算長度
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
#獲得每個基因下所有外顯子的總長度後，就可以利用上述公式計算FPKM了。
exons_gene_lens <- as.matrix(exons_gene_lens)
write.table(exons_gene_lens, 'length.txt', sep = '\t', col.names = NA, quote = FALSE)
#通過exonsBy獲取每個gene上的所有外顯子的起始位點和終止位點，然後用reduce去除掉重疊冗餘的部分，最後計算長度
RPKM標準化
#導入基因長度文件
leng <- read.delim('length.txt',header = T,stringsAsFactors = FALSE,check.names = FALSE)
#導入counts矩陣
count <- read.table('RPKM.txt',header = T,row.names = 1,stringsAsFactors = FALSE,check.names = FALSE)
#將其順序一致
count <- count[match(leng$Geneid,rownames(count)),]
length <- leng[,2]#提取長度
total_count<- colSums(count)#計算各樣本總數
#RPKM標準化
rpkm <- t(do.call( rbind,
                   lapply(1:length(total_count),
                          function(i){
                            10^9*count[,i]/length/total_count[i]
                          }) ))
dim(na.omit(rpkm))
colnames(rpkm) <- colnames(count)
rownames(rpkm) <- leng$Geneid
rpkm <- as.data.frame(rpkm)
熱圖
#相關性熱圖  
library(RColorBrewer)
rpkm <- log2(rpkm+1) #log一下
cor_pearson <- cor(rpkm, method = 'pearson')
coul <- colorRampPalette(brewer.pal(8, "PiYG"))(25)
pheatmap::pheatmap(cor_pearson,display_numbers = T,color = coul)
圖片.png箱線圖
#做各個樣本箱線圖
tiff('boxplot.tiff', width = 1000, height = 600)
#設置分組信息 方便作圖
group_list <- c(rep('T',2),rep('N',2),rep('a',2),rep('b',2),rep('c',2),rep('d',2))
group_list <- factor(group_list,levels = c('T','N','a','b','c','d'),ordered=TRUE) #設置一下分組，因為我有12個樣品，6個組.
boxplot(rpkm,outline = FALSE,border  = group_list,las =2,varwidth = T,
        main = 'RPKM counts',ylab = 'log(RPKM+1)')
dev.off()
圖片.png線性回歸擬合
對兩兩樣品之間我們可以做線性回歸散點圖，之前都是一個個的去提取兩組畫圖，感覺頗為麻煩，然後看到了一個帖子get到了新方法，還是用函數去解決更為方便~
scaplot_r2 <- function(indata,inx,iny){
  nms <- names(indata)
  x <- nms[inx]
  y <- nms[iny]
  regression <- paste0(x,"~",y)
  dat.lm <- lm(as.formula(regression),data = indata)
  r2 <- sprintf("italic(r) == %.4f",sqrt(summary(dat.lm)$r.squared)) 
  labels <-  data.frame(r = r2,stringsAsFactors = FALSE)
  #注意此處需要加上 !!ensym 函數才可以
  p3 <- ggplot(indata, aes(x = !!ensym(x), y = !!ensym(y))) +
    geom_point(colour = "black", size = 1) +
    theme_light() +
    stat_smooth(method = lm, se = FALSE,colour = 'pink') +
    labs(x=paste0(x," (log2 intensity)"),y=paste0(y," (log2 intensity)"))+
    geom_text(data=labels,mapping=aes(x = 6,y=1,label=r2),parse = TRUE,inherit.aes = FALSE,size = 6)
  #  annotate("text", label = r2, x = 5.0, y = 1)
  ggsave(paste(x,'vs',y,'.tiff',sep = ''),p3, width = 6, height = 5)
}
scaplot_r2(rpkm,3,7)#指定兩組
除了!!ensym 函數這種形式，下方參考連結2中還介紹了另一種方法用到everything函數，大家可以去看原貼學習~~其實轉錄組還有很多方法直接用原始矩陣去做差異分析，Deseq2,edgeR等等都有自己的標準化方法，大家也可學習！！今天就先介紹到這吧~
參考連結
1.https://www.bioinfo-scrounger.com/archives/342/ 
2.https://www.jianshu.com/p/393e92aaca68