整理分享最近转录组下游分析用到的代码，包括差异基因、GO/KEGG富集、GSEA、GSVA和常用的作图，主要是在生信技能树Jimmy老师代码的基础上根据自己的习惯修改了一些，感谢Jimmy老师的海量资源和无私奉献，我的分享也是希望自己记录的同时能帮助到有需要的小伙伴。
原始counts数据的检查，这一步主要是画各种图检查定量后样本counts分布，通过柱状图、小提琴图、PCA图多角度判断，数据的整体、组间、批次效应的情况。
加载示例数据
library(airway)
data(airway)
exprset <- assay(airway)
group_list <- c(rep(‘tr’,4),rep('cont',4))
数据变形
me <- function(exprset,group){ ####melt expression####
library(reshape2)
exprsetl <- log2(exprset+1) #####取log值归一化####
exprsetl <- melt(exprsetl)
colnames(exprsetl) <- c('probe','sample','value')
exprsetl$group <- rep(group,each=nrow(exprset))
return(exprsetl)
}
exprsetl <- me(exprset,group_list)
柱状图等
library(ggplot2)
ggplot(exprsetl,aes(x=sample,y=value,fill=group))+geom_boxplot()
#ggplot(exprsetl,aes(x=sample,y=value,fill=group))+geom_violin()
ggplot(exprsetl,aes(value,col=group))+geom_density()
ggplot(exprsetl,aes(value,col=sample))+geom_density()
ggplot(exprsetl,aes(value,col=group))+geom_histogram()
对比组间个别基因表达 如GAPDH等管家基因
gene_boxplot <- function(exprn,gene){
library(ggpubr)
library(reshape2)
exprt <- exprn[grep(gene,exprn$probe),] ####前面取过log值######
exprtg <- data.frame(exprt)
rownames(exprtg) <- rownames(exprt)
exprtg$group <- pd2$type
colnames(exprtg)[1] <- gene
ggboxplot(exprtg, x = "group", y = 'value',fill='group',
palette = "jco",color = 'black',
add = "jitter") +
stat_compare_means(method = 't.test') + ggtitle('GAPDH') ##加p值##
}
gene_boxplot(exprsetl,'GAPDH') ###GAPDH####
hcplot <- function(expr,group_list){
exprt <- expr
colnames(exprt) <- paste0(group_list,seq(1:ncol(exprt)),sep='')
hc <- hclust(dist(t(exprt)))
nodePar <- list(lab.cex=0.6,pch=c(NA,19),cex=1,col='blue')
#plot(as.dendrogram(hc))
plot(as.dendrogram(hc),nodePar=nodePar,horiz=T)
}
p <- hcplot(exprset,group_list)
PCA 主成分分析，如检查批次或组别
point PCA 简易版
pcap <- function(expr,group_list){
library(ggfortify)
df <- as.data.frame(t(exprset))
df$group <- as.character(group_list)
autoplot(prcomp(df[,1:(ncol(df)-1)]),data=df,colour='group')
}
pcap(exprset,group_list)
sample PCA 显示分组
pcae <- function(expr,group_list,pdata){
library("FactoMineR")
library("factoextra")
names(group_list) <- rownames(pdata) ###or pdata$batch####
df.pca <- PCA(t(exprset),graph = T)
fviz_pca_ind(df.pca,addEllipses = T,col.ind =group_list)
#fviz_pca_ind(df.pca,addEllipses = T,col.ind =factor(pdata$batch))
}
pcae(exprset,group_list,pdata)
if(T){
library(corrplot)
library(pheatmap)
corrplot(cor(log2(exprset+1)))
c <- cor(log2(exprset+1))
pheatmap(scale(c))
pheatmap(c,scale = 'none')
}
#######or filter top500 cpm#####
if(T){
library(pheatmap)
exprset=exprset[apply(exprset,1, function(x) sum(x>1) > 1),] ##filter####
dim(exprset)
exprset=log(edgeR::cpm(exprset)+1)
dim(exprset)
exprset=exprset[names(sort(apply(exprset, 1,mad),decreasing = T)[1:500]),]
dim(exprset)
M=cor(log2(exprset+1))
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(exprset)
corf=pheatmap(M,annotation_col = tmp)
exf=pheatmap(scale(cor(log2(exprset+1))))
pic <- list(corf=corf,exf=exf)
}
