search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
120 stars 30 forks source link

使用R而不是igv软件来展示bw文件 #3266

Closed ixxmu closed 1 year ago

ixxmu commented 1 year ago

https://mp.weixin.qq.com/s/UB6H3QS0dDpf7n1uVVSR1w

ixxmu commented 1 year ago

使用R而不是igv软件来展示bw文件 by 东林的扯淡小屋

###PLotEpiTrackByR###Author:wlt###Date:20211025###Replicate Figure 1b. in nature 2016 The landscape of accessible chromatin in###mammalian preimplantation embryos###Using GEO data GSE66581 for demonstration###Extra: using mESC from GSM1014187

###----0. Basic configuration----library(ggplot2)library(patchwork)library(cowplot)library(tidyverse)library(GEOquery)library(GenomicRanges)library(GenomicFeatures)library(org.Mm.eg.db)library(TxDb.Mmusculus.UCSC.mm10.knownGene)library(clusterProfiler)library(future.apply)library(rtracklayer)library(GENOVA)library(sigminer)library(RColorBrewer)

setwd("/data2/yudonglin/code/PlotEpiTrackByR-main")source("myfunction_lib.R")

####define colorbarsolarExtra <- colorRampPalette(c("#3361A5","#248AF3","#14B3FF","#88CEEF","#C1D5DC","#EAD397","#FDB31A","#E42A2A","#A31D1D"))hic.rdbu <- colorRampPalette(c( "#0065b2","#2399de","#67cdfe", "#f5f5f5","#fcb09d","#ed6855","#be2a21"))hic.sadle.rdbu <- colorRampPalette(colors = c("#1f67d1","#5299f2","#69b2fd","lightgrey","#ffc3aa","#ff876b","red1"))hic.contact.hot <- colorRampPalette(colors = GENOVA:::bezier_corrected_hot) hic.red <- colorRampPalette(c("white", "red"))hcl.pals(type = NULL)hcl.colors(200,palette = "YlGnBu")blueteal <- ggthemes::tableau_color_pal(palette = "Blue-Teal",type = "ordered-sequential")blueteal <- colorRampPalette(colors = blueteal(20))
scales::show_col(colours = hcl.colors(10,"YlGnBu"))
###---1. download data-----------###don't run it again
####you can donwload GEO supp data by getGEOSuppFilestmp.id <- "GSE66581"tmp.dir <- "data"dir.create(path = tmp.dir)getGEOSuppFiles(GEO = tmp.id,makeDirectory = T,baseDir = "data")#### unzip datauntar(tarfile="data/GSE66581/GSE66581_RAW.tar",exdir="data/GSE66581/byReplicate/bedgraph")
#### or directly donwload urltmp.input <- "https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1014nnn/GSM1014187/suppl/GSM1014187_mm9_wgEncodeUwDnaseEscj7S129ME0SigRep2.bigWig"tmp.output <- "data/ENCODE_mESC_DNaseI.bw" download.file(url=tmp.input,destfile=tmp.output)

安装软件和下载数据。

###----2. Plot gene Track-------library("TxDb.Mmusculus.UCSC.mm9.knownGene")###----2.1 Prepare gene track data------mm9 <- TxDb.Mmusculus.UCSC.mm9.knownGenemm9.tx <- myGetGRangesFromsTxDb(txdb = mm9,                                standard.chromosomes = T,                                verbose = T)tmp <- bitr(mm9.tx$gene_id,            fromType = "ENTREZID",            toType = "SYMBOL",            OrgDb = org.Mm.eg.db)tmp <- tmp[!duplicated(tmp$SYMBOL),]ttt.mm9.tx <- mm9.tx[mm9.tx$gene_id %in% tmp$ENTREZID]ttt.mm9.tx$gene_id <- plyr::mapvalues(ttt.mm9.tx$gene_id,                                      from = tmp$ENTREZID,                                      to = tmp$SYMBOL)mm9.tx <- ttt.mm9.tx
tmp.region <- "chr10:128222630-128407404"p_gene <- myGenePlot(annotation = mm9.tx,                     region = tmp.region,                     arrow_sbreaks = 600,                     font_size = 18,                     label_size = 3)+  ylab("Gene")+  theme(plot.margin = unit(c(0,0,0,0), "cm"),        axis.title.y = element_text(angle = 0,vjust = 0.5,hjust = 0.5),        axis.line = element_line(size = 1),        axis.ticks = element_line(size = 1))p_gene

mm10

mm9

###----2. plot DNase-seq track from bigwig---------
tmp.files <- "data/ENCODE_mESC_DNaseI.bw"tmp.colors <- "#446195"  tmp.levels <- "DNase-seq"  p_DNase <- myBigwigTrack(region = as(tmp.region,"GRanges"),                           bigwig = tmp.files,                           smooth = 100,                           lognorm = F,                           type = "coverage",                           y_label = tmp.levels,                           fontsize=18,                           track.color=tmp.colors,                           tmp.ylimits=c(0,100),                           max.downsample = 3000,                           downsample.rate = 0.1,                           tmp.seed=42)  p_DNase  ggsave(filename = "fig/DNase_coverage_plot.png",         width = 12,height = 8,dpi = 350,bg = "white")      p_DNase <- myBigwigTrack(region = as(tmp.region,"GRanges"),                           bigwig = tmp.files,                           smooth = 100,                           lognorm = F,                           type = "bar",                           y_label = tmp.levels,                           fontsize=18,                           track.color=tmp.colors,                           tmp.ylimits=c(0,100),                           max.downsample = 3000,                           downsample.rate = 0.1,                           tmp.seed=42)  p_DNase  ggsave(filename = "fig/DNase_bar_plot.png",         width = 12,height = 8,dpi = 350,bg = "white")  

  ####---------3.2 plot ATAC-seq track--------------    tmp.levels <- c("2cell","4cell","8cell","icm")  tmp.colors <- c("#3d7b42","#6a4181","#cf8e67","#a9477e")  scales::show_col(tmp.colors)    tmp.list.1 <- lapply(seq_along(tmp.levels),function(ii){    cat(ii,sep = "\n")    p1 <- myBigwigTrack(region = as(tmp.region,"GRanges"),                        bigwig = paste0("data/GSE66581/byReplicate/bigwig/",tmp.levels[ii],"_re1",".bw"),                        smooth = 100,                        lognorm = F,                        type = "bar",                        y_label = paste0(tmp.levels[ii]," rep1"),                        fontsize=18,                        track.color=tmp.colors[ii],                        tmp.ylimits=c(0,10),                        max.downsample = 3000,                        downsample.rate = 0.1,                        tmp.seed=42)    p2 <- myBigwigTrack(region = as(tmp.region,"GRanges"),                        bigwig = paste0("data/GSE66581/byReplicate/bigwig/",tmp.levels[ii],"_re2",".bw"),                        smooth = 100,                        lognorm = F,                        type = "bar",                        y_label =paste0(tmp.levels[ii]," rep2"),                        fontsize=18,                        track.color=tmp.colors[ii],                        tmp.ylimits=c(0,10),                        max.downsample = 3000,                        downsample.rate = 0.1,                        tmp.seed=42)    p <- list(p1,p2)    return(p)  })  tmp.list.1 <- Reduce(c,tmp.list.1)  p_ATAC_track_list <- tmp.list.1  plot.list <- c(p_ATAC_track_list,list(p_DNase))      (wrap_plots(plot.list,ncol = 1)  +       p_gene  + plot_layout(ncol = 1,heights = c(rep(1,12),3))) +     plot_annotation(title = "chr10:128,222,630–128,407,404") &     theme(plot.title = element_text(hjust = 1),          axis.line.x = element_blank(),          axis.text = element_blank(),          axis.ticks = element_blank(),          axis.title.x = element_blank())    ###create output directory  dir.create(path = "fig")  ggsave(filename = "fig/ATAC-seq_track_by_R.png",         width = 12,height = 8,         dpi = 350)      ####tmp.list.2  tmp.levels <- c("GSM1933924_2cell_re1","GSM1625847_4cell_re1","GSM1933928_8cell_re1")  tmp.colors <- c("#76492a","#76492a","#76492a")  scales::show_col(tmp.colors)    tmp.list.2 <- lapply(seq_along(tmp.levels),function(ii){    cat(ii,sep = "\n")    p <- myBigwigTrack(region = as(tmp.region,"GRanges"),                       bigwig = paste0("data/GSE66581/bigwig/",tmp.levels[ii],".bw"),                       smooth = 100,                       lognorm = F,                       type = "bar",                       y_label = tmp.levels[ii],                       fontsize=18,                       track.color=tmp.colors[ii],                       tmp.ylimits=c(0,10),                       max.downsample = 3000,                       downsample.rate = 0.1,                       tmp.seed=42)    return(p)  })    p_ATAC_track_list <- c(tmp.list.1,tmp.list.2)  plot.list <- c(p_ATAC_track_list,list(p_DNase))      (wrap_plots(plot.list,ncol = 1)  +       p_gene  + plot_layout(ncol = 1,heights = c(rep(1,12),3))) +     plot_annotation(title = "chr10:128,222,630–128,407,404") &     theme(plot.title = element_text(hjust = 1),          axis.line.x = element_blank(),          axis.text = element_blank(),          axis.ticks = element_blank(),          axis.title.x = element_blank())    ###create output directory  dir.create(path = "fig")  ggsave(filename = "fig/ATAC-seq_track_by_R.png",         width = 12,height = 8,         dpi = 350)

参考资料：

https://github.com/Landau1994/PlotEpiTrackByR