Unverified Commit eda69d4b authored by hy298's avatar hy298 Committed by GitHub
Browse files

Delete Volcanoplots_MAplots directory

parent b7c39adb

Volcanoplots_MAplots/ATAC_MAplots.R

deleted100644 → 0
+0 −32
Original line number Diff line number Diff line 
library(GenomicRanges)



### ATAC-seq MA plots



atac.new <- read.delim("ATAC_consensus_peaks.bed", sep="\t", header=FALSE,row.names=4)

getList <- function(fileIn,fileOut,plottitle){

    tmp0 <- read.delim(fileIn, sep="\t", header=TRUE, row.names=1)

    tmp1 = tmp0[row.names(tmp0) %in% row.names(atac.new),]

    pdf(fileOut, width = 4.5, height = 4.5)

    par(mar=c(4.5,4.8,2,1.7),mgp=c(2.8,1,0),bty="l",las=1)

    with(tmp1, plot(logCPM, logFC, pch=20, col="grey",xlim=c(0.5,8.5),

    ylim = c(-6, 6), cex=0.5, cex.lab = 1.5, cex.main=1.5,font.main=2,

    cex.axis=1.2, font.axis=2, font.lab=2, ylab="log2 FC",

    xlab="log2 CPM ATAC-seq",main=plottitle, axes=FALSE))

    box(lwd=3)

    axis(1,at=c(0,2,4,6,8),labels=NULL,lwd=3,font=2,cex.axis=1.5)

    axis(2,at=c(-6,-4,-2,0,2,4,6),labels=NULL,lwd=3,font=2,cex.axis=1.5)

    with(subset(tmp1, logFC >= 1), points(logCPM, logFC, pch=20, col="#db3236",cex=0.5))

    with(subset(tmp1, logFC <= -1), points(logCPM, logFC, pch=20, col="#4885ed",cex=0.5))

    up <- subset(tmp1, logFC >= 1)

    down <- subset(tmp1, logFC <= -1)

    Utext <- NROW(up)

    Dtext <- NROW(down)

    abline(h=c(1,-1),lty=2,lwd=2)

    legend("topright", x.intersp=-0.5, y.intersp=0.2, inset=0.01,legend=Utext, text.col="#db3236", cex=1.5, text.font=2,box.lty=0)

    legend("bottomright", x.intersp=-0.5, y.intersp=0.2, inset=0.01,legend=Dtext, text.col="#4885ed", cex=1.5, text.font=2,box.lty=0)

    dev.off()

}

getList("WT.NPM1.diff_ATAC.edgeR.FDR0.05.csv","WN.ATAC.diff_MAplot.pdf","Npm1c vs WT")

getList("WT.FLT3.diff_ATAC.edgeR.FDR0.05.csv","WF.ATAC.diff_MAplot.pdf","Npm1c vs WT")

getList("WT.DM.diff_ATAC.edgeR.FDR0.05.csv","WD.ATAC.diff_MAplot.pdf","Npm1c vs WT")

Volcanoplots_MAplots/Chromatin_modifications_MAplots.R

deleted100644 → 0
+0 −97
Original line number Diff line number Diff line 
library(GenomicRanges)



######################## Differential chromatin marks

## define atac peaks and enhancer H3K4me1 peaks

atac.summit_co20 <- read.table("ATAC.consensus.peak.summit.bed", sep="\t",header=F)

atac <- read.delim("ATAC_consensus_peaks_filter.bed", sep="\t", header=TRUE)

atac$id <- paste0(atac$chr,"_",atac$start,"_",atac$end)

atac.peakset=GRanges(seqnames=atac$chr, IRanges(atac$start,atac$end))



atac.ov.atac.summit_co20 = atac[atac$id %in% atac.summit_co20$id,]



atac.summit_peak = GRanges(seqnames=atac.ov.atac.summit_co20$chr,IRanges(atac.ov.atac.summit_co20$start,atac.ov.atac.summit_co20$end),id=atac.ov.atac.summit_co20$id)



### H3K4me1

getList <- function(fileIn,fileOut1,group){

    tmp=read.delim(fileIn, sep="\t", header=TRUE, row.names=1)

    tmp1=strsplit(as.character(row.names(tmp)), "_") ### starsplit can only split character factors

    tmp2=GRanges(seqnames=sapply(tmp1, "[", 1), IRanges(as.numeric(sapply(tmp1, "[", 2)), as.numeric(sapply(tmp1, "[", 3))),

    logFC=tmp$logFC,logCPM=tmp$logCPM,FDR=tmp$FDR)

    tmp3=tmp2[tmp2 %over% atac.summit_peak,]

    tmp4=data.frame(seqnames(tmp3),start(tmp3),end(tmp3),tmp3$logFC,tmp3$logCPM,tmp3$FDR,

    stringsAsFactors = FALSE)

    colnames(tmp4)=c("chr","start","end","logFC","logCPM","FDR")

    ## MA plots

    pdf(fileOut1, width = 4.5, height = 4.5)

    par(mar=c(4.5,4.5,3,1),mgp=c(2.8,1,0),bty="l",las=1)

    with(tmp4, plot(logCPM, logFC, pch=20, col="grey",xlim=c(1,7),

    ylim = c(-4, 4),cex=1, cex.lab = 1.5, cex.main=1.5,font.main=2,

    cex.axis=1.2, font.axis=2,font.lab=2,

    xlab="log2 CPM H3K4me1",ylab="log2 FC",main=group,axes=FALSE))

    box(lwd=3)

    axis(1,at=c(2,4,6),labels=NULL,lwd=3,font=2,cex.axis=1.5)

    axis(2,at=c(-4,-2,0,2,4),labels=NULL,lwd=3,font=2,cex.axis=1.5)

    with(subset(tmp4, logFC >= log(1.5,2)), points(logCPM, logFC, pch=20, col="#db3236",cex=1))

    with(subset(tmp4, logFC <= -log(1.5,2)), points(logCPM, logFC, pch=20, col="#4885ed",cex=1))

    up <- subset(tmp4, logFC >= log(1.5,2))

    down <- subset(tmp4, logFC <= -log(1.5,2))

    Utext <- NROW(up)

    Dtext <- NROW(down)

    legend("topright", x.intersp=-0.5, y.intersp=0.2, inset=0.01,legend=Utext, text.col="#db3236", cex=1.5, text.font=2,box.lty=0)

    legend("bottomright", x.intersp=-0.5, y.intersp=0.2, inset=0.01,legend=Dtext, text.col="#4885ed", cex=1.5, text.font=2,box.lty=0)

    abline(h=c(-log(1.5,2),log(1.5,2)), lty=2,lwd=2)

    dev.off()

}

getList("WT.NPM1.diff_H3K4me1.edgeR.FDR0.05.csv","WN.H3K4me1.diff_MAplot.pdf","NPM1 vs WT")

getList("WT.FLT3.diff_H3K4me1.edgeR.FDR0.05.csv","WF.H3K4me1.diff_MAplot.pdf","FLT3 vs WT")

getList("WT.DM.diff_H3K4me1.edgeR.FDR0.05.csv","WD.H3K4me1.diff_MAplot.pdf","DM vs WT")







## define atac peaks and enhancer H3K27ac peaks

enh.all <- read.delim("enhancerPeaks.bed", sep="\t", header=T)

enh.all.peakset <- GRanges(seqnames = enh.all$chr, ranges = IRanges(enh.all$start,enh.all$end))



atac.summit.ov.enh = atac.summit_co20[atac.summit_co20$id %in% atac.summit_peak.ov.enh$id, ]

atac.summit.ov.enh_peak = GRanges(seqnames = atac.summit.ov.enh$chr,ranges = IRanges(atac.summit.ov.enh$start,atac.summit.ov.enh$end),id = atac.summit.ov.enh$id)



atac.summit500bp.ov.enh = resize(atac.summit.ov.enh_peak, width = 500, fix = "center")



enh.ov.atac.summit500 = enh.all.peakset[enh.all.peakset %over% atac.summit500bp.ov.enh, ]



#### H3K27ac

getList <- function(fileIn,fileOut1,group){

  tmp=read.delim(fileIn, sep="\t", header=TRUE, row.names=1)

  tmp1=strsplit(as.character(row.names(tmp)), "_") ### starsplit can only split character factors

  tmp2=GRanges(seqnames=sapply(tmp1, "[", 1), IRanges(as.numeric(sapply(tmp1, "[", 2)), as.numeric(sapply(tmp1, "[", 3))),

               logFC=tmp$logFC,logCPM=tmp$logCPM,FDR=tmp$FDR)

  tmp3=tmp2[tmp2 %over% enh.ov.atac.summit500,]

  tmp4=data.frame(seqnames(tmp3),start(tmp3),end(tmp3),tmp3$logFC,tmp3$logCPM,tmp3$FDR,

                  stringsAsFactors = FALSE)

  colnames(tmp4)=c("chr","start","end","logFC","logCPM","FDR")

  ## MA plots

  pdf(fileOut1, width = 4.5, height = 4.5)

  par(mar=c(4.5,4.5,3,1),mgp=c(2.8,1,0),bty="l",las=1)

  with(tmp4, plot(logCPM, logFC, pch=20, col="grey",xlim=c(1.5,8.5), 

                  ylim = c(-4.5, 4.5),cex=1, cex.lab = 1.5, cex.main=1.5,font.main=2,

                  cex.axis=1.2, font.axis=2,font.lab=2, 

                  xlab="log2 CPM H3K27ac Enh",ylab="log2 FC",main=group,axes=FALSE))

  box(lwd=3)

  axis(1,at=c(2,4,6,8),labels=NULL,lwd=3,font=2,cex.axis=1.5)

  axis(2,at=c(-4,-2,0,2,4),labels=NULL,lwd=3,font=2,cex.axis=1.5)

  with(subset(tmp4, logFC >= log(1.5,2)), points(logCPM, logFC, pch=20, col="#db3236",cex=1))

  with(subset(tmp4, logFC <= -log(1.5,2)), points(logCPM, logFC, pch=20, col="#4885ed",cex=1))

  up <- subset(tmp4, logFC >= log(1.5,2))

  down <- subset(tmp4, logFC <= -log(1.5,2))

  Utext <- NROW(up)

  Dtext <- NROW(down)

  legend("topright", x.intersp=-0.5, y.intersp=0.2, inset=0.01,legend=Utext, text.col="#db3236", cex=1.5, text.font=2,box.lty=0)

  legend("bottomright", x.intersp=-0.5, y.intersp=0.2, inset=0.01,legend=Dtext, text.col="#4885ed", cex=1.5, text.font=2,box.lty=0)

  abline(h=c(-log(1.5,2),log(1.5,2)), lty=2,lwd=2)

  dev.off()

}

getList("WT.NPM1.diff_H3K27ac.edgeR.FDR0.05.csv","WN.H3K27ac.diff_MAplot.pdf","NPM1 vs WT")

getList("WT.FLT3.diff_H3K27ac.edgeR.FDR0.05.csv","WF.H3K27ac.diff_MAplot.pdf","FLT3 vs WT")

getList("WT.DM.diff_H3K27ac.edgeR.FDR0.05.csv","WD.H3K27ac.diff_MAplot.pdf","DM vs WT")

Volcanoplots_MAplots/README.md

deleted100644 → 0
+0 −14
Original line number Diff line number Diff line 
## Volcano plots for RNA-seq: input files were converted from DEseq2 differential anlaysis

RNAseq_volcanoplots.R





## MA plots for ChIPseq and ATACseq data: input files were converted from edgeR differential anlaysis

Chromatin_modification_MAplots.R

ATAC_MAplots.R

Volcanoplots_MAplots/RNAseq_volcanoplots.R

deleted100644 → 0
+0 −33
Original line number Diff line number Diff line 
setwd("~/Documents/20180920_manuscript/RNA-seq")

library(ggplot2)



getList <- function(fileIn, fileOut1, maintitle){

  tmp <- read.delim(fileIn, sep="\t",header=TRUE)

  tmp.up1=subset(tmp, padj < 0.05 & log2FC >= log(1.5,2))

  tmp.down1=subset(tmp, padj < 0.05 & log2FC <= -log(1.5,2))

  ### scatter plot

  tmp$padj[tmp$padj < 1e-100 ] <- 1e-100

  tmp[tmp == NA ] <- 1

  tmp.up=subset(tmp, padj < 0.05 & log2FC >= log(1.5,2))

  tmp.down=subset(tmp, padj < 0.05 & log2FC <= -log(1.5,2))

  pdf(fileOut1, width = 4.5, height = 4.5)

  #png(fileOut1, res=300, units = "in", width = 4.5, height = 4.5)

  par(mar=c(4.5,5,3,1),mgp=c(2.8,1,0),bty="l",las=1)

  with(tmp, plot(log2FC, -log10(padj), pch=20, col="grey", xlim=c(-10, 10), ylim=c(0,115),

                     cex=1, cex.lab = 1.5, cex.axis=1.5, font.axis=2,font.lab=2, 

                     xlab="log2 FC RNA-seq", ylab="-log10 adjP",axes=FALSE, main=maintitle, cex.main=1.6, font.main=2))

  box(lwd=3)

  axis(1,labels=TRUE,lwd=3,font=2,cex.axis=1.5)

  axis(2,labels=TRUE,lwd=3,font=2,cex.axis=1.5)

  with(tmp.up, points(log2FC, -log10(padj), pch=20, col="#db3236",cex=1))

  with(tmp.down, points(log2FC, -log10(padj), pch=20, col="#4885ed",cex=1))

  Utext <- nrow(tmp.up)

  Dtext <- nrow(tmp.down)

  legend("toprigh", x.intersp=-0.5, y.intersp=0.2, inset=0.01,legend=Utext, text.col="#db3236", cex=1.5, text.font=2,box.lty=0)

  legend("topleft", x.intersp=-0.5, y.intersp=0.2, inset=0.01,legend=Dtext, text.col="#4885ed", cex=1.5, text.font=2,box.lty=0)

  abline(v=c(log(1.5,2),-log(1.5,2)), h=-log(0.05,10), lty=2,lwd=2)

  dev.off()

}

getList("WT_NPM1_de_forvolcano.txt", "WN.diffExp.volcano.pdf", "NPM1 vs WT")

getList("WT_FLT3_de_forvolcano.txt", "WF.diffExp.volcano.pdf", "FLT3 vs WT")

getList("WT_DM_de_forvolcano.txt", "WD.diffExp.volcano.pdf", "DM vs WT")
 
 No newline at end of file

CRA Git | Maintained and supported by SUSTech CRA and CCSE