Notebook

Tissue Specific Motif Networks¶

In [38]:

net <- read.table('/public/workspace202011/encode/zhutao/analysis/gimmemotif/top2_networks.csv',header=T)
net <-  net[!duplicated(net[c(1,3)]),]
net$type <- net$Tissue
net <- net[,c("Tissue","Factor","zscore","type")]
head(net)

A data.frame: 6 × 4
	Tissue	Factor	zscore	type
	<chr>	<chr>	<dbl>	<chr>
1	AM1	LOC_Os02g01380	4.906071	AM1
2	AM1	OsLFL1	4.277486	AM1
3	AM1	OsERF117	3.451406	AM1
4	AM1	EIL1	3.415054	AM1
5	AM1	LOC_Os07g07974	3.330266	AM1
6	AM2	OsLFL1	4.454432	AM2

In [41]:

convert_ID <- read.table('/public/workspace202011/encode/zhutao/result/TF_MSU_RAP_converID.csv',head=T)
aa <- merge(net, convert_ID, by.x="Factor",by.y="gene",all=F)
net <- net[net$type!="AM2" & net$type!="Panicle2" &net$type!="Panicle3" & net$type!="Panicle4" & net$type!="SAM2" & net$type!="Seed1" & net$type!="Seed3" ,]
anno = data.frame(a=unique(net$Tissue),b='Tissues')
tf = data.frame(a=unique(net$Factor),b='TF')
c = rbind(anno,tf)
head(c)

A data.frame: 6 × 2
	a	b
	<chr>	<chr>
1	AM1	Tissues
2	Callus	Tissues
3	DBuds	Tissues
4	Leaf	Tissues
5	Lemma	Tissues
6	Palea	Tissues

In [43]:

library(igraph)
library(paletteer)
library(qgraph)
net3 <- graph_from_data_frame(d=net, vertices=c, directed=F) 
set.seed(2021)
mycol <- paletteer_d("ggthemes::Classic_20")[1:17]
V(net3)$size <- ifelse(V(net3)$b=="Tissues",8, 25)
V(net3)$shape <- ifelse(V(net3)$b=="Tissues", "circle","rectangle")
V(net3)$frame.color <- ifelse(V(net3)$b=="Tissues", NA,"black")
E(net3)$color <- 'grey'
options(repr.plot.width = 25, repr.plot.height = 25, repr.plot.res = 100)
e <- get.edgelist(net3,names=FALSE)
#e <- get.edgelist(net3,names=FALSE)
l <- qgraph.layout.fruchtermanreingold(e,vcount=vcount(net3),area=8*(vcount(net3)^2),repulse.rad=(vcount(net3)^3.1))
#l <-layout_with_dh(net3) 
#pdf("/public/workspace202011/encode/zhutao/figure/tf_net_works_top5.pdf",width=16,height=16)
plot(net3,vertex.color=mycol,edge.width=(E(net3)$zscore-1),edge.curved = F,layout = l,vertex.label.cex=1.5)
#dev.off()

Young Panicle Dynamic OCRs¶

In [ ]:

quant <- read.table("/public/workspace202011/encode/zhutao/analysis/Figure/Fig4/panicle_dynamic/quant/peak_quant_norm.csv",head=T,row.names=1)
library(philentropy)
entropy <- apply(quant, 1, function(x){H(x/sum(x))})
aa <- as.data.frame(H_specificity)
head(data.frame(aa[order(H_specificity),]),5001)
specificty_score <- data.frame(score=H_specificity, peak = names(H_specificity))
spec_peak <- specificty_score[specificty_score$score <=0.0008074347,]
quant[rownames(quant) %in% spec_peak$peak,]
write.table(spec_peak, "/public/workspace202011/encode/zhutao/analysis/Figure/Fig4/panicle_dynamic/cluster1_6_motif/5000_Non_specific_peak.txt",sep='\t',quote=F,col.names=F,row.names=F)

In [1]:

res <- read.table("/public/workspace202011/encode/zhutao/analysis/Figure/Fig4/panicle_dynamic/clutser1_6_merged_peak_TF_CHIP_enrichment.txt")
colnames(res) <- c("TF","left","right","two_tail","ratio")
res$two_tail <- ifelse(res$two_tail==0, 10e-300,res$two_tail)
res$score <- -log10(res$two_tail)
res$name <- substring(res$TF,11)
head(res)

A data.frame: 6 × 7
	TF	left	right	two_tail	ratio	score	name
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>
1	DRP000207_OSH1	0.00000	1.0000e+00	1.0000e-299	7.618	299.0000000	OSH1
2	DRP001345_OsTGAP1	1.00000	3.1699e-100	3.1699e-100	2.934	99.4989544	OsTGAP1
3	ERP106693_OsIDS1	0.90409	1.1514e-01	2.0542e-01	1.139	0.6873573	OsIDS1
4	ERP106918_NGR2	0.18932	9.3243e-01	3.3582e-01	0.455	0.4738934	NGR2
5	ERP108685_bZIP73	1.00000	4.4624e-10	4.4624e-10	3.720	9.3504315	bZIP73
6	SRP014840_OsSRT1	1.00000	1.6307e-08	1.7340e-08	2.518	7.7609509	OsSRT1

In [2]:

library(ggplot2)
library(ggpubr)
library(paletteer)
library(ggrepel)
options(repr.plot.width = 8, repr.plot.height = 5, repr.plot.res = 100)
p1 <- ggplot(res, aes(x=ratio, y=score,label=name))+theme_pubr()+geom_point(aes(size= ratio, fill = score),shape=21,alpha=0.9)+
scale_fill_gradientn(colours = paletteer_d("ggsci::pink_material")[1:6], limits=c(0, 300))+theme_pubr(base_size = 12, base_family = "", border = TRUE, margin = TRUE, 
                      legend = c("right"), x.text.angle = 0)+scale_size_continuous(range = c(0, 8))+
ylab('-log10(P-value)') +xlab("Fold enrichment")+geom_text_repel(data=subset(res, score > 2 & ratio > 2),
            aes(ratio,score,label=name))+ggtitle("56 ChIP-seq Enrchment")
p1

In [3]:

library(stringr)
sea <- read.table("/public/workspace202011/encode/zhutao/analysis/Figure/Fig4/panicle_dynamic/cluster1_6_motif/merged_res/sea.tsv",head=T)
sea <- sea[,c(4,10,15)]
sea$LOG_EVALUE <- -sea$LOG_EVALUE
colnames(sea) <- c("name","ratio","score")
sea$motif <- str_sub(sea$name,1,16)
sea$motif <- gsub("[()]","",sea$motif) 
sea <- sea[sea$score>0,]
# convert_ID <- read.table('/public/workspace202011/encode/zhutao/result/TF_MSU_RAP_converID.csv', head=T)
# merge(sea,convert_ID, by.x="motif",by.y="MSU", all=F)
library(ggplot2)
library(ggpubr)
library(paletteer)
library(ggrepel)
options(repr.plot.width = 8, repr.plot.height = 6, repr.plot.res = 100)
p2 <- ggplot(sea, aes(x=ratio, y=score,label=motif))+theme_pubr()+geom_point(aes(size= ratio, fill = score),shape=21,alpha=0.9)+
scale_fill_gradientn(colours = paletteer_d("ggsci::teal_material")[1:6], limits=c(0, 300))+theme_pubr(base_size = 12, base_family = "", border = TRUE, margin = TRUE, 
                      legend = c("right"), x.text.angle = 0)+scale_size_continuous(range = c(0, 8))+
ylab('-log10(P-value)') +xlab("Fold enrichment")+geom_text_repel(data=subset(sea, score > 2 & ratio > 1.5),
            aes(ratio,score,label=motif))+ggtitle("Motif Enrichment")
p2

Tissue Specific Motif Family Analysis¶

In [5]:

setwd("/public/workspace202011/encode/zhutao/analysis/tf_enrichment/top2500_peaks/sea_out_new/")
file_list = c(
    "AM1" = "AM/sea.tsv",
    "Callus" = "NIP_Callus/sea.tsv",
    "Dbuds" = "NIP_Dormant_buds_umi/sea.tsv",
    "Leaf" = "NIP_Leaf_4_5/sea.tsv",
    "Sheath" = "NIP_Leaf_sheath_4_5/sea.tsv",
    "Lemma"= "NIP_Lemma_umi/sea.tsv",
    "Palea" ="NIP_Palea_umi/sea.tsv",
    "PNN" = "NIP_Panicle_neck_node_umi/sea.tsv",
    "Pistil"= "NIP_Pistil_umi/sea.tsv",
    "Plumule" = "NIP_Plumule/sea.tsv",
    "Radicle" = "NIP_Radicle_umi/sea.tsv",
    "Root" = "NIP_Root_4_5/sea.tsv",
    "Stamen"="NIP_Stamen_umi/sea.tsv",
    "Stem" ="NIP_Stem_umi/sea.tsv",
    "Panicle1"="Panicle/sea.tsv",
    "SAM"="SAM/sea.tsv",
    "Seed"="Seed/sea.tsv"
)
library(GenomicRanges)
peak_list = lapply(file_list, function(f) {
    df = read.table(f,head=T)
    df <- df[df$QVALUE < 1e-10,]
    a = unique(df$ALT_ID)
})

#lengths(peak_list)

library(UpSetR)
library(ggpubr)
# options(repr.plot.width = 16, repr.plot.height = 8, repr.plot.res = 100)
# upset(fromList(peak_list), order.by = "freq",sets = names(lengths(peak_list)),point.size = 3, line.size = 1, 
#     mainbar.y.label = "Intersectionh size", sets.x.label = "Target genes")
#`````````````````````````````````````````My Custom Color Upset Plot`````````````````````````


tissue_col <- c(Leaf="#2ca02c",Sheath="#98df8a",Stem="#32a251",PNN="#bdbdbd",Lemma= "#ff9896",Palea="#9467bd",Pistil="#c5b0d5",
               Stamen="#8c564b",Seed="#e377c2",SAM1="#ff7f0e",AM1="#ff7f0f",
               DBuds="#b85a0d",Panicle1="#1f77b4",Plumule="#86b4a9",Radicle="#82853b",
               Callus="#ccc94d",Root ="#7f7f7f")

options(repr.plot.width = 12, repr.plot.height = 12, repr.plot.res = 100)
p3 <- myupset(fromList(peak_list), order.by = "freq",sets = names(lengths(peak_list)),point.size = 3, line.size = 1, 
    mainbar.y.label = "Intersectionh size", sets.x.label = "Enrichment motifs",sets.bar.color=tissue_col,
    main.bar.color= "#1170AA",matrix.color ="#FC7D0B", number.angles = 0)
p3

In [6]:

setwd("/public/workspace202011/encode/zhutao/analysis/tf_enrichment/top2500_peaks/sea_out_new/")
file_list = c(
    "AM1" = "AM/sea.tsv",
    "Callus" = "NIP_Callus/sea.tsv",
    "Dbuds" = "NIP_Dormant_buds_umi/sea.tsv",
    "Leaf" = "NIP_Leaf_4_5/sea.tsv",
    "Sheath" = "NIP_Leaf_sheath_4_5/sea.tsv",
    "Lemma"= "NIP_Lemma_umi/sea.tsv",
    "Palea" ="NIP_Palea_umi/sea.tsv",
    "PNN" = "NIP_Panicle_neck_node_umi/sea.tsv",
    "Pistil"= "NIP_Pistil_umi/sea.tsv",
    "Plumule" = "NIP_Plumule/sea.tsv",
    "Radicle" = "NIP_Radicle_umi/sea.tsv",
    "Root" = "NIP_Root_4_5/sea.tsv",
    "Stamen"="NIP_Stamen_umi/sea.tsv",
    "Stem" ="NIP_Stem_umi/sea.tsv",
    "Panicle1"="Panicle/sea.tsv",
    "SAM"="SAM/sea.tsv",
    "Seed"="Seed/sea.tsv"
)
library(GenomicRanges)
peak_list = lapply(file_list, function(f) {
    df = read.table(f,head=T)
    df <- df[1:100,]#df[df$QVALUE < 1e-25,]
    a = unique(df$ALT_ID)
})
motif <- stack(peak_list)
motif$motif <- sapply(strsplit(motif$values,")"), `[`, 1)
motif$moitf <- gsub("[()]","",motif$motif)
head(motif)

A data.frame: 6 × 4
	values	ind	motif	moitf
	<chr>	<fct>	<chr>	<chr>
1	(LOC_Os01g68370)_(Arabidopsis_thaliana)_(DBD_0.87)	AM1	(LOC_Os01g68370	LOC_Os01g68370
2	LOC_Os01g68370	AM1	LOC_Os01g68370	LOC_Os01g68370
3	LOC_Os01g51610	AM1	LOC_Os01g51610	LOC_Os01g51610
4	LOC_Os07g39320	AM1	LOC_Os07g39320	LOC_Os07g39320
5	(LOC_Os03g20780)_(Arabidopsis_thaliana)_(DBD_0.60)	AM1	(LOC_Os03g20780	LOC_Os03g20780
6	LOC_Os02g01380	AM1	LOC_Os02g01380	LOC_Os02g01380

In [8]:

fam <- read.table("/public/workspace202011/encode/zhutao/Osj_TF_list.txt",head=T)
sanky <- merge(motif, fam, by.x="motif",by.y="Gene_ID",all=F)
sanky$fre <- 1
sanky <- sanky[,c(3,6,7)]
#sanky$Family <- gsub("BBR-BPC|BES1|CAMTA|E2F/DP|GeBP|SRS|ZF-HD|C3H|EIL|RAV|LFY|LBD|AP2|B3|GATA|WOX|CPP|TALE|ARR-B|Trihelix","Others",sanky$Family)
#sanky$Family <- gsub("M-type_MADS|MIKC_MADS","MADS",sanky$Family)
sanky$tissue <- as.character(sanky$ind)
library(data.table)
library(tidyverse)
res <- table(sanky$tissue,sanky$Family) %>% as.data.frame()%>% reshape(idvar = "Var2", timevar = "Var1", direction = "wide")%>% 
       remove_rownames %>% column_to_rownames(var="Var2")
colnames(res) <- gsub("Freq.","",colnames(res))

In [11]:

library(pheatmap)
p43 <- pheatmap(res/do.call(pmax, res), show_rownames=F)
library(magrittr)
library(ComplexHeatmap)
library(pheatmap)
p3 = draw(p43)
r.dend <- row_dend(p3)

Warning message:
“The input is a data frame, convert it to the matrix.”

In [14]:

rcl.list <- row_order(p3)
#rownames(res[rcl.list,])
res2 <- table(sanky$tissue,sanky$Family) %>% as.data.frame()
res2 <- res/do.call(pmax, res)
res2$fam <- rownames(res2)
res3 <- melt(res2)
res3$fam <- factor(res3$fam, levels=rev(rownames(res[rcl.list,])))
res3$variable <- factor(res3$variable, levels=c("Lemma","Palea","AM1","Seed","Radicle","Root","Pistil","Dbuds","SAM","Stamen",
                                                "PNN","Stem","Leaf","Panicle1","Sheath","Callus","Plumule"))
head(res3)

Warning message in melt(res2):
“The melt generic in data.table has been passed a data.frame and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(res2). In the next version, this warning will become an error.”
Using fam as id variables

A data.frame: 6 × 3
	fam	variable	value
	<fct>	<fct>	<dbl>
1	AP2	AM1	1
2	ARF	AM1	0
3	ARR-B	AM1	0
4	B3	AM1	1
5	BES1	AM1	0
6	bHLH	AM1	0

In [15]:

myPalette <- colorRampPalette(as.character(paletteer_d("RColorBrewer::OrRd")[2:6]))
options(repr.plot.width = 8, repr.plot.height = 10, repr.plot.res = 100)
p4 <- ggplot(res3, aes(x=variable, y=fam)) + geom_point(aes(size= value,color=value),shape=15,alpha=0.9)+
scale_color_gradientn(colours =rev(paletteer_d("RColorBrewer::Spectral")[2:10]))+theme_bw()+theme(axis.text.x = element_text(angle = 270, vjust = 1, hjust=1))+scale_size_continuous(range = c(0, 7))+
ylab('') +xlab("")
p4

Tissue Specific Motif Analysis¶

In [35]:

test <- read.table('/public/workspace202011/encode/zhutao/ATACorrect/maelstrom/output_new/jaspar_aggregated_zscore_threshold2_cisbp.csv',sep='\t',head=T)
rownames(test) <- make.names(test[,1], unique = TRUE)
test <- subset(test,select=-c(factors))
test <- test[rownames(test) !="X",]
colnames(test) <- sub("z.score.","",colnames(test))
test

A data.frame: 96 × 24
	AM1	AM2	Callus	DBuds	Leaf	Lemma	Palea	Panicle1	Panicle2	Panicle3	⋯	Radicle	Root	SAM1	SAM2	Seed1	Seed2	Seed3	Sheath	Stamen	Stem
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
LOC_OS06G49830	-0.4666709	-0.32219985	4.17406645	-2.49179593	4.10970587	-0.603527936	-0.4312044	1.733376793	2.3203368	1.6043384	⋯	-4.37523237	-2.67189131	-3.16711542	-2.86476495	-0.42654862	-2.862051027	-1.56968550	2.93446106	3.61664281	3.66614293
HSFB2C.HSFB4B.OSHSFA6A.OSHSFC1B	-1.6975657	-0.78494324	0.41648706	-1.36736578	1.73501270	1.545557797	2.2428570	-0.217915047	0.6795691	-0.1775205	⋯	-1.40586964	-0.31356343	-2.10352840	-2.20378887	1.18905282	1.276868717	1.27585308	0.06600621	1.55461464	0.89141070
ONAC075	-1.9485590	-1.31263865	-0.04508732	-1.76837722	2.62086614	0.425970326	-0.4815454	1.959411072	2.1116965	1.8846729	⋯	-2.35857934	-1.05539112	-2.28674351	-2.21126628	-1.45764823	-2.055007637	2.82449225	2.78766586	-0.20824467	-0.67836387
JMJ702	0.7464345	-0.42476239	-0.06052989	1.99687707	-1.37598287	-1.758379793	-2.3085303	0.647483448	-0.2602739	0.9708447	⋯	-0.39206728	-1.08470494	1.62893684	1.70256419	1.39936631	0.797305110	-1.56726482	-1.68234407	-1.92095555	-1.29703262
LOC_OS05G03020	0.7483268	1.08851936	1.78285642	-0.07019153	0.38238018	-2.415411551	-1.9127050	1.238260664	1.9143577	1.7336910	⋯	-1.56949933	-0.61260507	0.20104603	0.62338756	0.51794631	0.600054522	0.02217347	-0.19305973	-1.53122402	-0.48429320
OSSTA129	2.3582932	1.56869258	-1.29558168	2.78037733	-3.20472685	-0.727564740	-2.4286669	0.713779172	-0.8492254	0.9548959	⋯	2.36398996	1.44319500	3.31756510	3.35983002	-1.11809977	-0.931829859	-2.49288110	-2.67469240	-2.75917541	-2.56371035
LOC_OS04G32590	-1.1249430	-2.03500098	3.55145143	-1.12602926	1.97510365	-3.799647939	-3.2538903	2.954976043	3.3391403	3.0690486	⋯	-2.82542127	-0.17006796	-0.25123694	0.42259885	-0.68337519	-1.428503130	0.16850372	1.70465150	-1.74924262	-1.38677549
OSNTL5	-1.8640346	0.23589583	-1.09305872	0.09043577	-1.15356276	0.070574031	-0.2075436	-1.513828174	-1.1002998	0.3939523	⋯	-0.01576561	-1.12404087	-0.04587536	0.25864348	2.43128986	2.227916549	2.78583162	-0.70585195	0.45052086	0.10134456
OSEIL3.LOC_OS02G36510	-0.9639781	-0.57176677	0.28568334	-0.66082366	2.49510717	-0.112963252	0.2577214	2.218951194	2.6059637	2.3366873	⋯	-3.03933781	-2.07658528	-1.90833175	-1.49599330	-1.01370814	-1.330259088	1.54379426	2.52002731	1.01676606	0.77486085
RF2B.OSBZIP84.RF2A	1.3905016	1.07847919	-1.55260507	1.12785709	-1.43094003	0.658336094	0.3582394	-0.859022330	-0.9765860	-1.2188317	⋯	1.03077945	-0.72342894	1.02588377	0.42064693	0.82847023	-0.091394969	-0.60743898	-0.77495321	0.88194864	0.82190362
AP2.EREBP127	-1.6415369	-2.18273340	-0.81518665	-0.76307332	1.20316351	-1.217828572	-1.6896608	1.310356751	1.7108577	1.7636549	⋯	-0.21589939	1.22036685	-0.67897140	-0.01158834	-0.56553122	-1.231709263	1.60857264	1.51179199	-1.55618018	-2.01780016
LOC_OS03G56090	-1.2780617	-1.75450438	1.06239763	-0.85840684	2.10580939	-3.391789245	-3.3130619	1.726341997	2.3250966	2.2123740	⋯	-0.66462898	-0.58226263	-1.08722255	-0.84095871	0.02596418	-0.274380132	2.00032354	1.11911087	-0.47588857	-0.47983846
LOC_OS02G09480.LOC_OS09G24800	0.4586454	-1.00875647	0.46823122	-0.16267339	1.00944515	0.009744425	0.4692638	1.933275938	1.3244135	1.2845366	⋯	-1.51096100	-0.25688999	-0.41284377	-0.28317151	-3.10412146	-2.310826219	0.37509773	1.61375974	1.00642083	0.04358908
OSZFP	2.4662854	0.79924007	-0.16718311	3.06712059	-2.01215348	-2.447444128	-2.8022736	3.561131140	3.1241544	4.3752324	⋯	-2.11642688	-1.21814982	3.19827190	3.35319931	-3.18964537	-2.844969215	-1.13003656	-1.56272425	-2.66998286	-2.53752342
OSMYB3R.2.CSA.OSMYB52.LOC_OS01G12860	0.4603955	-0.65036125	0.12285236	0.06469490	1.05567800	-0.275757068	0.2882994	1.838096706	1.2920673	0.9017782	⋯	-0.95363094	-1.07276032	-0.38178767	-0.27752293	0.12323221	-0.403783621	-0.52131560	0.45876819	0.17014591	0.09273642
LOC_OS05G27980	-1.3914821	-1.55168775	1.89766336	-3.59990905	4.16110750	1.094513176	1.8621461	0.290001938	0.5056365	-0.7355728	⋯	-0.83217608	0.01155583	-4.63137781	-3.32828732	0.16183183	-1.448298925	1.38881496	3.84417938	1.24433784	1.06478194
LOC_OS03G25430	-4.4544317	-4.10861556	2.45831618	-3.26205344	4.73829111	-1.433563785	-1.1776592	3.627623039	3.5363532	3.0341850	⋯	-1.29172176	0.82119278	-3.20922477	-3.81638170	-2.17545167	-4.105460643	2.57809118	5.08301691	-1.20301024	-0.30377132
AP2.EREBP033.RSR1	-0.3139671	-0.95656632	1.79414849	0.81366996	1.86726115	-1.390406867	-1.0154210	3.421535879	3.2854432	3.2538820	⋯	-2.05796786	-0.11144806	0.08873796	-0.52514659	-3.45746814	-3.188878618	0.92155811	2.10752553	-0.30115217	-0.26058734
OSBHLH033	-1.0436265	-0.94690779	1.06128218	-0.19037933	1.29618531	-1.000437359	-0.5789327	0.191894439	1.1230644	0.7385293	⋯	-0.52652937	1.88407252	-1.46562807	-0.82088092	0.64639990	0.461641461	0.75314280	1.16968572	-1.01241885	0.58964816
AP2.EREBP129.AP2.EREBP127.AP2.EREBP096	3.0785596	2.37226135	-4.90607078	2.84875403	-3.38379898	1.121760247	0.9996293	-1.599008797	-1.3780338	-0.2829992	⋯	0.14691951	-1.99470058	2.67526582	2.88170019	1.31503926	2.423361751	-1.72490266	-3.42408519	0.94042886	0.23218876
OSFBH1	-3.2791869	-3.38828367	3.61239758	-4.02941626	3.21151482	-2.163296656	-2.4299452	-2.327670702	-2.6716399	-2.3404084	⋯	2.76413341	4.03967304	-4.45443168	-4.08691312	0.55626121	0.128614376	0.80429876	3.92864970	0.91547947	0.26451476
OSWRKY63	-0.7742476	-0.47057619	0.37247377	-2.06709079	0.64076463	3.073008846	3.1582465	-0.001065989	0.4407639	-0.6988458	⋯	0.89090195	2.60851608	-1.87911027	-1.01919738	-2.51629179	-3.143907931	-0.50445567	2.37457209	-1.13797731	0.02189033
OSERF117	3.4514059	2.26022600	0.28352795	3.36223058	-2.85678820	-2.440241109	-2.4761480	1.914413324	1.2847878	1.5239338	⋯	0.14722684	-0.22915662	2.91622218	2.62645416	-0.70640142	-0.661415745	-2.94199578	-2.95322234	-1.68698662	-2.26900272
OSBZIP40.OSBZIP72.LOC_OS02G10860.BZIP23	-3.4259069	-2.85978946	2.50531738	-0.68167736	0.25477552	-0.852949166	-1.9607446	1.984235267	2.4086726	2.2875339	⋯	1.13340964	1.83094666	-0.61244821	-0.38252658	-0.18220713	0.583648960	1.12232093	1.06882922	-3.23290058	-3.30671892
LOC_OS01G51154.LOC_OS01G40670	0.7465798	-0.07306627	1.46704364	-0.78269911	0.39098849	-0.184305075	0.2557395	-0.093957545	0.5456451	0.3986240	⋯	-1.27451435	-1.24342787	-0.82082125	-0.33050103	-0.01894375	0.009306786	-1.09743302	0.48092263	3.56113114	3.22386930
OSIDD1.OSIDD3.OSIDD2.OSIDD8	-2.4889667	-3.19219233	0.65640177	-0.86448858	-0.04547556	-0.675587472	-1.1934329	0.940206839	1.0786540	1.2262931	⋯	-0.59046264	1.02863424	-0.70611365	-0.26352571	4.73829111	4.459352914	1.19740414	-0.03331836	-3.40631042	-1.90483743
OSHB4	0.1187012	0.99862340	-0.72317511	1.26171674	-0.18914960	-1.096726287	-0.8484058	0.339030973	-0.4199199	0.6433346	⋯	-2.03474933	-1.71343200	1.27931413	1.54690069	-0.45621812	0.912772659	0.49544472	0.26160097	0.70700644	0.11126771
LOC_OS02G03690.OSBHLH113.BZR1.LOC_OS06G35900.OSBHLH009...	-0.6316541	-0.23699156	-0.29093388	0.73726024	0.78261983	-1.471374970	-0.9134907	0.652070747	0.5612703	0.8106614	⋯	0.30188019	-0.79177301	0.76717711	0.52696583	1.68947275	1.484037672	0.27890084	-0.02458057	-0.11262519	-1.53005569
OSMYB60.LOC_OS08G33940	-0.4374302	2.04197599	-2.12381599	-0.22352102	-2.06135156	0.428539612	-0.7407595	-2.845150705	-2.6125566	-2.6858368	⋯	2.51413378	0.26899765	-0.87707972	-0.80811792	3.31574486	4.356684849	2.44717226	-1.90935784	0.98517926	1.10144403
ONAC052	-1.9657302	-1.97174601	2.44847368	-0.95195186	0.67565802	-1.055366471	-1.2490403	1.935927847	2.2720826	1.9596132	⋯	-2.05810907	0.72779617	-1.46719122	-1.53049671	-0.88147626	-0.793482856	2.26221126	0.37202383	0.07781289	-0.98837106
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋱	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
ONAC015.NAC2.ONAC007.ONAC020.ONAC006...	0.93032840	2.343887524	0.31707686	0.3775908	-2.6388757	-0.081231512	-0.58462737	-1.7449803702	-1.95142112	-1.2927832	⋯	1.7990726	1.7141026	0.67053166	0.937460754	2.4784194	2.81268452	1.19497766	-2.28157500	-0.31754744	-2.0603525
LOC_OS12G01490	-3.10542573	-2.238916437	1.06881200	-3.0704562	2.8289928	-0.067515608	-0.24243608	0.4271979473	0.50406027	0.9324045	⋯	-0.4226213	0.8658668	-2.14219861	-1.675636454	-0.9441167	-1.05604115	2.77132707	2.24580557	0.33250002	-0.6657075
OSH42.OSKN2.HB365.OSH6.OSH15	-2.43670228	-2.638875692	-2.57776197	-1.8301494	2.7504047	3.177388910	2.85960022	0.0008800928	-0.31431380	-1.1225061	⋯	-1.4034796	1.2146787	-1.91813615	-2.307567913	-1.6039481	-1.71850310	3.16640737	2.25771791	0.25462855	1.5993663
OSMYBS3.MID1.LOC_OS01G09280.LOC_OS05G10690	-1.48507971	2.135759616	0.59234414	1.5302700	-0.6102463	0.755575997	1.69761528	0.6678913470	0.23837817	0.5193551	⋯	-1.1085107	-0.5642341	0.41324330	-2.221831745	0.7159034	0.28492738	0.01316692	-0.78941572	0.96187952	1.7910363
WOX1	0.31413217	0.300264285	-1.25688308	1.1288473	-1.5660192	1.691490569	1.83434205	-1.1784419853	-1.43416110	-0.7407892	⋯	0.1492971	-0.7833970	0.93767012	1.097108369	-0.4348394	0.91312156	-0.29346806	-1.18625485	-0.33573073	0.6102295
CCA1.LOC_OS02G45670.LOC_OS06G51260.OSMYB1R	-1.22312513	-0.032700798	-1.52483073	-0.1888229	1.1112853	3.351839761	3.53425548	-0.7065354864	-0.66145748	-0.9422159	⋯	-1.6560535	-0.8416171	-1.39195553	-1.355066881	0.3691453	0.27401870	0.96189214	1.17007918	0.83811287	1.6108615
OSPCF5.OSTCP21.OSTCP5	0.28352830	-0.193169282	-2.30020175	3.5375787	-1.0412777	-1.359691017	-1.81120044	-0.4635384262	-1.28297698	0.9851249	⋯	-2.1186590	-3.4654470	2.18150749	0.885491746	1.8634010	2.66437465	-1.21279661	0.66885189	-2.72001848	-1.7249163
OSDREB1G.OSERF035.OSARF7.OSERF018.OSERF007...	-3.06333411	-2.705595063	2.53297122	-2.6111596	2.3228171	0.001142265	-0.15591249	0.5474183308	0.68264787	0.4901331	⋯	1.4799065	1.8001250	-1.50574242	-1.398827232	-3.1357094	-3.13946573	3.29598139	3.17346216	-2.52002731	-1.3560995
OSSPL10	-1.71760051	-2.569321207	1.75381084	1.9808478	-0.4312404	-3.011570510	-3.06712059	3.3195438174	3.74413627	3.5638141	⋯	-0.3812606	1.2457183	2.08939665	2.073399810	-2.1455132	-1.60024755	2.29959760	0.02950445	-4.68277944	-4.1519410
OSHOX6.OSWOX8.OSHOX27	-3.03490965	-3.383388159	-0.25191704	-1.4822966	0.1539104	1.476001482	1.36202150	-0.8232932598	-0.84083632	-1.0730108	⋯	1.9045177	1.4855859	-1.23587678	-0.001948446	-1.7604704	-1.22085895	2.33613163	1.19057096	0.01172795	1.7887176
SPW1.RAG.MFO1.OSMADS55.OSMADS47...	-3.74473005	-3.632002314	-0.60499810	-3.4705656	0.5795229	-1.150902453	-0.96862209	0.0223327113	-0.06702799	0.2442982	⋯	-1.1810885	-0.3886685	-1.84920360	-1.673609444	1.7098710	4.28665201	5.08301691	-0.21290490	1.09170269	-1.7378830
OSMYB1R1.LOC_OS05G37730.OSMYBS3	-0.05544739	0.617341910	0.99423554	0.4040051	1.5561457	-0.610156167	-0.18807099	1.2773320270	1.29658285	0.2853807	⋯	-0.5330069	-0.3929779	-0.62128900	-1.769997817	0.3915312	-0.63065008	0.36950037	-0.13751059	-0.09811881	-0.2763111
EHD2.OSDOF16.OSDOF21.OSDOF7.OSMADS2	1.67736195	1.329282424	-1.82821505	1.9281510	-2.0236262	-0.439577831	-0.37788195	-1.0954256614	-0.92749766	-1.2744817	⋯	1.2200221	0.2884176	2.13134116	1.932704174	2.3781292	1.20161477	-2.92735298	-2.17672179	0.16311807	0.3175815
OSE2F1	0.77327345	0.325581142	0.68967153	-0.2112799	-1.0658642	2.293437317	1.41640137	-2.0937634373	-1.78301757	-2.5312563	⋯	1.8036983	1.0516836	0.81295613	0.911951095	1.7923290	1.49911624	-1.84460259	-0.78224950	-0.15725044	0.1713545
OSGATA6.OSGATA28.OSGATA2.OSGATA7.OSGATA3...	1.09409386	-0.792913710	-1.96554250	1.2340171	-0.8987149	-0.047514097	-0.08887095	-0.2504788887	-0.08944652	0.3158579	⋯	1.1883064	-0.3860344	0.55395705	0.837585494	-1.9924906	-0.18967743	-0.52630500	-1.58706888	1.43577455	-0.3110270
OSTCP7.OSTCP6.OSTCP19.OSTCP17.PCF1...	1.00535146	-1.248990575	0.52203553	-2.2534421	2.5951543	-3.228670200	-1.96621088	-1.5544828557	-0.51423404	-1.5303730	⋯	-3.1653262	-4.2774855	-0.21876838	-1.217809126	4.1797387	2.57288107	-1.64663092	1.57473419	2.62546247	4.6313778
OSHOX11.OSHOX17.DL	1.40479056	0.728433359	0.20444625	1.1106895	-1.6011388	-2.714430463	-2.61713472	0.6983380345	0.78981608	1.2212926	⋯	1.7519198	0.2500246	1.60110511	1.221368373	-2.7812869	-0.71902219	1.23741683	-1.29733049	-2.28137959	-1.7010967
OSCBT.LOC_OS03G09100.LOC_OS07G43030.LOC_OS10G22950	-0.41332567	-0.578677170	-0.24184343	-0.7697487	1.4951110	-0.178328628	-0.19183371	0.1021339907	0.34609918	-0.2977743	⋯	1.3823272	1.6827414	-0.76746417	-0.526893314	1.2944880	0.68412776	-0.11082932	2.69030694	0.01345706	-0.7056466
MADS63.OSMADS24.DEP.MADS14	1.30116052	0.692806585	0.34066174	-1.0716221	-0.6566346	0.719975801	1.23734058	-1.3278943858	-1.38707111	-1.4394319	⋯	-1.4414770	-2.4588587	-1.25167276	-1.400339957	2.0139566	3.87189337	2.47643176	-0.92373751	2.25606094	0.9759496
OSSPL9.OSSPL15.GW8.SPL6.WFP...	0.33402703	-0.002474019	-0.11154909	1.5372909	-0.4069566	-0.943488228	-1.10968865	3.0758538066	2.63703670	2.4000338	⋯	1.9448966	1.4461081	1.54558886	1.729839499	-4.7291246	-3.86968078	0.16645174	0.78443228	-3.93166943	-3.1505534
OSBZIP08.BZIP50.OSLG2.OSBZIP03.OSBZIP63...	-1.14597053	-2.509298733	-2.10035062	-1.1072140	0.3950251	-0.448258041	-0.10928916	0.3281068477	0.72230601	0.2359297	⋯	0.6771072	2.6268135	-0.55959940	-0.137361843	-0.4481356	-0.69209447	1.98922977	0.60238884	-1.32279133	-0.3265279
OSZHD7	-3.20333659	-2.030059162	1.99192520	-3.3622306	-1.5468512	0.410229165	0.80767842	-4.5521785091	-4.08199188	-4.0819919	⋯	4.9060708	4.5960128	-0.42461219	0.033857239	0.7688426	1.80341590	3.37675175	-0.84034472	-0.61267206	0.6511621
LOC_OS08G33050.LOC_OS02G47190.LOC_OS01G08160.OSMYBC.OSKANADI4...	0.88200297	0.067289692	0.42617587	-0.2623315	-0.4273818	0.853189917	0.69126637	0.8787743733	0.55584293	0.2449258	⋯	1.7854475	1.6649500	-0.04446175	-0.141394242	-0.7294843	-2.56692549	-2.04657480	0.17832863	-0.56314163	-0.3368944
OSERF044.OSERF008.LOC_OS06G09717.OSERF011.OSERF049...	-0.13412349	0.552569021	-0.01557689	-1.7991527	1.4316346	1.508998431	1.63071908	-0.0965769517	-0.26060396	-1.0603460	⋯	-0.1529669	-1.7920008	-2.84092075	-2.452280727	0.3258687	0.69963778	-0.31250632	1.02028163	2.37093969	1.6805057
ROC8.T..ROC4.ROC7.T..ROC5.T.	0.91583009	-1.534094337	0.75429614	-0.3565842	0.1556203	-1.348815343	-1.72902239	-0.1197560144	-0.38214585	0.1002826	⋯	1.1730987	0.6294062	-0.28363963	-1.497441553	1.2828239	1.42336776	-0.02832301	1.14083202	-0.10196114	-1.4200758
PHR1.LOC_OS07G48596.LOC_OS01G08160.LOC_OS06G40710.OSMYBC...	-0.20470671	-0.262954288	0.94210293	1.0090629	-1.0663063	-0.866673181	-0.43544093	-0.4513390614	-0.27301894	-0.5869105	⋯	1.4168838	-1.4339118	1.67311145	2.045698568	-1.3104017	0.44926620	-0.14695066	-1.39860527	0.75963535	1.3458452
OSHOX16.OSHOX4.OSHOX14.OSHOX20.OSHOX21...	3.21995232	4.384398840	-3.57635660	2.4189643	-1.0384231	-2.071435725	-2.45521598	-1.5012505568	-1.39690440	-0.3154339	⋯	2.1042018	-1.2293532	1.48518408	1.784670031	-3.7844225	-1.24767073	-1.06891033	-0.88421110	-1.10404203	-0.1630412
VP1..OSLFL1	4.27748554	4.454431679	-1.57749339	5.0830169	-5.0830169	-2.821030015	-3.92864970	-0.7552370968	-1.22002210	1.6283352	⋯	1.5559358	-2.1395951	5.08301691	5.083016914	0.8180333	2.45719658	-3.81222015	-5.08301691	-3.05900805	-4.9060708
LOC_OS06G04010.OSH1.LOC_OS10G02584	1.91918190	0.956924005	1.37933054	2.8409208	-1.1985844	-3.014577925	-2.72622925	2.2184317187	0.93834169	2.9992037	⋯	-1.9353700	-0.7501255	4.17973871	4.375232374	-2.0211927	-0.37925142	-1.76071973	-2.38888497	-2.10189420	-2.5463891
ERF130.OSERF123.DERF3.OSERF101.OSERF079...	0.61954023	0.839741042	1.68692938	-3.3319733	1.2674287	1.967696123	2.65542660	-0.6313982369	-1.61325570	-1.7612768	⋯	-0.7347624	0.6224690	-3.96136861	-4.242120573	0.5223799	0.03345495	0.04850089	1.73566235	1.07954679	1.4209759

In [36]:

options(repr.plot.width = 15, repr.plot.height =25, repr.plot.res = 100)
library(pheatmap)
blueYellow <- c("1"="#352A86","2"="#343DAE","3"="#0262E0","4"="#1389D2","5"="#2DB7A3","6"="#A5BE6A","7"="#F8BA43","8"="#F6DA23","9"="#F8FA0D")
solarExtra <- c("5"='#3361A5', "7"='#248AF3', "1"='#14B3FF', "8"='#88CEEF', "9"='#C1D5DC', "4"='#EAD397', "3"='#FDB31A',"2"= '#E42A2A', "6"='#A31D1D')
horizonExtra <- c("1"="#000436","4"="#021EA9","6"="#1632FB","8"="#6E34FC","3"="#C732D5","9"="#FD619D","7"="#FF9965","5"="#FFD32B","2"="#FFFC5A")
p5<- pheatmap(test,cellwidth=20,cellheight=15)
p5

Warning message:
“The input is a data frame, convert it to the matrix.”

Footprint R script¶

In [26]:

#!/usr/bin/env Rscript
# options(warn=-1)
# args = commandArgs(trailingOnly=TRUE)
# motif <- args[1]
# save <- args[2]
plot_footprint <- function(motif, save){
    library(genomation)
    library(EnrichedHeatmap)
    library(GenomicRanges)
    cal_footprint <- function(gene_bed, coverage_bed, group_name) { 
    gene <- readBed(gene_bed)
    coverage <- read.table(coverage_bed)
    colnames(coverage) <- c("chr","start","end","coverage")
    cov <- makeGRangesFromDataFrame(coverage,keep.extra.columns=T)
    mat1 = normalizeToMatrix(cov, gene, value_column = "coverage", 
                             extend = 200, mean_mode = "w0", w = 1,smooth = TRUE,background = NA,trim=c(0.01, 0.01))
    # group_name <- gsub(".bed","", basename(args[2]))
    # gene_name <- gsub(".bed","",basename(args[1]))
    quant_smo <- data.frame(x = c(1:length(colnames(mat1))), density=colMeans(mat1,na.rm=T),group=group_name)
    return(quant_smo)
}
    library(parallel)
    a <- list(sprintf("/public/workspace202011/encode/zhutao/analysis/TOBIAS/ATACorrect/panicle_correct/bindetect/%s/beds/%s_all.bed",motif,motif),
         sprintf("/public/workspace202011/encode/zhutao/analysis/TOBIAS/ATACorrect/panicle_correct/bindetect/%s/beds/%s_all.bed",motif,motif),
         sprintf("/public/workspace202011/encode/zhutao/analysis/TOBIAS/ATACorrect/panicle_correct/bindetect/%s/beds/%s_all.bed",motif,motif),
         sprintf("/public/workspace202011/encode/zhutao/analysis/TOBIAS/ATACorrect/panicle_correct/bindetect/%s/beds/%s_all.bed",motif,motif))

    b <- list("/public/workspace202011/encode/zhutao/analysis/TOBIAS/ATACorrect/panicle_correct/NIP_YP4_final.bed",
             "/public/workspace202011/encode/zhutao/analysis/TOBIAS/ATACorrect/panicle_correct/NIP_YP3_final.bed",
             "/public/workspace202011/encode/zhutao/analysis/TOBIAS/ATACorrect/panicle_correct/NIP_YP2_final.bed",
             "/public/workspace202011/encode/zhutao/analysis/TOBIAS/ATACorrect/panicle_correct/NIP_YP1_final.bed")
    c <- list("YP4","YP3","YP2","YP1")
    res <- list()
    res <- parallel::mcmapply(function(gene_bed, coverage_bed, group_name){
        return(cal_footprint(gene_bed, coverage_bed, group_name))
    }, gene_bed=a,coverage_bed=b,group_name=c,mc.cores=length(b))
    tt <- data.frame(x=c(res[[1]],res[[4]],res[[7]],res[[10]]), density=c(res[[2]],res[[5]],res[[8]],res[[11]]), group=c(res[[3]], res[[6]],res[[9]],res[[12]]))
    tt <- tt[tt$group==c("YP1","YP2","YP3","YP4"),]
    options(repr.plot.width = 7, repr.plot.height = 6, repr.plot.res = 100)
    p3 <- ggplot(tt, aes(x=x, y=density, color=group)) +
            geom_line(size=0.5)+#+theme(panel.background = element_rect(fill = "white", colour = "grey50"))+
            scale_color_manual(values = c("#ABD9E9","#74ADD1","#4575B4","#313695"))+scale_linetype_manual(values=c(rep("solid",5)))+
            theme_pubr(base_size = 12,base_family = "", border = TRUE, margin = TRUE, legend = c("right"), x.text.angle = 0)+
            geom_vline(xintercept = c(200-8,200+8), linetype="dashed", color = "#5C6068", size=0.5)+
            coord_cartesian(expand = T, ylim = c(0, 1.2*quantile(tt$density, 0.999)))+ scale_x_continuous(breaks=c(0,100,200,300,400),
   labels=c("-200", "-100", "0","100","200"))+labs(x = "Distance to motif center (bp)", y = 'Tn5 bias-corrected normalized insertions')
   #ggtitle(sprintf("%s_%s",unlist(strsplit(motif,"_"))[1],unlist(strsplit(motif,"_"))[2]))+theme(plot.title = element_text(hjust = 0.5))
   pdf(save,width=7,height=6)
   print(p3)
   dev.off()
   #p3
   return
}
# plot_footprint(motif, save)