拿到基因表达矩阵之后的那点事(一)
做转录组一般拿到基因表达矩阵之后工作即可开始做差异分析，在此之前还有一步就是对矩阵做标准化，常见的几种RPKM、FPKM、TMM等，虽然RPKM、FPKM方法被吐槽的尤为厉害，
但是大多数测序公司给出的结果依然还是很多在使用这种方法，这里我还是以RPKM作为演示
公式如下：
RPKM= read counts / (mapped reads (Millions) * exon length(KB))
其实标准化就是俩点 基因的长度和所有基因上计数后的reads总数
之前看到了一个帖子，感觉说的不错，对基因长度和mapped reads做了介绍，并推荐了一种方法计算这个基因长度点Here
计算基因长度需要用到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)
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)
箱线图看一下表达量分布
#做各个样本箱线图 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()
接下来绘制两两样本散点图
之前都是一个个的去提取两组画图，感觉颇为麻烦，然后看到了一个帖子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 函数这种形式，帖子中还介绍了另一种方法用到everything函数，大家可以去看原贴学习~~
其实转录组还有很多方法直接用原始矩阵去做差异分析，Deseq2,edgeR等等都有自己的标准化方法，大家也可学习！！今天就先介绍到这吧~
