拿到基因表达矩阵之后的那点事(一)

做转录组一般拿到基因表达矩阵之后工作即可开始做差异分析，在此之前还有一步就是对矩阵做标准化，常见的几种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)

图片.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)#指定两组

图片.png

除了!!ensym 函数这种形式，帖子中还介绍了另一种方法用到everything函数，大家可以去看原贴学习~~

其实转录组还有很多方法直接用原始矩阵去做差异分析，Deseq2,edgeR等等都有自己的标准化方法，大家也可学习！！今天就先介绍到这吧~

下面是刚创建个人公众号，会定时更新R、linux、python，组学方面的学习内容，请多多支持呦~~

image