技巧系列
截?cái)?types<-c("protein_coding","lncRNA","processed_pseudogene","others")
p1<-ggboxplot(human_new_df, x = "type", y = "polya_length",
palette ="npg",fill = "type",legend = "none",
ylab = "Poly(A) Length")+coord_cartesian(ylim = c(0,160))+
theme(axis.text.x = element_text( vjust = 1, hjust = 1, angle = 45))
my_comparisons <- list(combn(types,2)[,1],combn(types,2)[,2],combn(types,2)[,3],combn(types,2)[,4],combn(types,2)[,5],combn(types,2)[,6] )
p2<-ggboxplot(human_new_df, x = "type", y = "polya_length",
palette ="npg",fill = "type")+
labs(x=NULL,y=NULL,fill=NULL) +
theme(axis.text.x = element_blank(),axis.ticks.x = element_blank(),axis.line.x= element_blank())+
coord_cartesian(ylim = c(600,800))+scale_y_continuous(breaks = c(600,800,100))
# p2<-p2+stat_compare_means(comparisons = my_comparisons,
# paired = F,method = "t.test",position = "identity",label = "p.signif")
ggarrange(p2,p1,heights=c(1/8, 7/8),ncol = 1, nrow = 2,common.legend = TRUE,legend="top",align = "v")
1、成對(duì)分析
常規(guī)用法
rm(list=ls())
suppressMessages(library(ggpubr))
suppressMessages(library(tidyverse))
suppressMessages(library(vcd))
my_comparisons <- list(c("control", "treatment"))
ggpaired(comb_fre_new, x = "xj", y ="percentage",ylab ="Gene 3' UTR Modification Frequence\n(% of Total Modificantion Events)",
line.color = "gray", line.size = 0.05,
palette = "npg"#c("#0A5EB9","#DF3D8C"
)
p<-p+stat_compare_means(comparisons = my_comparisons,
paired = T,method = "t.test",position = "identity",label = "p.signif")
做完是一根線，所以進(jìn)階改良（還是ggplot2 靠譜）
主要就是定義一個(gè)數(shù)xj 讓它實(shí)現(xiàn)偏移
comb_fre_new$Group <- factor(comb_fre_new$Group, levels=unique(comb_fre_new$Group))
comb_fre_new$"x"<-c(rep(1,nrow(fed_fre)),rep(2,nrow(fasting_fre)))
comb_fre_new$xj <- jitter(comb_fre_new$x, amount=.04)
p<-ggplot(data=comb_fre_new, aes(y=percentage)) +
geom_boxplot(aes(x=Group, group=Group), width=0.2, outlier.shape = NA) +
geom_point(aes(x=xj,colour = factor(Group))) + scale_colour_manual(name="",values = c("control"="#0A5EB9", "treatment"="#DF3D8C"))+
geom_line(aes(x=xj, group=gene_name),size=0.05,color="gray") +theme_bw()
2、相關(guān)性分析
suppressMessages(library(ggpubr))
sp<-ggscatter(utr3_mean, x = "control", y = "treatment",
add.params = list(color = "black",linetype="dashed"), # Customize reg. line
add = "reg.line", # Add regressin line
conf.int = F,color ="#54C6DC",
ellipse.alpha=0.5,
xlab = "control",
ylab="treatment",
label = "Gene_symbol",
label.select = targets,
repel = T,
font.label = c(11, "bold", "#f16446")
)
sp<-sp + stat_cor(method = "pearson")
#填充
sp<-ggscatter(det_df, x = "DETs_logFC", y = "m6A_logFC",fill="divergence",
palette=c("#4DBBD54C","#E64B354C"),llipse.alpha=0.3,
size=3,shape=21,color="black",
xlab="log2(Fold Change) of Transcript Expression Levels",
ylab="log2(Fold Change) of m6A Modificantion Frequence", # Customize reg. line
label="new_name",
label.select=labels,
repel=T,label.rectangle=F,
font.label=c(4,"#3C5488FF")
)
sp<-sp + stat_cor(method = "pearson")
#不同顏色的線
sp<-ggscatter(polysome_data_filter_df, x = "polysome_f_a", y = "m6a_log2FC",
fill="group",
palette=c("#4DBBD54C","#E64B354C"),
#llipse.alpha=0.3,
size=3,
shape=21,
color="black",
fullrange = TRUE,
add = "reg.line",
add.params = list(color = "group",linetype="dashed"),
xlab="log2(Fold Change) of Polysomal Transcripts",
ylab="log2(Fold Change) of m6A Modificantions", # Customize reg. line
label="new_name",
label.select=labels,
repel=T,label.rectangle=F,
font.label=c(6,"#3C5488FF")
)
sp<-sp + stat_cor(aes(color = group), label.x = 3)
sp
3、 boxplot
#例子1
new_data<-new_ccle_data %>% gather(key=Gene_name, value=expression_levels,-pick_levels)
colnames(new_data)<-c("Levels","Gene_name","Relative_expression_levels")
filename=paste0(data_dir,"/",j,"_expression_levels_in_",i,"_Median.pdf")
pdf(filename,width=3.5,height=5)
p <- ggboxplot(new_data, x = "Gene_name", y = "Relative_expression_levels",
color = "Levels", palette ="jco",
shape = "Levels",group="Levels")
p<-p + stat_compare_means(aes(group = Levels),label = "p.signif",hide.ns=F, paired=F,method = "t.test")
plot(p)
dev.off()
#例子2
p<-ggboxplot(target_df, x = "group", y = pick_target, add = c("jitter"),add.params = list(color="group"),
color = "black", palette =c("#B3B3B3","#CE86AF"),ylab="Relative Expression Levels",
shape = "group",group="group")
p<-p + stat_compare_means(aes(group = group),label = "p.signif",hide.ns=F, paired=F,method = "t.test")
ggsave(p,filename = paste0("~/Desktop/",pick_target,"_expression_TCGA.pdf"),height = 5,width = 2.5,units = "in")
4、小提琴圖
my_comparisons<-list(c("control","treatment"))
p<-ggviolin(human, x = "Group", y = "ELAVL1",
palette = c("#54C6DC","#F16446"),fill = "Group",
add = "boxplot", add.params = list(fill = "white"),
ylab = "Relative Expression Levels\n(TPM)")
p<-p + stat_compare_means(comparisons = my_comparisons,label = "p.signif",method = "t.test")
ggsave(p,filename ="~/Desktop/hur_kd/human_hur_kd_violin.pdf",width = 3,height = 5,units = "in" )
5、直方圖
#例子1
p <- ggbarplot(data, x = "New_name", y = "Relative_Exp",
color = "Treatment", palette = c("#4DBBD5","#E64B35"),
ylab = "Relative Expression Levels",position = position_dodge(0.9),
fill = "Treatment",alpha=0.2,add = c("mean_se", "jitter"),add.params = list(size=0.5))
p<-p + stat_compare_means(aes(group = Treatment),label = "p.signif",hide.ns=F, paired=F,method = "t.test")
p<-p+theme(axis.text.x = element_text(size = 7, vjust = 0.7, hjust =0.7, angle = 45))
#例子2 橫向
ggbarplot(df_top20, "pathway", "Frequence", orientation = "horiz",
fill = "Treatment", color = "Treatment", palette = c("#DF3D8C","#0A5EB9"),
label = TRUE,
position = position_dodge(0.95))+ scale_x_discrete(labels = wrap_format(30))
#例子3 橫向占比
dis_rt$distribution <- factor(dis_rt$distribution, levels=unique(dis_rt$distribution))
ggbarplot(dis_rt, "donors", "Freq",
fill = "distribution", color = "black", palette = "npg",orientation = "horiz",
label = TRUE,lab.col = "white", lab.pos = "in",ylab = "Percentage (%)")
#例子4 計(jì)算各組比例
suppressMessages(library(tidyverse))
suppressMessages(library(ggpubr))
Kid.combined@meta.data %>%
group_by(seurat_clusters,orig.ident) %>%
count() %>%
group_by(seurat_clusters) %>%
mutate(percent=100*n/sum(n)) %>%
ungroup() %>%
ggbarplot(x="seurat_clusters",y="percent", fill="orig.ident",color = "black",palette = "jco")
#例子5
data<-degs %>% left_join(pcr,by="New_name")
data<-data[complete.cases(data$Relative_Exp),]
p <- ggbarplot(data, x = "New_name", y = "Relative_Exp",
color = "Treatment", palette = c("#4DBBD5","#E64B35"),
ylab = "Relative Expression Levels",position = position_dodge(0.9),
fill = "Treatment",alpha=0.2,add = c("mean_se", "jitter"),add.params = list(size=0.5))
p<-p + stat_compare_means(aes(group = Treatment),label = "p.signif",hide.ns=F, paired=F,method = "t.test")
p<-p+theme(axis.text.x = element_text(size = 7, vjust = 0.7, hjust =0.7, angle = 45))
#例子6
#帶點(diǎn)的bar圖
p<-ggbarplot(df, x = "group", y = "exp",add = c("mean_se", "dotplot"),
color = "group", palette = c("#0A5EB9","#EBB208"),
fill = "group",alpha=0.2,
ylab = "Plasma TG Levels (mg/dL)")+coord_cartesian(ylim = c(120,165))
my_comparisons<-list(c("V5","hAS1"))
p + stat_compare_means(comparisons = my_comparisons,label = "p.signif",method = "t.test")
6、散點(diǎn)圖
#外圈為黑色的圖
#fill 為固定顏色
library(ggplot2)
ggplot(data=new_data, aes(x =groups , y = Motif)) +
geom_point(aes(size = value),shape = 21,fill="#b2eb08", colour = "black")+theme_bw()+theme_classic()
#fill為變量
plot_pathway<-function(new_data){
suppressMessages(library(ggplot2))
suppressMessages(library(scales))
suppressMessages(library(gridExtra))
suppressMessages(library(ggthemes))
suppressMessages(library(stringr))
new_data$Group <- factor(new_data$Group, levels=unique(new_data$Group))
new_data$Description<-str_to_title(as.character(new_data$Description), locale = "")
p<-ggplot(data=new_data, # you can replace the numbers to the row number of pathway of your interest
aes(x =Group , y = Description)) +
geom_point(aes(size = Count,fill = -log10(pvalue)),shape = 21, colour = "black")+scale_fill_gradient2(low = "blue", high = "red",midpoint = 1.3,limit=c(min(new_data$'-log10(pvalue)'), max(new_data$'-log10(pvalue)')))+
theme_bw()+theme_classic()+
theme(axis.title.x = element_text(size=12,face="bold",colour = "black"),
axis.text.x = element_text(size=12,face="bold",colour = "black"),
axis.title.y = element_text(size=12,colour = "black",face = "bold"),
axis.text.y= element_text(size=12,face="bold",colour = "black"))+
theme(axis.text.x = element_text( face = "bold", vjust = 1, hjust = 1, angle = 45))+scale_y_discrete(labels = wrap_format(40))
return(p)
}
#火山圖
##例子一，常規(guī)火山圖
rm(list = ls())
suppressMessages(library(ggpubr))
suppressMessages(library(tidyverse))
suppressMessages(library(ggrepel))
plot_degs<-function(degs_dir,fc,padj,top){
degs<-read.csv(degs_dir)
degs$sigORnot = as.factor(ifelse(degs$padj < padj & abs(degs$log2FoldChange) > fc,
ifelse(degs$log2FoldChange > fc ,'Increased','Decreased'),'Not Significant'))
degs_down<-degs %>% filter(sigORnot=="down")%>% arrange(log2FoldChange)
degs_down_list<-as.character(degs_down$Gene)[1:top]
degs_up<-degs %>% filter(sigORnot=="up")%>% arrange(desc(log2FoldChange))
degs_up_list<-as.character(degs_up$Gene)[1:top]
degs_not<-degs %>% filter(sigORnot=="Not Significant")
degs<-as.data.frame(rbind(degs_down,degs_not,degs_up))
degs$log_p<-(-log10(degs$padj))
leables<-c(degs_down_list,degs_up_list)
degs$sigORnot <- factor(degs$sigORnot, levels=unique(degs$sigORnot))
p<-ggscatter(degs,x="log2FoldChange",y="log_p",fill="sigORnot",
palette=c("#4DBBD5FF","grey70","#E64B35FF"),size=3,
ellipse.alpha=0.3,shape=21,color="black",
xlab="log2(FoldChange)",
ylab="-log10(Adjustedp-value)",
label="Gene",
label.select=leables,
repel=T,label.rectangle=T,
font.label=c(8,"bold","#3C5488FF")
)
p=p+geom_hline(yintercept=(-log10(padj)),linetype=3)+geom_vline(xintercept=c(-(fc),fc),linetype=3)
return(p)
}
##例子2 突出重點(diǎn)的點(diǎn)
###構(gòu)建第二個(gè)矩陣，用來(lái)重點(diǎn)標(biāo)記
rm(list = ls())
suppressMessages(library(ggplot2))
suppressMessages(library(ggrepel))
pvalue=0.05
fc_value=0.5
degs<-read.csv(deg_dir)
out_tab<-read.csv(paste0(out_dir,"conserved_PRJNA523510_NALF_NASH_overlap.csv"))
out_tab<-out_tab %>% dplyr::select("human_gene_id","new_name")
leables<-as.character(out_tab$new_name)
deg2<-degs %>% dplyr::rename(human_gene_id=X) %>% left_join(out_tab,by="human_gene_id")
deg2[is.na(deg2)]<-""
deg2$sigORnot<-ifelse(deg2$new_name=="","NA","FCLs")
deg2$sigORnot <- factor(deg2$sigORnot, levels=unique(deg2$sigORnot))
deg2$log_p<-(-log10(deg2$pvalue))
deg2$fc<-deg2$log2FoldChang
#選取部分作為標(biāo)注矩陣deg3
deg3<-deg2 %>% filter(sigORnot!="NA")
ggplot(deg2,aes(x=fc,y=log_p))+geom_point(shape=21,color="#CCCCCC")+
#coord_equal()+
geom_point(data=deg3,aes(x=fc,y=log_p),colour="#800080")+
geom_text_repel(data=deg3,aes(x=fc,y=log_p),label=leables,size=2)+
theme_bw()+theme_classic()+
xlim(-7.5,7.5)+ylim(0,38)+
xlab("log2(FoldChange)")+ylab("-log10(pvalue)")+
geom_hline(yintercept=(-log10(pvalue)),linetype=3)+geom_vline(xintercept=c(-(fc_value),fc_value),linetype=3)
#例子3 普通點(diǎn)圖
ggplot(data=cv_data, # you can replace the numbers to the row number of pathway of your interest
aes(x =rank , y = cv)) +
geom_point(aes(fill = group,alpha=0.7),shape = 21,colour = "black" )+scale_fill_manual(values =c("#DF3D8C","#0A5EB9"))+
theme_bw()+theme_classic()
#點(diǎn)圖 label
ggplot(data=df_pick,
aes(x =log2(tpm) , y = log2(n_c_ratio),label = treatment))+
geom_point(aes(fill = gene_name),shape = 21, colour = "black",size=2)+
scale_fill_manual(values =c("#F16446","#54C6DC"))+theme_bw()+theme_classic()+
xlim(c(-2,12))+ylim(c(-4,5.5))+
xlab("Log2(Whole Cell TPM)")+ylab("Log2(nucl. TPM/cyto. TPM)")+
geom_text(aes(colour = factor(gene_name)),hjust = 0, nudge_x = 0.2,size=3)
重點(diǎn)標(biāo)記火山圖
7、分布圖
ggdensity(rt, x = "polya_length",
palette = c("#54C6DC","#F16446"), add = "mean", color = "doners",
ggtheme = theme_light(), legend = "top",xlab = "Poly(A) length")
8、畫(huà)餅圖
rm(list = ls())
suppressMessages(library(tidyverse))
suppressMessages(library(ggpubr))
suppressMessages(library(vcd))
distribution<-function(x){
if (x=="protein_coding"){
y="protein_coding"
}else if(x=="lncRNA"){
y="lncRNA"
}else if(x=="processed_pseudogene"){
y="processed_pseudogene"
} else {
y="others"
}
return(y)
}
data_1st$type<-unlist(lapply(as.character(data_1st$gene_type), FUN = distribution))
signature<-function(rt){
signature<-with(rt,table(type))
signature<-as.data.frame(prop.table(signature)*100)
signature$labs<-paste0(as.character(round(as.numeric(signature$Freq),2)),"%")
signature_rt=data.frame()
for (i in c("lncRNA","protein_coding","processed_pseudogene","others")){
signature_rrt<-signature[signature$type==i,]
signature_rt<-rbind(signature_rt,signature_rrt)
}
return(signature_rt)
}
data_1st_type<-signature(data_1st)
plot<-function(signature){signature$type <- factor(signature$type, levels=unique(signature$type))
p<-ggpie(signature,"Freq",
label = "labs",
lab.pos = "out", lab.font = "white",
fill = "type",
color = "black",
palette = "npg")
return(p)
}
data_1st_plot<-plot(data_1st_type)
data_1st_plot
折線圖
ggline(human_tpm_treatment_clean, "human_gene_id", "cv",
linetype = "treatment", shape ="treatment",ylab = "Coefficient of Variation",
color = "treatment", palette ="npg",size=0.1,plot_type="b")+scale_shape_manual(values=c(19,1,2,3))+
theme(
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
點(diǎn)線圖
ggdotchart(human_mouse_mean_rt, x = "human_name", y = "logFC",
color = "group", # Color by groups
palette = c("#4DBBD5FF","#E64B35FF"), # Custom color palette
sorting = "descending", # Sort value in descending order
add = "segments", # Add segments from y = 0 to dots
add.params = list(color = "lightgray", size = 0.2), # Change segment color and size
group = "group", # Order by groups
dot.size = 2.5,
shape = "group",ylab = "log2FoldChange",
ggtheme = theme_pubr(),
label="human_name",
label.select=human_mouse_combine_label,
label.rectangle=F
# ggplot2 theme
)+geom_hline(yintercept=0,linetype=2)
#折線分位圖
gene_dist <- dist(df_zscore)
gene_hclust <- hclust(gene_dist, method = "complete")
plot(gene_hclust, labels = FALSE)
gene_cluster <- cutree(gene_hclust, k = 5) %>%
# turn the named vector into a tibble
enframe() %>%
# rename some of the columns
rename(gene = name, cluster = value)
df_zscore$gene<-row.names(df_zscore)
df_zscore_spread<-df_zscore %>% gather(key = "sample",value="frequence",-gene)
gene_cluster_zscore<-gene_cluster %>% left_join(df_zscore_spread,by="gene") %>%
mutate(group=unlist(lapply(as.character(sample), FUN = function(x) {return(strsplit(x, split = "_",fixed = T)[[1]][1])}))) %>%
mutate(treatment=unlist(lapply(as.character(sample), FUN = function(x) {return(strsplit(x, split = "_",fixed = T)[[1]][2])})))
gene_cluster_zscore$"group"<-factor(gene_cluster_zscore$group, levels=unique(gene_cluster_zscore$group))
gene_cluster_zscore$"treatment"<-factor(gene_cluster_zscore$treatment, levels=unique(gene_cluster_zscore$treatment))
gene_cluster_zscore$"cluster"<-factor(gene_cluster_zscore$cluster, levels=unique(gene_cluster_zscore$cluster))
gene_cluster_zscore %>%
ggplot(aes(group, frequence)) +
geom_line(aes(group = gene,colour = cluster),size=0.1) +
facet_grid(rows = vars(treatment), cols = vars(cluster)) +
scale_colour_manual(values =c("#E64B35FF", "#4DBBD5FF","#00A087FF","#3C5488FF","#F39B7FFF"))+
theme_bw()
折線分位圖
9、相關(guān)性分析
10、韋恩圖
11、個(gè)性化繪制Pathway
12 、Heatmap
13、sample distance heatmap
suppressMessages(library(DESeq2))
suppressMessages(library("gplots"))
sampleDists <- dist( t(tpm_zscore_tissues) )
sampleDistMatrix <- as.matrix( sampleDists )
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
hc <- hclust(sampleDists)
heatmap.2( sampleDistMatrix, Rowv=as.dendrogram(hc),
symm=TRUE, trace="none", col=colors,
margins=c(2,10), labCol=FALSE )
14、雷達(dá)密度圖
ggplot(df, aes(x = log2(all_tpm_mean), y = log2(total_n_c))) +
stat_density2d(geom="density2d", aes(color = type,alpha = ..level..),contour=T,position = "identity")+
scale_color_manual(values =c("#54C6DC","#F16446"))+theme_bw()+theme_classic()+
geom_hline(yintercept=0,linetype=3)+xlim(c(-2,12))+ylim(c(-4,5.5))+
xlab("Log2(Whole Cell TPM)")+ylab("Log2(nucl. TPM/cyto. TPM)")
ggplot(df, aes(x=log2(al_n_c), y=log2(fast_n_c))) +
stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
scale_fill_distiller(palette=4, direction=-1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme(
legend.position='none'
)
ggplot(df, aes(x=log2(al_n_c), y=log2(fast_n_c)) ) +
geom_hex(bins = 60) +
scale_fill_continuous(type = "viridis") +
theme_bw()
ggplot(df, aes(x=log2(al_n_c), y=log2(fast_n_c)) ) +
stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white")
15、GSEA自定義做圖
16、qPCR的擴(kuò)增曲線
#用excel自己準(zhǔn)備每個(gè)孔，每個(gè)cycle對(duì)應(yīng)的熒光值、Rn值或者delta Rn值，對(duì)應(yīng)函數(shù)里面的名字也要改
#字體之類(lèi)的自行調(diào)整
rm(list = ls())
suppressMessages(library(ggplot2))
suppressMessages(library(tidyverse))
plot_amp<-function(df){
df$log_d_rn<-(log10(df$Delta.Rn))
p<-ggplot(df, aes(x = Cycle, y = log_d_rn,color = Sample,linetype = Group)) +
stat_smooth(aes(x = Cycle, y = log_d_rn),
method = "loess", span = 0.6, se = FALSE, fullrange = F) +
scale_x_continuous(limits = c(1, 40)) +
labs(x = "Cycles", y = "log 10 (Delta Rn)", title = "qPCR Amplification Curve")+
scale_linetype_manual(values=c("twodash", "dotted","solid"))+
theme_bw()+theme_classic()
return(p)
}
fcl1<-read.csv("~/Desktop/amp_hfcl1.csv")
fcl1<-fcl1 %>% filter(Target.Name=="XX-Q1")
fcl1_plot<-plot_amp(fcl1)
#調(diào)整下輸出圖像的范圍
fcl1_plot_adjust<-fcl1_plot +coord_cartesian(ylim = c(-2.8,0.5),xlim=c(15,40))
例子
