做轉錄組一般拿到基因表達矩陣之后工作即可開始做差異分析鲫趁，在此之前還有一步就是對矩陣做標準化确沸，常見的幾種RPKM丧凤、FPKM募闲、TMM等，雖然RPKM愿待、FPKM方法被吐槽的尤為厲害浩螺，但是大多數(shù)測序公司給出的結果依然還是很多在使用這種方法，這里我還是以RPKM作為演示
公式如下：
RPKM= read counts / (mapped reads (Millions) * exon length(KB))
其實標準化就是倆點 基因的長度和所有基因上計數(shù)后的reads總數(shù)
之前看到了一個帖子仍侥，感覺說的不錯要出，對基因長度和mapped reads做了介紹，并推薦了一種方法計算這個基因長度點Here
計算基因長度需要用到gtf文件,代碼如下:
# 根據(jù)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)#計算各樣本總數(shù)
#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)
OK,簡單做個熱圖瞅一瞅
#相關性熱圖
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到了新方法否淤，還是用函數(shù)去解決，參照原貼
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 函數(shù)才可以
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)#指定兩組
圖片.png
除了!!ensym 函數(shù)這種形式棠隐，帖子中還介紹了另一種方法用到everything函數(shù)石抡，大家可以去看原貼學習~~
其實轉錄組還有很多方法直接用原始矩陣去做差異分析，Deseq2,edgeR等等都有自己的標準化方法助泽，大家也可學習浮！今天就先介紹到這吧~
下面是剛創(chuàng)建個人公眾號嗡贺，會定時更新R隐解、linux良风、python檬洞，組學方面的學習內容，請多多支持呦~~
image
