From d14dd1a290866e35694babdf3815e1ef16b00be0 Mon Sep 17 00:00:00 2001 From: Source Codes for Computational and Theoretical Life Science Date: Thu, 4 Aug 2022 01:36:33 -0500 Subject: [PATCH] Add files via upload --- Call_Mutations/figures/2-getFraction.pl | 176 ++ Call_Mutations/figures/3-DensityPlot.R | 380 ++++ Call_Mutations/figures/5-number.sh | 5 + Call_Mutations/figures/6-overlap.sh | 76 + Call_Mutations/figures/7-getMotifs.pl | 72 + Call_Mutations/figures/7-histogram.R | 60 + Call_Mutations/figures/Anno_v20210106.R | 2180 ++++++++++++++++++++++ Call_Mutations/figures/Anno_v20210106.sh | 14 + Call_Mutations/figures/figures.pl | 118 ++ Call_Mutations/figures/heatmap.R | 157 ++ Call_Mutations/figures/histogram.R | 39 + 11 files changed, 3277 insertions(+) create mode 100644 Call_Mutations/figures/2-getFraction.pl create mode 100644 Call_Mutations/figures/3-DensityPlot.R create mode 100644 Call_Mutations/figures/5-number.sh create mode 100644 Call_Mutations/figures/6-overlap.sh create mode 100644 Call_Mutations/figures/7-getMotifs.pl create mode 100644 Call_Mutations/figures/7-histogram.R create mode 100644 Call_Mutations/figures/Anno_v20210106.R create mode 100644 Call_Mutations/figures/Anno_v20210106.sh create mode 100644 Call_Mutations/figures/figures.pl create mode 100644 Call_Mutations/figures/heatmap.R create mode 100644 Call_Mutations/figures/histogram.R diff --git a/Call_Mutations/figures/2-getFraction.pl b/Call_Mutations/figures/2-getFraction.pl new file mode 100644 index 0000000..db095b2 --- /dev/null +++ b/Call_Mutations/figures/2-getFraction.pl @@ -0,0 +1,176 @@ +#!/usr/bin/env perl +use strict; +use warnings; +use v5.12; +## Perl5 version >= 5.12 +## Suffixes of all self-defined global variables must be "_g". +################################################################################################################################################################################################### + + + + + +################################################################################################################################################################################################### +my $input_g = ''; ## such as "1-finalSites/HEK293.bed", global variable. +my $output_g = ''; ## such as "2-getFraction/HEK293.bed", global variable. +my $file_g = ''; ## such as "A_C_merge.spikeins.fit.txt", global variable. + +{ +## Help Infromation +my $HELP = ' + ------------------------------------------------------------------------------------------------------------------------------------------------------ + ------------------------------------------------------------------------------------------------------------------------------------------------------ + Usage: + perl 2-getFraction.pl [-version] [-help] [-f file] [-in inputFile] [-out outFile] + For instance: + perl 2-getFraction.pl -f A_C_merge.spikeins.fit.txt -in 1-finalSites/HEK293.bed -out 2-getFraction/HEK293.bed + + Optional arguments: + -version Show version number of this program and exit. + -help Show this help message and exit. + + Required arguments: + -f file Linear regression results. + -in inputFile "inputFile" is the name of your input file. (no default) + -out outFile "outFile" is the name of output file. (no default) + + ------------------------------------------------------------------------------------------------------------------------------------------------------ + ------------------------------------------------------------------------------------------------------------------------------------------------------ +'; + +## Version Infromation +my $version = " Version 0.2, 2021-04-01."; + +## Keys and Values +if ($#ARGV == -1) { say "\n$HELP\n"; exit 0; } ## when there are no any command argumants. +if ($#ARGV%2 == 0) { @ARGV = (@ARGV, "-help") ; } ## when the number of command argumants is odd. +my %args = @ARGV; + +## Initialize Variables +$input_g = '1-finalSites/HEK293.bed'; ## This is only an initialization value or suggesting value, not default value. +$output_g = '2-getFraction/HEK293.bed'; ## This is only an initialization value or suggesting value, not default value. + +## Available Arguments +my $available = " -version -help -in -out -f "; +my $boole = 0; +while( my ($key, $value) = each %args ) { + if ( ($key =~ m/^\-/) and ($available !~ m/\s$key\s/) ) {say "\n\tCann't recognize $key"; $boole = 1; } +} +if($boole == 1) { + say "\tThe Command Line Arguments are wrong!"; + say "\tPlease see help message by using 'perl 2-getFraction.pl -help' \n"; + exit 0; +} + +## Get Arguments +if ( exists $args{'-version' } ) { say "\n$version\n"; exit 0; } +if ( exists $args{'-help' } ) { say "\n$HELP\n"; exit 0; } +if ( exists $args{'-f' } ) { $file_g = $args{'-f' }; }else{say "\n -f is required.\n"; say "\n$HELP\n"; exit 0; } +if ( exists $args{'-in' } ) { $input_g = $args{'-in' }; }else{say "\n -in is required.\n"; say "\n$HELP\n"; exit 0; } +if ( exists $args{'-out' } ) { $output_g = $args{'-out' }; }else{say "\n -out is required.\n"; say "\n$HELP\n"; exit 0; } + +## Conditions +$input_g =~ m/^\S+$/ || die "\n\n$HELP\n\n"; +$output_g =~ m/^\S+$/ || die "\n\n$HELP\n\n"; + +## Print Command Arguments to Standard Output +say "\n + ################ Arguments ############################### + File: $file_g + Input File: $input_g + Output File: $output_g + ########################################################## +\n"; +} +################################################################################################################################################################################################### + + + + + +################################################################################################################################################################################################### +say "\n\n\n\n\n\n##################################################################################################"; +say "Running ......"; + +sub myMakeDir { + my $path = $_[0]; + if ( !( -e $path) ) { system("mkdir -p $path"); } + if ( !( -e $path) ) { mkdir $path || die; } +} + + +&myMakeDir($output_g); +system( "rm -rf $output_g" ); + +my $file_a = $file_g; +open(fileIn, "<", $input_g ) or die; +open(fileAa, "<", $file_a ) or die; +open(fileO1, ">", $output_g) or die; + +my @lines1 = ; +my @lines2 = ; +################################################################################################################################################################################################### + + + + + +################################################################################################################################################################################################### +my @motifs = @lines2; +my @coefficients1 = @lines2; +my @coefficients2 = @lines2; +for(my $i=0; $i<=$#lines2; $i=$i+1) { + my $temp = $lines2[$i]; + $temp =~ m/^(\S+)\s+([e\-\.\d]+)\s+([e\-\.\d]+)\n/ or die "##$temp##"; + $motifs[$i] = $1; + $coefficients1[$i] = $2; + $coefficients2[$i] = $3; +} + + +for(my $i=0; $i<=$#lines1; $i=$i+1) { + my $temp = $lines1[$i]; + $temp =~ m/^(chr\S+)\s+(\d+)\s+(\d+)\s+(\S+)\s+([\.\d]+)\s+([-+])\s/ or die "##$temp##"; + my $chrom = $1; + my $start = $2; + my $end = $3; + my $name = $4; + my $mutation = $5; + my $strand = $6; + $mutation = $mutation/100; + next unless $chrom =~ m/^chr[\dXYM]+$/; + + $name =~ m/\.{6}([AGCT]{5})/ or die "##$name\n##"; + my $fiveMer = $1; + + my $bool = 0; + my $fraction = -1; + + for(my $j=0; $j<=$#lines2; $j=$j+1) { + if($motifs[$j] eq $fiveMer) {$fraction = ($mutation-$coefficients2[$j])/$coefficients1[$j]; $bool=$bool+1; } + } + + ($bool == 1) or die "##$bool##$fiveMer##\n"; + if($fraction < 0.01) {$fraction = 0.01}; + $fraction = $fraction*100; + $mutation = $mutation*100; + if($fraction > 100) {$fraction = 100; } + + print fileO1 "$chrom\t$start\t$end\t$name\t$mutation\t$strand\t$fraction\n"; +} + +################################################################################################################################################################################################### + + + + + + +say "Done ......"; + + + + + + +##### diff --git a/Call_Mutations/figures/3-DensityPlot.R b/Call_Mutations/figures/3-DensityPlot.R new file mode 100644 index 0000000..7e3ad9c --- /dev/null +++ b/Call_Mutations/figures/3-DensityPlot.R @@ -0,0 +1,380 @@ +################################################################################################################## +library(ggplot2) +library(randomcoloR) +################################################################################################################## + + + + + + +library("ggplot2") +library("reshape2") +library("psych") +library("minerva") +library("plyr") + +MyTheme_1_g <- function(textSize1=14, hjust1=NULL, vjust1=NULL, angle1=NULL) { # "hjust=1, vjust=1, angle=30" for some boxplots. + theme( + line = element_line(colour="black", size=0.1, linetype=1, lineend=NULL), ## all line elements. 局部优先总体,下面3个也是,只对非局部设置有效. 所有线属性. + rect = element_rect(colour="black", size=0.1, linetype=1, fill="transparent" ), ## all rectangluar elements. hjust=1: 靠右对齐. 所有矩形区域属性. + text = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## all text elements. "serif" for a serif font. 所有文本相关属性. + title = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## all title elements: plot, axes, legends. hjust:水平对齐的方向. 所有标题属性. + ## aspect.ratio = 1, ##高宽比 + + axis.title = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## label of axes (element_text; inherits from text). horizontal: 水平的, 水平线 + axis.title.x = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## x axis label (element_text; inherits from axis.title) + axis.title.y = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=90, lineheight=1.0, margin = NULL, debug = NULL), ## y axis label (element_text; inherits from axis.title) + axis.text = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## tick labels along axes (element_text; inherits from text). 坐标轴刻度的标签的属性. + axis.text.x = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=hjust1, vjust=vjust1, angle=angle1, lineheight=1.0, margin = NULL, debug = NULL), ## x axis tick labels (element_text; inherits from axis.text) + axis.text.y = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## y axis tick labels (element_text; inherits from axis.text) + + axis.ticks = element_line(colour="black", size=0.1, linetype=1, lineend=NULL), ## tick marks along axes (element_line; inherits from line). 坐标轴刻度线. + axis.ticks.x = element_line(colour="black", size=0.1, linetype=1, lineend=NULL), ## x axis tick marks (element_line; inherits from axis.ticks) + axis.ticks.y = element_line(colour="black", size=0.1, linetype=1, lineend=NULL), ## y axis tick marks (element_line; inherits from axis.ticks) + axis.ticks.length = grid::unit(1.0, "mm", data=NULL), ## length of tick marks (unit), ‘"mm"’ Millimetres. 10 mm = 1 cm. 刻度线长度 + axis.line = element_line(colour="transparent", size=0.1, linetype=1, lineend=NULL), ## lines along axes (element_line; inherits from line). 坐标轴线 + axis.line.x = element_line(colour="transparent", size=0.1, linetype=1, lineend=NULL), ## line along x axis (element_line; inherits from axis.line) + axis.line.y = element_line(colour="transparent", size=0.1, linetype=1, lineend=NULL), ## line along y axis (element_line; inherits from axis.line) + + legend.background = element_rect(colour="transparent", size=0.1, linetype=1, fill="transparent" ), ## background of legend (element_rect; inherits from rect) + legend.spacing = grid::unit(0.1, "mm", data=NULL), ## extra space added around legend (unit). linetype=1指的是矩形边框的类型. + legend.key = element_rect(colour="transparent", size=0.1, linetype=1, fill="transparent" ), ## background underneath legend keys. 图例符号. size=1指的是矩形边框的大小. + legend.key.size = grid::unit(1, "mm", data=NULL) , ## size of legend keys (unit; inherits from legend.key.size) + legend.key.height = grid::unit(1, "mm", data=NULL) , ## key background height (unit; inherits from legend.key.size) + legend.key.width = grid::unit(1, "mm", data=NULL) , ## key background width (unit; inherits from legend.key.size) + legend.text = element_text(family="serif", face=NULL, colour="black", size=textSize1, hjust=NULL, vjust=NULL, angle=NULL, lineheight=NULL), ##legend item labels. 图例文字标签. + legend.text.align = 0, ## alignment of legend labels (number from 0 (left) to 1 (right)) + legend.title = element_blank(), ## title of legend (element_text; inherits from title) + legend.title.align = 0, ## alignment of legend title (number from 0 (left) to 1 (right)) + legend.position = "right", ## the position of legends. ("left", "right", "bottom", "top", or two-element numeric vector) + legend.direction = "vertical", ## layout of items in legends ("horizontal" or "vertical") 图例排列方向 + legend.justification = "center", ## anchor point for positioning legend inside plot ("center" or two-element numeric vector) 图例居中方式 + legend.box = NULL, ## arrangement of multiple legends ("horizontal" or "vertical") 多图例的排列方式 + legend.box.just = NULL, ## justification of each legend within the overall bounding box, when there are multiple legends ("top", "bottom", "left", or "right") 多图例的居中方式 + + panel.background = element_rect(colour="transparent", size=0.0, linetype=1, fill="transparent" ), ## background of plotting area, drawn underneath plot (element_rect; inherits from rect) + panel.border = element_rect(colour="black", size=0.1, linetype=1, fill=NA ), ## border around plotting area, drawn on top of plot so that it covers tick marks and grid lines. This should be used with fill=NA (element_rect; inherits from rect) + panel.spacing = grid::unit(1, "mm", data=NULL) , ## margin around facet panels (unit) 分面绘图区之间的边距 + panel.spacing.x = grid::unit(1, "mm", data=NULL) , + panel.spacing.y = grid::unit(1, "mm", data=NULL) , + panel.grid = element_blank(), ## grid lines (element_line; inherits from line) 绘图区网格线 + panel.grid.major = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## major grid lines (element_line; inherits from panel.grid) 主网格线 + panel.grid.minor = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## minor grid lines (element_line; inherits from panel.grid) 次网格线 + panel.grid.major.x = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## vertical major grid lines (element_line; inherits from panel.grid.major) + panel.grid.major.y = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## horizontal major grid lines (element_line; inherits from panel.grid.major) + panel.grid.minor.x = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## vertical minor grid lines (element_line; inherits from panel.grid.minor) + panel.grid.minor.y = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## horizontal minor grid lines (element_line; inherits from panel.grid.minor) + + plot.background = element_rect(colour="transparent", size=NULL, linetype=NULL, fill="transparent" ), ## background of the entire plot (element_rect; inherits from rect) 整个图形的背景 + plot.title = element_text(family="serif", face=NULL, colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=NULL, lineheight=NULL), ## plot title (text appearance) (element_text; inherits from title) 图形标题 + plot.margin = grid::unit(c(0.2, 0.2, 0.2, 0.2), "mm", data=NULL), ## margin around entire plot (unit with the sizes of the top, right, bottom, and left margins) + + strip.background = element_rect(colour=NULL, size=NULL, linetype=NULL, fill=NULL ), ## background of facet labels (element_rect; inherits from rect) 分面标签背景 + strip.text = element_text(family="serif", face=NULL, colour=NULL, size=NULL, hjust=NULL, vjust=NULL, angle=NULL, lineheight=NULL), ## facet labels (element_text; inherits from text) + strip.text.x = element_text(family="serif", face=NULL, colour=NULL, size=NULL, hjust=NULL, vjust=NULL, angle=NULL, lineheight=NULL), ## facet labels along horizontal direction (element_text; inherits from strip.text) + strip.text.y = element_text(family="serif", face=NULL, colour=NULL, size=NULL, hjust=NULL, vjust=NULL, angle=NULL, lineheight=NULL) ## facet labels along vertical direction (element_text; inherits from strip.text) + ) +} + + +MySaveGgplot2_1_g <- function(ggplot2Figure1, path1, fileName1, height1, width1) { + SVG1 <- paste(path1, "/", "SVG", sep = "", collapse = NULL) + PNG1 <- paste(path1, "/", "PNG", sep = "", collapse = NULL) + PDF1 <- paste(path1, "/", "PDF", sep = "", collapse = NULL) + #EPS1 <- paste(path1, "/", "EPS", sep = "", collapse = NULL) + if( ! file.exists(SVG1) ) { dir.create(SVG1) } + if( ! file.exists(PNG1) ) { dir.create(PNG1) } + if( ! file.exists(PDF1) ) { dir.create(PDF1) } + #if( ! file.exists(EPS1) ) { dir.create(EPS1) } + ggsave( filename = paste(SVG1, "/", fileName1, ".svg", sep="", collapse=NULL), height=height1, width=width1, dpi = 1200 ) + ggsave( filename = paste(PNG1, "/", fileName1, ".png", sep="", collapse=NULL), height=height1, width=width1, dpi = 1200 ) + ggsave( filename = paste(PDF1, "/", fileName1, ".pdf", sep="", collapse=NULL), height=height1, width=width1, dpi = 1200 ) + #ggsave( filename = paste(EPS1, "/", fileName1, ".eps", sep="", collapse=NULL), height=height1, width=width1, dpi = 1200, device=cairo_ps) +} + + +## df contains two columns, the first column (cond_col=1) is sample type, the second column (val_col=2) is value. (must be). +whisk_1_g <- function(df, cond_col=1, val_col=2) { + require(reshape2) + condname <- names(df)[cond_col] ## save the name of the first column. + names(df)[cond_col] <- "cond" + names(df)[val_col] <- "value" + b <- boxplot(value~cond, data=df, plot=FALSE) + df2 <- cbind(as.data.frame(b$stats), c("min","lq","m","uq","max")) + names(df2) <- c(levels(df$cond), "pos") + df2 <- melt(df2, id="pos", variable.name="cond") + df2 <- dcast(df2, cond~pos) + names(df2)[1] <- condname + print(df2) + df2 +} + + +MyBoxViolinPlot_1_g <- function(vector2, sampleType2, sampleRank2, colours2, path2, fileName2, title2, xLab2, yLab2, height2=4, width2=4, Ymin2=0, Ymax2=3) { + vector2[vector2>Ymax2] <- Ymax2 + vector2[vector2Xmax2] <- Xmax2 + vector2[vector2xMax2] <- xMax2 + vector2[vector2 $out/number.txt + diff --git a/Call_Mutations/figures/6-overlap.sh b/Call_Mutations/figures/6-overlap.sh new file mode 100644 index 0000000..7521af9 --- /dev/null +++ b/Call_Mutations/figures/6-overlap.sh @@ -0,0 +1,76 @@ +A1_1="2-getFraction/sc0.5ng.bed"; + +A1_3="2-getFraction/sc5ng.bed"; + +A1_5="2-getFraction/sc50ng.bed"; + + +A1_1a="2-getFraction/HeLa.m6A.Peaks.bed"; +A1_2a="2-getFraction/Lulu_Hela.bed"; +A1_3a="2-getFraction/Lulu_Hela_RIP.bed"; +A1_4a="2-getFraction/A1_11_common.bed"; +A1_5a="2-getFraction/A2_11_common.bed"; +A1_6a="2-getFraction/As1_11_common.bed"; +A1_7a="2-getFraction/As2_11_common.bed"; +A1_8a="2-getFraction/Ass1_11_common.bed"; +A1_9a="2-getFraction/Ass2_11_common.bed"; + + +################################## + + + +outpath1="6-overlap/pairwise.fraction" +mkdir -p $outpath1 + +intervene pairwise \ +--input \ +$A1_1 \ +$A1_3 \ +$A1_5 \ +$A1_1a \ +$A1_2a \ +$A1_3a \ +$A1_4a \ +$A1_5a \ +$A1_6a \ +$A1_7a \ +$A1_8a \ +$A1_9a \ +--compute frac --htype color --output $outpath1 \ + --dpi 1200 --figtype svg --figsize 10 10 --fontsize 10 + + + + + + + + + +outpath2="6-overlap/pairwise.count" +mkdir -p $outpath2 + +intervene pairwise \ +--input \ +$A1_1 \ +$A1_3 \ +$A1_5 \ +$A1_1a \ +$A1_2a \ +$A1_3a \ +$A1_4a \ +$A1_5a \ +$A1_6a \ +$A1_7a \ +$A1_8a \ +$A1_9a \ +--compute count --htype color --output $outpath2 \ + --dpi 1200 --figtype svg --figsize 10 10 --fontsize 10 + + + + + + + diff --git a/Call_Mutations/figures/7-getMotifs.pl b/Call_Mutations/figures/7-getMotifs.pl new file mode 100644 index 0000000..a90245b --- /dev/null +++ b/Call_Mutations/figures/7-getMotifs.pl @@ -0,0 +1,72 @@ +#!/usr/bin/env perl +use strict; +use warnings; +use v5.12; +## Perl5 version >= 5.12 +## Suffixes of all self-defined global variables must be "_g". +################################################################################################################################################################################################### + + + +my $input1="2-getFraction"; +my $out="7-getMotifs"; + + +################################################################################################################################################################################################### +say "\n\n\n\n\n\n##################################################################################################"; +say "Running ......"; + +sub myMakeDir { + my $path = $_[0]; + if ( !( -e $path) ) { system("mkdir -p $path"); } + if ( !( -e $path) ) { mkdir $path || die; } +} +&myMakeDir($out); + +opendir(my $FH_input1_g, $input1) || die; +my @files = readdir($FH_input1_g); +################################################################################################################################################################################################### + + + + + +################################################################################################################################################################################################### + +for(my $i1=0; $i1<=$#files; $i1=$i1+1) { + my $temp1 = $files[$i1]; + next unless $temp1 =~ m/\.bed$/; + print "$temp1\n"; + $temp1 =~ s/\.bed$// or die; + + open(my $FH1_g, "$input1/$temp1.bed" ) || die; + my @lines1 = <$FH1_g> ; + open(fileA, ">", "$out/$temp1.txt" ) or die; + + for(my $i=0; $i<=$#lines1; $i=$i+1) { + my $temp = $lines1[$i]; + $temp =~ m/\.\.\.\.\.\.([AGCT]{5})\s/i or die "##$temp##"; + my $motif = $1; + print fileA "$motif\n"; + } + + close(fileA); + system( "sort $out/$temp1.txt | uniq -c | sort -n -k 1 -r -o $out/$temp1.number" ); + system( "Rscript 7-histogram.R -i $out/$temp1.number -o $out/$temp1.pdf" ); +} + +################################################################################################################################################################################################### + + + + + + +say "Done ......"; + + + + + + +##### diff --git a/Call_Mutations/figures/7-histogram.R b/Call_Mutations/figures/7-histogram.R new file mode 100644 index 0000000..27eb0eb --- /dev/null +++ b/Call_Mutations/figures/7-histogram.R @@ -0,0 +1,60 @@ + +############################################################################################################################################################################################## +suppressPackageStartupMessages( library(optparse) ) ## To run the script in command lines. + +getParameters_f <- function() { + option_list_Local <- list( ## Options list with associated default value. + optparse::make_option(opt_str=c("-i", "--inputFile"), + default="7-getMotifs/Ig1.bed.number", + type="character", dest="inputFile", + help="Input file. [default: %default]."), + + optparse::make_option(opt_str=c("-o", "--outFile"), + default="7-getMotifs/Ig1.bed.pdf", + type="character", dest="outFile", + help="Output file. [default: %default].") + ) + + ## Now parse the command line to check which option is given and get associated values. + parser_Local <- optparse::OptionParser(usage="usage: Rscript %prog [options]", + option_list=option_list_Local, + description="May 22, 2021.", + epilogue="For comments, bug reports etc..., please contact Yong Peng ." + ) + opt_Local <- optparse::parse_args(parser_Local, args=commandArgs(trailingOnly=TRUE), positional_arguments=0)$options + return(opt_Local) +} +############################################################################################################################################################################################## + + + + + +############################################################################################################################################################################################## +opt_g = getParameters_f() + +inputFile_g <- as.character(opt_g$inputFile) +outFile_g <- as.character(opt_g$outFile) + + + +############################################# +matrix1 = read.table(file= inputFile_g, header = F, sep = character(0), quote = "\"'", dec = ".") +dim(matrix1) + + +ratio1 = as.numeric( matrix1[,1] ) +col_names1 = matrix1[,2] + + +pdf( outFile_g , height=5, width=10) + +xx1 <- barplot( ratio1, xaxt = 'n', xlab = 'motifs', ylab = "number", col="red", width = 1 ) +text(x = xx1, y = ratio1, label = ratio1, pos = 3, cex = 0.9, col = "black") +axis(1, at=xx1, labels=col_names1, tick=FALSE, las=2, line=0, cex.axis=0.6, cex=5 ) + +dev.off() + + + + diff --git a/Call_Mutations/figures/Anno_v20210106.R b/Call_Mutations/figures/Anno_v20210106.R new file mode 100644 index 0000000..184ab32 --- /dev/null +++ b/Call_Mutations/figures/Anno_v20210106.R @@ -0,0 +1,2180 @@ +############################################################################################################################################################################################## +## Author: Yong Peng, yongp@outlook.com. +## Run " Rscript Anno_v20210106.R -h" to get some help. +############################################################################################################################################################################################## + + + + + +############################################################################################################################################################################################## +suppressPackageStartupMessages( library(optparse) ) ## To run the script in command lines. + +getParameters_f <- function() { + option_list_Local <- list( ## Options list with associated default value. + optparse::make_option(opt_str=c("-i", "--inputDir"), + default="1-BED", + type="character", dest="inputDir", + help="Path to the directory containing BED files with genomic regions. The 1st to 6th columns of BED files, all of them are required. [default: %default]."), + + optparse::make_option(opt_str=c("-o", "--outDir"), + default="2-Annotation", + type="character", dest="outDir", + help="Path to the directory containing all the analysis results. [default: %default]."), + + optparse::make_option(opt_str=c("-r", "--refGenome"), + default="hg38", + type="character", dest="refGenome", + help="A string that defines the genome assembly or reference genome for your input BED files, such as hg38, hg19, mm10 and mm9. [default: %default]."), + + optparse::make_option(opt_str=c("-s", "--stranded"), + default=FALSE, + type="logical", dest="stranded", + help="Logical, whether find nearest/overlap gene in the same strand. [default: %default]."), + + optparse::make_option(opt_str=c("-f", "--fast"), + default=FALSE, + type="logical", dest="fast", + help="Logical, whether use fast mode. If it is TRUE, some steps will be skipped. [default: %default]."), + + optparse::make_option(opt_str=c("-U", "--upstream"), + default=3000, + type="integer", dest="upstream", + help="Number of upstream bases of TSS to define promoter. [default: %default]."), + + optparse::make_option(opt_str=c("-D", "--downstream"), + default=3000, + type="integer", dest="downstream", + help="Number of downstream bases of TSS to define promoter. [default: %default].") + + ) + + ## Now parse the command line to check which option is given and get associated values. + parser_Local <- optparse::OptionParser(usage="usage: Rscript %prog [options]", + option_list=option_list_Local, + description="ChIPseeker v6, February 22, 2020.", + epilogue="For comments, bug reports etc..., please contact Yong Peng ." + ) + opt_Local <- optparse::parse_args(parser_Local, args=commandArgs(trailingOnly=TRUE), positional_arguments=0)$options + return(opt_Local) +} +############################################################################################################################################################################################## + + + + + +############################################################################################################################################################################################## +opt_g = getParameters_f() + +inputDir_g <- as.character(opt_g$inputDir) +outDir_g <- as.character(opt_g$outDir) +refGenome_g <- as.character(opt_g$refGenome) +stranded_g <- as.logical(opt_g$stranded) +fast_g <- as.logical(opt_g$fast) +upstream_g <- as.integer(opt_g$upstream) +downstream_g <- as.integer(opt_g$downstream) + +rm(getParameters_f) +rm(opt_g) + +options(digits=10) +continue_on_error_g <- function() { + print( "NOTE: THERE WAS AN ERROR HERE. We are continuing because we have set 'options(error=continue_on_error())'. " ) +} +options( error=continue_on_error_g ) ## This option is very important. + +# inputDir_g <- "1-BED" +# outDir_g <- "2-Annotation" +# refGenome_g <- "hg38" +# stranded_g <- TRUE +# fast_g <- FALSE +# upstream_g <- 1000 +# downstream_g <- 0 + +if( ! file.exists(outDir_g) ) { dir.create(outDir_g, recursive = TRUE) } +###################################################################################################################################################### + + + + + +###################################################################################################################################################### +continue_on_error <- function() { + print(" NOTE: THERE WAS AN ERROR HERE. Yong Peng. We are continuing because we have set 'options(error=continue_on_error())'. ") +} +# This is the key option +options(error=continue_on_error) + +suppressPackageStartupMessages( library(GenomicFeatures) ) +suppressPackageStartupMessages( library(GenomicRanges) ) +suppressPackageStartupMessages( library(clusterProfiler) ) +suppressPackageStartupMessages( library(ReactomePA) ) +suppressPackageStartupMessages( library(ChIPseeker) ) +suppressPackageStartupMessages( library(DOSE) ) +suppressPackageStartupMessages( library(ggplot2) ) +suppressPackageStartupMessages( library(topGO) ) +suppressPackageStartupMessages( library(KEGG.db) ) +suppressPackageStartupMessages( library(enrichplot) ) +suppressPackageStartupMessages( library(ggupset) ) +suppressPackageStartupMessages( library(ggimage) ) + +my_txdb_g <- "" +my_orgdb_g <- "" +my_organism_g <- "" +my_organism2_g <- "" +if(refGenome_g == "hg38") { + suppressPackageStartupMessages( library(TxDb.Hsapiens.UCSC.hg38.knownGene) ) + suppressPackageStartupMessages( library(org.Hs.eg.db) ) + my_txdb_g <- TxDb.Hsapiens.UCSC.hg38.knownGene + my_orgdb_g <- "org.Hs.eg.db" + my_organism_g <- "human" + my_organism2_g <- "hsa" +} +if(refGenome_g == "mm10") { + suppressPackageStartupMessages( library(TxDb.Mmusculus.UCSC.mm10.knownGene) ) + suppressPackageStartupMessages( library(org.Mm.eg.db) ) + my_txdb_g <- TxDb.Mmusculus.UCSC.mm10.knownGene + my_orgdb_g <- "org.Mm.eg.db" + my_organism_g <- "mouse" + my_organism2_g <- "mmu" +} +if(refGenome_g == "danRer11") { + suppressPackageStartupMessages( library(TxDb.Drerio.UCSC.danRer11.refGene) ) + suppressPackageStartupMessages( library(org.Dr.eg.db) ) + my_txdb_g <- TxDb.Drerio.UCSC.danRer11.refGene + my_orgdb_g <- "org.Dr.eg.db" + my_organism_g <- "zebrafish" + my_organism2_g <- "dre" +} + +sink( paste(outDir_g, "Parameters.log.txt", sep="/") ) + print( "inputDir_g:" ) + print( inputDir_g ) + cat("\n\n\n") + print( "outDir_g: ") + print( outDir_g ) + cat("\n\n\n") + print( "refGenome_g: ") + print( refGenome_g ) + cat("\n\n\n") + print( "stranded_g: ") + print( stranded_g ) + cat("\n\n\n") + print( "fast_g: ") + print( fast_g ) + cat("\n\n\n") + print( "upstream_g: ") + print( upstream_g ) + cat("\n\n\n") + print( "downstream_g: ") + print( downstream_g ) + cat("\n\n\n") + print( "my_organism_g:" ) + print( my_organism_g ) + cat("\n\n\n") + print( "my_organism2_g:" ) + print( my_organism2_g ) + cat("\n\n\n") +sink() + + +## read the input files. +peakFiles <- list.files(path = inputDir_g, pattern = ".bed$", full.names = TRUE ) +peakFiles_onlyName <- list.files(path = inputDir_g, pattern = ".bed$", full.names = FALSE ) +cat("####################\n\n\n") +print(peakFiles_onlyName) +print(length(peakFiles_onlyName)) +cat("####################\n\n\n") + +sink( paste(outDir_g, "All_the_input_BED_files.txt", sep="/") ) +print(peakFiles) +print(length(peakFiles)) +cat("####################\n\n\n") +print(peakFiles_onlyName) +print(length(peakFiles_onlyName)) +sink() +###################################################################################################################################################### + + + + + +###################################################################################################################################################### +MyTheme_1 <- function(textSize1=14, hjust1=NULL, vjust1=NULL, angle1=NULL) { # "hjust=1, vjust=1, angle=30" for some boxplots. + theme( + line = element_line(colour="black", size=1.0, linetype=1, lineend=NULL), ## all line elements. 局部优先总体,下面3个也是,只对非局部设置有效. 所有线属性. + rect = element_rect(colour="black", size=1.0, linetype=1, fill="transparent" ), ## all rectangluar elements. hjust=1: 靠右对齐. 所有矩形区域属性. + text = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## all text elements. "serif" for a serif font. 所有文本相关属性. + title = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## all title elements: plot, axes, legends. hjust:水平对齐的方向. 所有标题属性. + ## aspect.ratio = 1, ##高宽比 + + axis.title = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## label of axes (element_text; inherits from text). horizontal: 水平的, 水平线 + axis.title.x = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## x axis label (element_text; inherits from axis.title) + axis.title.y = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=90, lineheight=1.0, margin = NULL, debug = NULL), ## y axis label (element_text; inherits from axis.title) + axis.text = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## tick labels along axes (element_text; inherits from text). 坐标轴刻度的标签的属性. + axis.text.x = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=hjust1, vjust=vjust1, angle=angle1, lineheight=1.0, margin = NULL, debug = NULL), ## x axis tick labels (element_text; inherits from axis.text) + axis.text.y = element_text(family="serif", face="plain", colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=0, lineheight=1.0, margin = NULL, debug = NULL), ## y axis tick labels (element_text; inherits from axis.text) + + axis.ticks = element_line(colour="black", size=0.5, linetype=1, lineend=NULL), ## tick marks along axes (element_line; inherits from line). 坐标轴刻度线. + axis.ticks.x = element_line(colour="black", size=0.5, linetype=1, lineend=NULL), ## x axis tick marks (element_line; inherits from axis.ticks) + axis.ticks.y = element_line(colour="black", size=0.5, linetype=1, lineend=NULL), ## y axis tick marks (element_line; inherits from axis.ticks) + axis.ticks.length = grid::unit(2.0, "mm", data=NULL), ## length of tick marks (unit), ‘"mm"’ Millimetres. 10 mm = 1 cm. 刻度线长度 + axis.line = element_line(colour="transparent", size=0.3, linetype=1, lineend=NULL), ## lines along axes (element_line; inherits from line). 坐标轴线 + axis.line.x = element_line(colour="transparent", size=0.3, linetype=1, lineend=NULL), ## line along x axis (element_line; inherits from axis.line) + axis.line.y = element_line(colour="transparent", size=0.3, linetype=1, lineend=NULL), ## line along y axis (element_line; inherits from axis.line) + + legend.background = element_rect(colour="transparent", size=1, linetype=1, fill="transparent" ), ## background of legend (element_rect; inherits from rect) + legend.spacing = grid::unit(1, "mm", data=NULL), ## extra space added around legend (unit). linetype=1指的是矩形边框的类型. + legend.key = element_rect(colour="transparent", size=2, linetype=1, fill="transparent" ), ## background underneath legend keys. 图例符号. size=1指的是矩形边框的大小. + legend.key.size = grid::unit(6, "mm", data=NULL) , ## size of legend keys (unit; inherits from legend.key.size) + legend.key.height = grid::unit(6.5, "mm", data=NULL) , ## key background height (unit; inherits from legend.key.size) + legend.key.width = grid::unit(8, "mm", data=NULL) , ## key background width (unit; inherits from legend.key.size) + legend.text = element_text(family="serif", face=NULL, colour="black", size=textSize1, hjust=NULL, vjust=NULL, angle=NULL, lineheight=NULL), ##legend item labels. 图例文字标签. + legend.text.align = 0, ## alignment of legend labels (number from 0 (left) to 1 (right)) + legend.title = element_blank(), ## title of legend (element_text; inherits from title) + legend.title.align = 0, ## alignment of legend title (number from 0 (left) to 1 (right)) + legend.position = "right", ## the position of legends. ("left", "right", "bottom", "top", or two-element numeric vector) + legend.direction = "vertical", ## layout of items in legends ("horizontal" or "vertical") 图例排列方向 + legend.justification = "center", ## anchor point for positioning legend inside plot ("center" or two-element numeric vector) 图例居中方式 + legend.box = NULL, ## arrangement of multiple legends ("horizontal" or "vertical") 多图例的排列方式 + legend.box.just = NULL, ## justification of each legend within the overall bounding box, when there are multiple legends ("top", "bottom", "left", or "right") 多图例的居中方式 + + panel.background = element_rect(colour="transparent", size=0.0, linetype=1, fill="transparent" ), ## background of plotting area, drawn underneath plot (element_rect; inherits from rect) + panel.border = element_rect(colour="black", size=0.5, linetype=1, fill=NA ), ## border around plotting area, drawn on top of plot so that it covers tick marks and grid lines. This should be used with fill=NA (element_rect; inherits from rect) + panel.spacing = grid::unit(1, "mm", data=NULL) , ## margin around facet panels (unit) 分面绘图区之间的边距 + panel.spacing.x = grid::unit(1, "mm", data=NULL) , + panel.spacing.y = grid::unit(1, "mm", data=NULL) , + panel.grid = element_blank(), ## grid lines (element_line; inherits from line) 绘图区网格线 + panel.grid.major = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## major grid lines (element_line; inherits from panel.grid) 主网格线 + panel.grid.minor = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## minor grid lines (element_line; inherits from panel.grid) 次网格线 + panel.grid.major.x = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## vertical major grid lines (element_line; inherits from panel.grid.major) + panel.grid.major.y = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## horizontal major grid lines (element_line; inherits from panel.grid.major) + panel.grid.minor.x = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## vertical minor grid lines (element_line; inherits from panel.grid.minor) + panel.grid.minor.y = element_line(colour="transparent", size=NULL, linetype=NULL, lineend=NULL) , ## horizontal minor grid lines (element_line; inherits from panel.grid.minor) + + plot.background = element_rect(colour="transparent", size=NULL, linetype=NULL, fill="transparent" ), ## background of the entire plot (element_rect; inherits from rect) 整个图形的背景 + plot.title = element_text(family="serif", face=NULL, colour="black", size=textSize1, hjust=0.5, vjust=0.5, angle=NULL, lineheight=NULL), ## plot title (text appearance) (element_text; inherits from title) 图形标题 + plot.margin = grid::unit(c(5, 5, 5, 5), "mm", data=NULL), ## margin around entire plot (unit with the sizes of the top, right, bottom, and left margins) + + strip.background = element_rect(colour=NULL, size=NULL, linetype=NULL, fill=NULL ), ## background of facet labels (element_rect; inherits from rect) 分面标签背景 + strip.text = element_text(family="serif", face=NULL, colour=NULL, size=NULL, hjust=NULL, vjust=NULL, angle=NULL, lineheight=NULL), ## facet labels (element_text; inherits from text) + strip.text.x = element_text(family="serif", face=NULL, colour=NULL, size=NULL, hjust=NULL, vjust=NULL, angle=NULL, lineheight=NULL), ## facet labels along horizontal direction (element_text; inherits from strip.text) + strip.text.y = element_text(family="serif", face=NULL, colour=NULL, size=NULL, hjust=NULL, vjust=NULL, angle=NULL, lineheight=NULL) ## facet labels along vertical direction (element_text; inherits from strip.text) + ) +} + + +MySaveGgplot2_1 <- function(path1, fileName1, height1, width1) { + SVG1 <- paste(path1, "/", "SVG", sep = "", collapse = NULL) + PNG1 <- paste(path1, "/", "PNG", sep = "", collapse = NULL) + PDF1 <- paste(path1, "/", "PDF", sep = "", collapse = NULL) + EPS1 <- paste(path1, "/", "EPS", sep = "", collapse = NULL) + if( ! file.exists(SVG1) ) { dir.create(SVG1) } + if( ! file.exists(PNG1) ) { dir.create(PNG1) } + if( ! file.exists(PDF1) ) { dir.create(PDF1) } + if( ! file.exists(EPS1) ) { dir.create(EPS1) } + ggsave( filename = paste(SVG1, "/", fileName1, ".svg", sep="", collapse=NULL), height=height1, width=width1, dpi = 1200 ) + ggsave( filename = paste(PNG1, "/", fileName1, ".png", sep="", collapse=NULL), height=height1, width=width1, dpi = 1200 ) + ggsave( filename = paste(PDF1, "/", fileName1, ".pdf", sep="", collapse=NULL), height=height1, width=width1, dpi = 1200 ) + ggsave( filename = paste(EPS1, "/", fileName1, ".eps", sep="", collapse=NULL), height=height1, width=width1, dpi = 1200, device=cairo_ps) +} +###################################################################################################################################################### + + + + + +###################################################################################################################################################### +############# For annotating one file or one group. + + +## +MyPeaksAnno_OneGroup_1_g <- function(myPeak_t1, myPath_t1) { + myTempFunction1 <- function() { + if( ! file.exists(myPath_t1) ) { dir.create(myPath_t1, recursive = TRUE) } + pdf(file=paste(myPath_t1, "1_GenomicRegions_over_Chromosomes.pdf", sep="/"), width=20, height=10) + print( covplot(peak=myPeak_t1, weightCol = "V5", xlab = "Chromosome Size (bp)", ylab = "strength", title = "Genomic Regions over Chromosomes", chrs = NULL, xlim = NULL, lower = 0.01) ) + dev.off() + pdf(file=paste(myPath_t1, "2_GenomicRegions_over_only-chromXY.pdf", sep="/"), width=20, height=5) + print( covplot(peak=myPeak_t1, weightCol = "V5", xlab = "Chromosome Size (bp)", ylab = "strength", title = "Genomic Regions over Chromosome X and Y", chrs = c("chrX", "chrY"), xlim = NULL, lower = 0.01) ) + dev.off() + pdf(file=paste(myPath_t1, "3_chromLength.pdf", sep="/"), width=20, height=10) + print( covplot(peak=myPeak_t1, weightCol = "V5", xlab = "Chromosome Size (bp)", ylab = "strength", title = "Length of Each Chromosome", chrs = NULL, xlim = NULL, lower = -1) ) + dev.off() + } + tryCatch( + myTempFunction1(), + error = function(err){"MyPeaksAnno_OneGroup_1_g:error_1. Yong Peng."} + ) +} + + +## +MyPeaksSignal_OneGroup_1_g <- function(myPeak_t1, myPath_t1, up1, down1) { + myTempFunction2 <- function() { + if( ! file.exists(myPath_t1) ) { dir.create(myPath_t1, recursive = TRUE) } + promoter2 <- getPromoters(TxDb=my_txdb_g, upstream=up1, downstream=down1) + tagMatrix1 <- getTagMatrix(myPeak_t1, windows=promoter2) + write.table(x=tagMatrix1, file =paste(myPath_t1, "1_peakMatrix_TSS.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = TRUE, col.names = TRUE ) + + pdf(file = paste(myPath_t1, "2_Heatmap_of_GenomicRegions_binding_to_TSS_regions.pdf", sep="/"), width=12, height=10) + print( tagHeatmap(tagMatrix1, xlim=c(-up1, down1), color="red") ) + dev.off() + + pdf(file=paste(myPath_t1, "3_Average_Profile_of_GenomicRegions_binding_to_TSS_region.pdf", sep="/"), width=12, height=10) + print( plotAvgProf(tagMatrix1, xlim=c(-up1, down1), xlab="Genomic Region (5'->3')", ylab = "Peak Count Frequency" ) ) + dev.off() + } + tryCatch( + myTempFunction2(), + error = function(err){"MyPeaksSignal_OneGroup_1_g:error_2."} + ) +} + + +## +MyPeaksDistribution_OneGroup_1_g <- function(myPeak_t1, myPath_t1, up1, down1) { + myTempFunction3 <- function() { + if( ! file.exists(myPath_t1) ) { dir.create(myPath_t1, recursive = TRUE) } + peakAnno <- annotatePeak( myPeak_t1, tssRegion=c(-up1, down1), TxDb=my_txdb_g, annoDb=my_orgdb_g, sameStrand = stranded_g , overlap="all" ) + sink( file = paste(myPath_t1, "1A_peakAnnotation_summary.txt", sep="/") ) + print( peakAnno ) + sink() + write.table(x=as.data.frame(peakAnno), file = paste(myPath_t1, "1B_peakAnnotation_details.txt", sep="/"), append = FALSE, + quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + pdf(file=paste(myPath_t1, "2A_GenomicDistribution_by_pieplot.pdf", sep="/"), width=7, height=5) + print( plotAnnoPie(peakAnno) ) + dev.off() + svg(filename=paste(myPath_t1, "2B_GenomicDistribution_by_pieplot.svg", sep="/"), width=7, height=5) + print( plotAnnoPie(peakAnno) ) + dev.off() + + pdf(file=paste(myPath_t1, "3A_GenomicDistribution_by_barplot.pdf", sep="/"), width=7, height=3) + print( plotAnnoBar(peakAnno) ) + dev.off() + svg(filename=paste(myPath_t1, "3B_GenomicDistribution_by_barplot.svg", sep="/"), width=7, height=3) + print( plotAnnoBar(peakAnno) ) + dev.off() + + pdf(file=paste(myPath_t1, "4A_GenomicDistribution_by_vennpie.pdf", sep="/"), width=7, height=5) + print( vennpie(peakAnno) ) + dev.off() + svg(filename=paste(myPath_t1, "4B_GenomicDistribution_by_vennpie.svg", sep="/"), width=7, height=5) + print( vennpie(peakAnno) ) + dev.off() + + pdf(file=paste(myPath_t1, "5A_GenomicDistribution_by_upsetplot.pdf", sep="/"), width=8, height=5) + print( upsetplot(peakAnno) ) + dev.off() + svg(filename=paste(myPath_t1, "5B_GenomicDistribution_by_upsetplot.svg", sep="/"), width=8, height=5) + print( upsetplot(peakAnno) ) + dev.off() + + pdf(file=paste(myPath_t1, "6A_GenomicDistribution_by_upsetplot_vennpie.pdf", sep="/"), width=8, height=5) + print( upsetplot(peakAnno, vennpie=TRUE) ) + dev.off() + svg(filename=paste(myPath_t1, "6B_GenomicDistribution_by_upsetplot_vennpie.svg", sep="/"), width=8, height=5) + print( upsetplot(peakAnno, vennpie=TRUE) ) + dev.off() + + pdf(file=paste(myPath_t1, "7A_Distribution_of_peaks_to_TSS.pdf", sep="/"), width=8, height=2) + print( plotDistToTSS(peakAnno, title="Distribution of peaks relative to TSS") ) + dev.off() + svg(filename=paste(myPath_t1, "7B_Distribution_of_peaks_to_TSS.svg", sep="/"), width=800, height=200) + print( plotDistToTSS(peakAnno, title="Distribution of peaks relative to TSS") ) + dev.off() + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + +############# +MyGeneOntology_BP_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- clusterProfiler::enrichGO(gene=as.data.frame(peakAnno2)$geneId, OrgDb=my_orgdb_g, keyType = "ENTREZID", ont = "BP", + pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.2, readable = TRUE, pool = FALSE) + + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 10) ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + # simplify + enrich_raw_simplify <- clusterProfiler::simplify(enrich_raw, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_raw_simplify , file = paste(MyFolder, "/", MyFileName, "4_enrich_raw_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_5_enrich_raw_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_6_enrich_raw_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 30) ) + dev.off() + enrich_raw_simplify_dataFrame <- as.data.frame( enrich_raw_simplify ) + if(nrow(enrich_raw_simplify_dataFrame) > 1) { + enrich_raw_simplify_dataFrame$GeneRatio <- sapply(enrich_raw_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_simplify_dataFrame$pvalue <- -log10(enrich_raw_simplify_dataFrame$pvalue) + enrich_raw_simplify_dataFrame$p.adjust <- -log10(enrich_raw_simplify_dataFrame$p.adjust) + enrich_raw_simplify_dataFrame$qvalue <- -log10(enrich_raw_simplify_dataFrame$qvalue) + enrich_raw_simplify_dataFrame_selected <- enrich_raw_simplify_dataFrame[1:20, -8] + if( nrow(enrich_raw_simplify_dataFrame) <20 ) {enrich_raw_simplify_dataFrame_selected <- enrich_raw_simplify_dataFrame[, -8]} + enrich_raw_simplify_dataFrame_selected <- transform(enrich_raw_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_simplify_dataFrame_selected$GeneRatio) + min(enrich_raw_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_simplify_dataFrame_selected$ID, labels=rev(enrich_raw_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_6_enrich_raw_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + + + ######################### + enrich_drop1 <- dropGO(enrich_raw, level = c(1:4), term = NULL) ## drop specific GO terms or level + + write.table(x= enrich_drop1 , file = paste(MyFolder, "/", MyFileName, "7_enrich_drop1.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_8_enrich_drop1.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop1, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_9_enrich_drop1.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1, showCategory = 30) ) + dev.off() + enrich_drop1_dataFrame <- as.data.frame( enrich_drop1 ) + if(nrow(enrich_drop1_dataFrame) > 1) { + enrich_drop1_dataFrame$GeneRatio <- sapply(enrich_drop1$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop1_dataFrame$pvalue <- -log10(enrich_drop1_dataFrame$pvalue) + enrich_drop1_dataFrame$p.adjust <- -log10(enrich_drop1_dataFrame$p.adjust) + enrich_drop1_dataFrame$qvalue <- -log10(enrich_drop1_dataFrame$qvalue) + enrich_drop1_dataFrame_selected <- enrich_drop1_dataFrame[1:20, -8] + if( nrow(enrich_drop1_dataFrame) <20 ) {enrich_drop1_dataFrame_selected <- enrich_drop1_dataFrame[, -8]} + enrich_drop1_dataFrame_selected <- transform(enrich_drop1_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop1_dataFrame_selected$GeneRatio) + min(enrich_drop1_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop1_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop1_dataFrame_selected$ID, labels=rev(enrich_drop1_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_9_enrich_drop1", sep=""), height1=6, width1=5) + } + + # simplify + enrich_drop1_simplify <- clusterProfiler::simplify(enrich_drop1, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_drop1_simplify , file = paste(MyFolder, "/", MyFileName, "10_enrich_drop1_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_11_enrich_drop1_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_12_enrich_drop1_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 30) ) + dev.off() + enrich_drop1_simplify_dataFrame <- as.data.frame( enrich_drop1_simplify ) + if(nrow(enrich_drop1_simplify_dataFrame) > 1) { + enrich_drop1_simplify_dataFrame$GeneRatio <- sapply(enrich_drop1_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop1_simplify_dataFrame$pvalue <- -log10(enrich_drop1_simplify_dataFrame$pvalue) + enrich_drop1_simplify_dataFrame$p.adjust <- -log10(enrich_drop1_simplify_dataFrame$p.adjust) + enrich_drop1_simplify_dataFrame$qvalue <- -log10(enrich_drop1_simplify_dataFrame$qvalue) + enrich_drop1_simplify_dataFrame_selected <- enrich_drop1_simplify_dataFrame[1:20, -8] + if( nrow(enrich_drop1_simplify_dataFrame) <20 ) {enrich_drop1_simplify_dataFrame_selected <- enrich_drop1_simplify_dataFrame[, -8]} + enrich_drop1_simplify_dataFrame_selected <- transform(enrich_drop1_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop1_simplify_dataFrame_selected$GeneRatio) + min(enrich_drop1_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop1_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop1_simplify_dataFrame_selected$ID, labels=rev(enrich_drop1_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_12_enrich_drop1_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + + + + ######################### + enrich_drop2 <- dropGO(enrich_raw, level = c(1:5), term = NULL) ## drop specific GO terms or level + + write.table(x= enrich_drop2 , file = paste(MyFolder, "/", MyFileName, "13_enrich_drop2.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_14_enrich_drop2.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop2, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_15_enrich_drop2.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2, showCategory = 30) ) + dev.off() + enrich_drop2_dataFrame <- as.data.frame( enrich_drop2 ) + if(nrow(enrich_drop2_dataFrame) > 1) { + enrich_drop2_dataFrame$GeneRatio <- sapply(enrich_drop2$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop2_dataFrame$pvalue <- -log10(enrich_drop2_dataFrame$pvalue) + enrich_drop2_dataFrame$p.adjust <- -log10(enrich_drop2_dataFrame$p.adjust) + enrich_drop2_dataFrame$qvalue <- -log10(enrich_drop2_dataFrame$qvalue) + enrich_drop2_dataFrame_selected <- enrich_drop2_dataFrame[1:20, -8] + if( nrow(enrich_drop2_dataFrame) <20 ) {enrich_drop2_dataFrame_selected <- enrich_drop2_dataFrame[, -8]} + enrich_drop2_dataFrame_selected <- transform(enrich_drop2_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop2_dataFrame_selected$GeneRatio) + min(enrich_drop2_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop2_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop2_dataFrame_selected$ID, labels=rev(enrich_drop2_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_15_enrich_drop2", sep=""), height1=6, width1=5) + } + + # simplify + enrich_drop2_simplify <- clusterProfiler::simplify(enrich_drop2, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_drop2_simplify , file = paste(MyFolder, "/", MyFileName, "16_enrich_drop2_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_17_enrich_drop2_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_18_enrich_drop2_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 30) ) + dev.off() + enrich_drop2_simplify_dataFrame <- as.data.frame( enrich_drop2_simplify ) + if(nrow(enrich_drop2_simplify_dataFrame) > 1) { + enrich_drop2_simplify_dataFrame$GeneRatio <- sapply(enrich_drop2_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop2_simplify_dataFrame$pvalue <- -log10(enrich_drop2_simplify_dataFrame$pvalue) + enrich_drop2_simplify_dataFrame$p.adjust <- -log10(enrich_drop2_simplify_dataFrame$p.adjust) + enrich_drop2_simplify_dataFrame$qvalue <- -log10(enrich_drop2_simplify_dataFrame$qvalue) + enrich_drop2_simplify_dataFrame_selected <- enrich_drop2_simplify_dataFrame[1:20, -8] + if( nrow(enrich_drop2_simplify_dataFrame) <20 ) {enrich_drop2_simplify_dataFrame_selected <- enrich_drop2_simplify_dataFrame[, -8]} + enrich_drop2_simplify_dataFrame_selected <- transform(enrich_drop2_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop2_simplify_dataFrame_selected$GeneRatio) + min(enrich_drop2_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop2_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop2_simplify_dataFrame_selected$ID, labels=rev(enrich_drop2_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_18_enrich_drop2_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + +############# +MyGeneOntology_MF_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- clusterProfiler::enrichGO(gene=as.data.frame(peakAnno2)$geneId, OrgDb=my_orgdb_g, keyType = "ENTREZID", ont = "MF", + pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.2, readable = TRUE, pool = FALSE) + + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 10) ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + # simplify + enrich_raw_simplify <- clusterProfiler::simplify(enrich_raw, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_raw_simplify , file = paste(MyFolder, "/", MyFileName, "4_enrich_raw_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_5_enrich_raw_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_6_enrich_raw_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 30) ) + dev.off() + enrich_raw_simplify_dataFrame <- as.data.frame( enrich_raw_simplify ) + if(nrow(enrich_raw_simplify_dataFrame) > 1) { + enrich_raw_simplify_dataFrame$GeneRatio <- sapply(enrich_raw_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_simplify_dataFrame$pvalue <- -log10(enrich_raw_simplify_dataFrame$pvalue) + enrich_raw_simplify_dataFrame$p.adjust <- -log10(enrich_raw_simplify_dataFrame$p.adjust) + enrich_raw_simplify_dataFrame$qvalue <- -log10(enrich_raw_simplify_dataFrame$qvalue) + enrich_raw_simplify_dataFrame_selected <- enrich_raw_simplify_dataFrame[1:20, -8] + if( nrow(enrich_raw_simplify_dataFrame) <20 ) {enrich_raw_simplify_dataFrame_selected <- enrich_raw_simplify_dataFrame[, -8]} + enrich_raw_simplify_dataFrame_selected <- transform(enrich_raw_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_simplify_dataFrame_selected$GeneRatio) + min(enrich_raw_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_simplify_dataFrame_selected$ID, labels=rev(enrich_raw_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_6_enrich_raw_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + + + ######################### + enrich_drop1 <- dropGO(enrich_raw, level = c(1:4), term = NULL) ## drop specific GO terms or level + + write.table(x= enrich_drop1 , file = paste(MyFolder, "/", MyFileName, "7_enrich_drop1.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_8_enrich_drop1.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop1, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_9_enrich_drop1.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1, showCategory = 30) ) + dev.off() + enrich_drop1_dataFrame <- as.data.frame( enrich_drop1 ) + if(nrow(enrich_drop1_dataFrame) > 1) { + enrich_drop1_dataFrame$GeneRatio <- sapply(enrich_drop1$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop1_dataFrame$pvalue <- -log10(enrich_drop1_dataFrame$pvalue) + enrich_drop1_dataFrame$p.adjust <- -log10(enrich_drop1_dataFrame$p.adjust) + enrich_drop1_dataFrame$qvalue <- -log10(enrich_drop1_dataFrame$qvalue) + enrich_drop1_dataFrame_selected <- enrich_drop1_dataFrame[1:20, -8] + if( nrow(enrich_drop1_dataFrame) <20 ) {enrich_drop1_dataFrame_selected <- enrich_drop1_dataFrame[, -8]} + enrich_drop1_dataFrame_selected <- transform(enrich_drop1_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop1_dataFrame_selected$GeneRatio) + min(enrich_drop1_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop1_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop1_dataFrame_selected$ID, labels=rev(enrich_drop1_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_9_enrich_drop1", sep=""), height1=6, width1=5) + } + + # simplify + enrich_drop1_simplify <- clusterProfiler::simplify(enrich_drop1, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_drop1_simplify , file = paste(MyFolder, "/", MyFileName, "10_enrich_drop1_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_11_enrich_drop1_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_12_enrich_drop1_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 30) ) + dev.off() + enrich_drop1_simplify_dataFrame <- as.data.frame( enrich_drop1_simplify ) + if(nrow(enrich_drop1_simplify_dataFrame) > 1) { + enrich_drop1_simplify_dataFrame$GeneRatio <- sapply(enrich_drop1_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop1_simplify_dataFrame$pvalue <- -log10(enrich_drop1_simplify_dataFrame$pvalue) + enrich_drop1_simplify_dataFrame$p.adjust <- -log10(enrich_drop1_simplify_dataFrame$p.adjust) + enrich_drop1_simplify_dataFrame$qvalue <- -log10(enrich_drop1_simplify_dataFrame$qvalue) + enrich_drop1_simplify_dataFrame_selected <- enrich_drop1_simplify_dataFrame[1:20, -8] + if( nrow(enrich_drop1_simplify_dataFrame) <20 ) {enrich_drop1_simplify_dataFrame_selected <- enrich_drop1_simplify_dataFrame[, -8]} + enrich_drop1_simplify_dataFrame_selected <- transform(enrich_drop1_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop1_simplify_dataFrame_selected$GeneRatio) + min(enrich_drop1_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop1_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop1_simplify_dataFrame_selected$ID, labels=rev(enrich_drop1_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_12_enrich_drop1_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + + + + ######################### + enrich_drop2 <- dropGO(enrich_raw, level = c(1:5), term = NULL) ## drop specific GO terms or level + + write.table(x= enrich_drop2 , file = paste(MyFolder, "/", MyFileName, "13_enrich_drop2.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_14_enrich_drop2.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop2, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_15_enrich_drop2.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2, showCategory = 30) ) + dev.off() + enrich_drop2_dataFrame <- as.data.frame( enrich_drop2 ) + if(nrow(enrich_drop2_dataFrame) > 1) { + enrich_drop2_dataFrame$GeneRatio <- sapply(enrich_drop2$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop2_dataFrame$pvalue <- -log10(enrich_drop2_dataFrame$pvalue) + enrich_drop2_dataFrame$p.adjust <- -log10(enrich_drop2_dataFrame$p.adjust) + enrich_drop2_dataFrame$qvalue <- -log10(enrich_drop2_dataFrame$qvalue) + enrich_drop2_dataFrame_selected <- enrich_drop2_dataFrame[1:20, -8] + if( nrow(enrich_drop2_dataFrame) <20 ) {enrich_drop2_dataFrame_selected <- enrich_drop2_dataFrame[, -8]} + enrich_drop2_dataFrame_selected <- transform(enrich_drop2_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop2_dataFrame_selected$GeneRatio) + min(enrich_drop2_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop2_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop2_dataFrame_selected$ID, labels=rev(enrich_drop2_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_15_enrich_drop2", sep=""), height1=6, width1=5) + } + + # simplify + enrich_drop2_simplify <- clusterProfiler::simplify(enrich_drop2, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_drop2_simplify , file = paste(MyFolder, "/", MyFileName, "16_enrich_drop2_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_17_enrich_drop2_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_18_enrich_drop2_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 30) ) + dev.off() + enrich_drop2_simplify_dataFrame <- as.data.frame( enrich_drop2_simplify ) + if(nrow(enrich_drop2_simplify_dataFrame) > 1) { + enrich_drop2_simplify_dataFrame$GeneRatio <- sapply(enrich_drop2_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop2_simplify_dataFrame$pvalue <- -log10(enrich_drop2_simplify_dataFrame$pvalue) + enrich_drop2_simplify_dataFrame$p.adjust <- -log10(enrich_drop2_simplify_dataFrame$p.adjust) + enrich_drop2_simplify_dataFrame$qvalue <- -log10(enrich_drop2_simplify_dataFrame$qvalue) + enrich_drop2_simplify_dataFrame_selected <- enrich_drop2_simplify_dataFrame[1:20, -8] + if( nrow(enrich_drop2_simplify_dataFrame) <20 ) {enrich_drop2_simplify_dataFrame_selected <- enrich_drop2_simplify_dataFrame[, -8]} + enrich_drop2_simplify_dataFrame_selected <- transform(enrich_drop2_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop2_simplify_dataFrame_selected$GeneRatio) + min(enrich_drop2_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop2_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop2_simplify_dataFrame_selected$ID, labels=rev(enrich_drop2_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_18_enrich_drop2_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + +############# +MyGeneOntology_CC_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- clusterProfiler::enrichGO(gene=as.data.frame(peakAnno2)$geneId, OrgDb=my_orgdb_g, keyType = "ENTREZID", ont = "CC", + pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.2, readable = TRUE, pool = FALSE) + + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 10) ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + # simplify + enrich_raw_simplify <- clusterProfiler::simplify(enrich_raw, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_raw_simplify , file = paste(MyFolder, "/", MyFileName, "4_enrich_raw_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_5_enrich_raw_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_raw_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_6_enrich_raw_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw_simplify, showCategory = 30) ) + dev.off() + enrich_raw_simplify_dataFrame <- as.data.frame( enrich_raw_simplify ) + if(nrow(enrich_raw_simplify_dataFrame) > 1) { + enrich_raw_simplify_dataFrame$GeneRatio <- sapply(enrich_raw_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_simplify_dataFrame$pvalue <- -log10(enrich_raw_simplify_dataFrame$pvalue) + enrich_raw_simplify_dataFrame$p.adjust <- -log10(enrich_raw_simplify_dataFrame$p.adjust) + enrich_raw_simplify_dataFrame$qvalue <- -log10(enrich_raw_simplify_dataFrame$qvalue) + enrich_raw_simplify_dataFrame_selected <- enrich_raw_simplify_dataFrame[1:20, -8] + if( nrow(enrich_raw_simplify_dataFrame) <20 ) {enrich_raw_simplify_dataFrame_selected <- enrich_raw_simplify_dataFrame[, -8]} + enrich_raw_simplify_dataFrame_selected <- transform(enrich_raw_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_simplify_dataFrame_selected$GeneRatio) + min(enrich_raw_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_simplify_dataFrame_selected$ID, labels=rev(enrich_raw_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_6_enrich_raw_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + + + ######################### + enrich_drop1 <- dropGO(enrich_raw, level = c(1:4), term = NULL) ## drop specific GO terms or level + + write.table(x= enrich_drop1 , file = paste(MyFolder, "/", MyFileName, "7_enrich_drop1.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_8_enrich_drop1.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop1, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop1, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop1, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_9_enrich_drop1.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1, showCategory = 30) ) + dev.off() + enrich_drop1_dataFrame <- as.data.frame( enrich_drop1 ) + if(nrow(enrich_drop1_dataFrame) > 1) { + enrich_drop1_dataFrame$GeneRatio <- sapply(enrich_drop1$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop1_dataFrame$pvalue <- -log10(enrich_drop1_dataFrame$pvalue) + enrich_drop1_dataFrame$p.adjust <- -log10(enrich_drop1_dataFrame$p.adjust) + enrich_drop1_dataFrame$qvalue <- -log10(enrich_drop1_dataFrame$qvalue) + enrich_drop1_dataFrame_selected <- enrich_drop1_dataFrame[1:20, -8] + if( nrow(enrich_drop1_dataFrame) <20 ) {enrich_drop1_dataFrame_selected <- enrich_drop1_dataFrame[, -8]} + enrich_drop1_dataFrame_selected <- transform(enrich_drop1_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop1_dataFrame_selected$GeneRatio) + min(enrich_drop1_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop1_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop1_dataFrame_selected$ID, labels=rev(enrich_drop1_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_9_enrich_drop1", sep=""), height1=6, width1=5) + } + + # simplify + enrich_drop1_simplify <- clusterProfiler::simplify(enrich_drop1, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_drop1_simplify , file = paste(MyFolder, "/", MyFileName, "10_enrich_drop1_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_11_enrich_drop1_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop1_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop1_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop1_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_12_enrich_drop1_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop1_simplify, showCategory = 30) ) + dev.off() + enrich_drop1_simplify_dataFrame <- as.data.frame( enrich_drop1_simplify ) + if(nrow(enrich_drop1_simplify_dataFrame) > 1) { + enrich_drop1_simplify_dataFrame$GeneRatio <- sapply(enrich_drop1_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop1_simplify_dataFrame$pvalue <- -log10(enrich_drop1_simplify_dataFrame$pvalue) + enrich_drop1_simplify_dataFrame$p.adjust <- -log10(enrich_drop1_simplify_dataFrame$p.adjust) + enrich_drop1_simplify_dataFrame$qvalue <- -log10(enrich_drop1_simplify_dataFrame$qvalue) + enrich_drop1_simplify_dataFrame_selected <- enrich_drop1_simplify_dataFrame[1:20, -8] + if( nrow(enrich_drop1_simplify_dataFrame) <20 ) {enrich_drop1_simplify_dataFrame_selected <- enrich_drop1_simplify_dataFrame[, -8]} + enrich_drop1_simplify_dataFrame_selected <- transform(enrich_drop1_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop1_simplify_dataFrame_selected$GeneRatio) + min(enrich_drop1_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop1_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop1_simplify_dataFrame_selected$ID, labels=rev(enrich_drop1_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_12_enrich_drop1_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + + + + ######################### + enrich_drop2 <- dropGO(enrich_raw, level = c(1:5), term = NULL) ## drop specific GO terms or level + + write.table(x= enrich_drop2 , file = paste(MyFolder, "/", MyFileName, "13_enrich_drop2.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_14_enrich_drop2.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop2, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop2, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop2, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_15_enrich_drop2.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2, showCategory = 30) ) + dev.off() + enrich_drop2_dataFrame <- as.data.frame( enrich_drop2 ) + if(nrow(enrich_drop2_dataFrame) > 1) { + enrich_drop2_dataFrame$GeneRatio <- sapply(enrich_drop2$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop2_dataFrame$pvalue <- -log10(enrich_drop2_dataFrame$pvalue) + enrich_drop2_dataFrame$p.adjust <- -log10(enrich_drop2_dataFrame$p.adjust) + enrich_drop2_dataFrame$qvalue <- -log10(enrich_drop2_dataFrame$qvalue) + enrich_drop2_dataFrame_selected <- enrich_drop2_dataFrame[1:20, -8] + if( nrow(enrich_drop2_dataFrame) <20 ) {enrich_drop2_dataFrame_selected <- enrich_drop2_dataFrame[, -8]} + enrich_drop2_dataFrame_selected <- transform(enrich_drop2_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop2_dataFrame_selected$GeneRatio) + min(enrich_drop2_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop2_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop2_dataFrame_selected$ID, labels=rev(enrich_drop2_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_15_enrich_drop2", sep=""), height1=6, width1=5) + } + + # simplify + enrich_drop2_simplify <- clusterProfiler::simplify(enrich_drop2, cutoff=0.7, by="p.adjust", select_fun=min) + write.table(x= enrich_drop2_simplify , file = paste(MyFolder, "/", MyFileName, "16_enrich_drop2_simplify.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_17_enrich_drop2_simplify.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 20 ) ) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 10) ) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 5) ) + print( enrichplot::emapplot(enrich_drop2_simplify, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_drop2_simplify, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::goplot(enrich_drop2_simplify, showCategory = 5) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_18_enrich_drop2_simplify.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_drop2_simplify, showCategory = 30) ) + dev.off() + enrich_drop2_simplify_dataFrame <- as.data.frame( enrich_drop2_simplify ) + if(nrow(enrich_drop2_simplify_dataFrame) > 1) { + enrich_drop2_simplify_dataFrame$GeneRatio <- sapply(enrich_drop2_simplify$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_drop2_simplify_dataFrame$pvalue <- -log10(enrich_drop2_simplify_dataFrame$pvalue) + enrich_drop2_simplify_dataFrame$p.adjust <- -log10(enrich_drop2_simplify_dataFrame$p.adjust) + enrich_drop2_simplify_dataFrame$qvalue <- -log10(enrich_drop2_simplify_dataFrame$qvalue) + enrich_drop2_simplify_dataFrame_selected <- enrich_drop2_simplify_dataFrame[1:20, -8] + if( nrow(enrich_drop2_simplify_dataFrame) <20 ) {enrich_drop2_simplify_dataFrame_selected <- enrich_drop2_simplify_dataFrame[, -8]} + enrich_drop2_simplify_dataFrame_selected <- transform(enrich_drop2_simplify_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_drop2_simplify_dataFrame_selected$GeneRatio) + min(enrich_drop2_simplify_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_drop2_simplify_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_drop2_simplify_dataFrame_selected$ID, labels=rev(enrich_drop2_simplify_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_18_enrich_drop2_simplify", sep=""), height1=6, width1=5) + } + + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + + +############# +MyReactome_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- enrichPathway(gene=as.data.frame(peakAnno2)$geneId, organism = my_organism_g, pvalueCutoff = 0.05, + pAdjustMethod = "BH", qvalueCutoff = 0.2, readable = TRUE) + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + +############# +MyKEGG_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- enrichKEGG(gene=as.data.frame(peakAnno2)$geneId, organism = my_organism2_g, keyType = "kegg", pvalueCutoff = 0.05, + pAdjustMethod = "BH", qvalueCutoff = 0.2, use_internal_data = TRUE ) + + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "Count", color = "p.adjust", font.size = 12, title = "") ) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "GeneRatio", color = "p.adjust", font.size = 12, title = "") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + + +############# +MyModuleKEGG_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- enrichMKEGG(gene=as.data.frame(peakAnno2)$geneId, organism = my_organism2_g, keyType = "kegg", pvalueCutoff = 0.05, + pAdjustMethod = "BH", qvalueCutoff = 0.2 ) + + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "Count", color = "p.adjust", font.size = 12, title = "") ) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "GeneRatio", color = "p.adjust", font.size = 12, title = "") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + +############# +MyEnrichDGN_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { ## Enrichment analysis based on the DisGeNET (http://www.disgenet.org/) + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- enrichDGN(gene=as.data.frame(peakAnno2)$geneId, pvalueCutoff = 0.05, pAdjustMethod = "BH", + qvalueCutoff = 0.2, readable = TRUE ) + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "Count", color = "p.adjust", font.size = 12, title = "") ) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "GeneRatio", color = "p.adjust", font.size = 12, title = "") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::Map(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::Map(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::Map(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + + +############# +MyEnrichDO_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { ## DO Enrichment Analysis + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- enrichDO(gene=as.data.frame(peakAnno2)$geneId, ont = "DO", pvalueCutoff = 0.05, pAdjustMethod = "BH", + qvalueCutoff = 0.2, readable = TRUE ) + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "Count", color = "p.adjust", font.size = 12, title = "") ) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "GeneRatio", color = "p.adjust", font.size = 12, title = "") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::emapplot(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + + + +############# +MyEnrichNCG_OneGroup_1_g <- function(MyMaxDis=900000000, MyPeaksAnno, MyFolder, MyFileName ) { ## Enrichment analysis based on the Network of Cancer Genes database (http://ncg.kcl.ac.uk/) + myTempFunction3 <- function() { + peakAnno2 <- as.data.frame(MyPeaksAnno) + peakAnno2 <- peakAnno2[ abs(peakAnno2$distanceToTSS) < MyMaxDis, ] + write.table(x=as.data.frame(peakAnno2), file = paste(MyFolder, "Genomic_regions_annotation_within-XXXkb.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + ## 只有这里需要改动 + enrich_raw <- enrichNCG(gene=as.data.frame(peakAnno2)$geneId, pvalueCutoff = 0.05, pAdjustMethod = "BH", + qvalueCutoff = 0.2, readable = TRUE ) + + # + write.table(x= enrich_raw , file = paste(MyFolder, "/", MyFileName, "1_enrich_raw.txt", sep=""), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + myTempFunction1234gty <- function() { + pdf(file=paste(MyFolder, "/", MyFileName, "_2_enrich_raw.pdf", sep=""), width=8, height=5) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "Count", color = "p.adjust", font.size = 12, title = "") ) + print( barplot(enrich_raw, drop=TRUE, showCategory=20, x = "GeneRatio", color = "p.adjust", font.size = 12, title = "") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 20 ) ) + print( enrichplot::Map(enrich_raw, showCategory = 20) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 20, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 10) ) + print( enrichplot::Map(enrich_raw, showCategory = 10) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 10, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichplot::dotplot(enrich_raw, showCategory = 5) ) + print( enrichplot::Map(enrich_raw, showCategory = 5) ) + print( enrichplot::cnetplot(enrich_raw, showCategory = 5, categorySize="qvalue", foldChange="GeneRatio") ) + print( enrichMap(enrich_raw, n = 10, fixed = TRUE, vertex.label.font = 1) ) + print( enrichMap(enrich_raw, n = 20, fixed = TRUE, vertex.label.font = 1) ) + dev.off() + } + tryCatch( + myTempFunction1234gty(), + error = function(err){"myTempFunction1234gty_444888"} + ) + svg(filename=paste(MyFolder, "/", MyFileName, "_3_enrich_raw.svg", sep=""), width=8, height=5) + print( enrichplot::dotplot(enrich_raw, showCategory = 30) ) + dev.off() + enrich_raw_dataFrame <- as.data.frame( enrich_raw ) + if(nrow(enrich_raw_dataFrame) > 1) { + enrich_raw_dataFrame$GeneRatio <- sapply(enrich_raw$GeneRatio, function(x) eval(parse(text=x))) * 100 + enrich_raw_dataFrame$pvalue <- -log10(enrich_raw_dataFrame$pvalue) + enrich_raw_dataFrame$p.adjust <- -log10(enrich_raw_dataFrame$p.adjust) + enrich_raw_dataFrame$qvalue <- -log10(enrich_raw_dataFrame$qvalue) + enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[1:20, -8] + if( nrow(enrich_raw_dataFrame) <20 ) {enrich_raw_dataFrame_selected <- enrich_raw_dataFrame[, -8]} + enrich_raw_dataFrame_selected <- transform(enrich_raw_dataFrame_selected, ID= rev( factor(ID, levels=unique(ID))) ) + myMidValue <- ( max(enrich_raw_dataFrame_selected$GeneRatio) + min(enrich_raw_dataFrame_selected$GeneRatio) )/2 + FigureTemp1 <- ggplot( data = enrich_raw_dataFrame_selected, aes(x = p.adjust, y = ID, size=Count , colour=GeneRatio ) ) + + geom_point( ) + scale_colour_gradient2( low = "blue", mid = "yellow2", high = "red" , midpoint = myMidValue ) + + scale_y_discrete( breaks=enrich_raw_dataFrame_selected$ID, labels=rev(enrich_raw_dataFrame_selected$ID) ) + + xlab("-log10(ajusted p-value)") + ylab("Ontology terms") + ggtitle("Ontology terms") + theme_bw() + MySaveGgplot2_1( path1=MyFolder, fileName1=paste(MyFileName, "_3_enrich_raw", sep=""), height1=6, width1=5) + } + + myTempFunction98769 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFunction98769(), + error = function(err){"myTempFunction98769_abhfd"} + ) + + + } + tryCatch( + myTempFunction3(), + error = function(err){"MyPeaksDistribution_OneGroup_1_g:error_3."} + ) +} + + +###################################################################################################################################################### + + + + + +###################################################################################################################################################### +MyAnnotation_OneGroup_1_g <- function(myFile_t1, myPath_t1) { + if( ! file.exists( myPath_t1 ) ) { dir.create(myPath_t1, recursive = TRUE) } + sink(file = paste(myPath_t1, "1_fileName.txt", sep="/") ) + print( myFile_t1 ) + sink() + myPeak_t1 = readPeakFile( myFile_t1 ) + + myPeak_t1_data.frame = as.data.frame(myPeak_t1) + sink(file = paste(myPath_t1, "2_dimension-of-inputFile.txt", sep="/") ) + print( dim(myPeak_t1_data.frame) ) + sink() + + write.table(x=myPeak_t1_data.frame, file = paste(myPath_t1, "3_Genomic-regions.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + #MyPeaksAnno_OneGroup_1_g(myPeak_t1 = myPeak_t1, myPath_t1 = paste(myPath_t1, "1_GenomicRegions-on-chromosomes", sep="/") ) + #MyPeaksSignal_OneGroup_1_g(myPeak_t1 = myPeak_t1, myPath_t1 = paste(myPath_t1, "2_GenomicRegions-on-TSSs-10kb", sep="/"), up1=upstream_g, down1=downstream_g) + MyPeaksDistribution_OneGroup_1_g(myPeak_t1 = myPeak_t1, myPath_t1 = paste(myPath_t1, "3_GenomicRegions-distribution", sep="/"), up1=upstream_g, down1=downstream_g) + #MyPeaksDistribution_OneGroup_1_g(myPeak_t1 = myPeak_t1, myPath_t1 = paste(myPath_t1, "4_GenomicRegions-distribution_x2", sep="/"), up1=upstream_g+1000, down1=downstream_g+1000) + + peakAnno <- annotatePeak( myPeak_t1, tssRegion=c(-upstream_g, downstream_g), TxDb=my_txdb_g, annoDb=my_orgdb_g, sameStrand = stranded_g , overlap="all" ) + write.table(x=as.data.frame(peakAnno), file = paste(myPath_t1, "4_Genomic_regions_annotation.txt", sep="/"), + append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE ) + + + if(fast_g == FALSE) { + my_folder_5A <- paste(myPath_t1, "5A_GeneOntology-BiologicalProcess-All", sep="/") + if( ! file.exists(my_folder_5A) ) { dir.create(my_folder_5A, recursive = TRUE) } + MyGeneOntology_BP_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_5A, MyFileName="5A_all-NearestGenes_") + + my_folder_5B <- paste(myPath_t1, "5B_GeneOntology-BiologicalProcess-within100kb", sep="/") + if( ! file.exists(my_folder_5B) ) { dir.create(my_folder_5B, recursive = TRUE) } + MyGeneOntology_BP_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_5B, MyFileName="5B_100kb-NearestGenes_") + + my_folder_5C <- paste(myPath_t1, "5C_GeneOntology-BiologicalProcess-within50kb", sep="/") + if( ! file.exists(my_folder_5C) ) { dir.create(my_folder_5C, recursive = TRUE) } + MyGeneOntology_BP_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_5C, MyFileName="5C_50kb-NearestGenes_") + + my_folder_5D <- paste(myPath_t1, "5D_GeneOntology-BiologicalProcess-within20kb", sep="/") + if( ! file.exists(my_folder_5D) ) { dir.create(my_folder_5D, recursive = TRUE) } + MyGeneOntology_BP_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_5D, MyFileName="5D_20kb-NearestGenes_") + + my_folder_5E <- paste(myPath_t1, "5E_GeneOntology-BiologicalProcess-within5kb", sep="/") + if( ! file.exists(my_folder_5E) ) { dir.create(my_folder_5E, recursive = TRUE) } + MyGeneOntology_BP_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_5E, MyFileName="5E_5kb-NearestGenes_") + + + + my_folder_6A <- paste(myPath_t1, "6A_GeneOntology-MolecularFunction-All", sep="/") + if( ! file.exists(my_folder_6A) ) { dir.create(my_folder_6A, recursive = TRUE) } + MyGeneOntology_MF_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_6A, MyFileName="6A_all-NearestGenes_") + + my_folder_6B <- paste(myPath_t1, "6B_GeneOntology-MolecularFunction-within100kb", sep="/") + if( ! file.exists(my_folder_6B) ) { dir.create(my_folder_6B, recursive = TRUE) } + MyGeneOntology_MF_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_6B, MyFileName="6B_100kb-NearestGenes_") + + my_folder_6C <- paste(myPath_t1, "6C_GeneOntology-MolecularFunction-within50kb", sep="/") + if( ! file.exists(my_folder_6C) ) { dir.create(my_folder_6C, recursive = TRUE) } + MyGeneOntology_MF_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_6C, MyFileName="6C_50kb-NearestGenes_") + + my_folder_6D <- paste(myPath_t1, "6D_GeneOntology-MolecularFunction-within20kb", sep="/") + if( ! file.exists(my_folder_6D) ) { dir.create(my_folder_6D, recursive = TRUE) } + MyGeneOntology_MF_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_6D, MyFileName="6D_20kb-NearestGenes_") + + my_folder_6E <- paste(myPath_t1, "6E_GeneOntology-MolecularFunction-within5kb", sep="/") + if( ! file.exists(my_folder_6E) ) { dir.create(my_folder_6E, recursive = TRUE) } + MyGeneOntology_MF_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_6E, MyFileName="6E_5kb-NearestGenes_") + + + + + my_folder_7A <- paste(myPath_t1, "7A_GeneOntology-CellularComponent-All", sep="/") + if( ! file.exists(my_folder_7A) ) { dir.create(my_folder_7A, recursive = TRUE) } + MyGeneOntology_CC_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_7A, MyFileName="7A_all-NearestGenes_") + + my_folder_7B <- paste(myPath_t1, "7B_GeneOntology-CellularComponent-within100kb", sep="/") + if( ! file.exists(my_folder_7B) ) { dir.create(my_folder_7B, recursive = TRUE) } + MyGeneOntology_CC_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_7B, MyFileName="7B_100kb-NearestGenes_") + + my_folder_7C <- paste(myPath_t1, "7C_GeneOntology-CellularComponent-within50kb", sep="/") + if( ! file.exists(my_folder_7C) ) { dir.create(my_folder_7C, recursive = TRUE) } + MyGeneOntology_CC_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_7C, MyFileName="7C_50kb-NearestGenes_") + + my_folder_7D <- paste(myPath_t1, "7D_GeneOntology-CellularComponent-within20kb", sep="/") + if( ! file.exists(my_folder_7D) ) { dir.create(my_folder_7D, recursive = TRUE) } + MyGeneOntology_CC_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_7D, MyFileName="7D_20kb-NearestGenes_") + + my_folder_7E <- paste(myPath_t1, "7E_GeneOntology-CellularComponent-within5kb", sep="/") + if( ! file.exists(my_folder_7E) ) { dir.create(my_folder_7E, recursive = TRUE) } + MyGeneOntology_CC_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_7E, MyFileName="7E_5kb-NearestGenes_") + + + + + my_folder_8A <- paste(myPath_t1, "8A_Reactome-All", sep="/") + if( ! file.exists(my_folder_8A) ) { dir.create(my_folder_8A, recursive = TRUE) } + MyReactome_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_8A, MyFileName="8A_all-NearestGenes_") + + my_folder_8B <- paste(myPath_t1, "8B_Reactome-within100kb", sep="/") + if( ! file.exists(my_folder_8B) ) { dir.create(my_folder_8B, recursive = TRUE) } + MyReactome_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_8B, MyFileName="8B_100kb-NearestGenes_") + + my_folder_8C <- paste(myPath_t1, "8C_Reactome-within50kb", sep="/") + if( ! file.exists(my_folder_8C) ) { dir.create(my_folder_8C, recursive = TRUE) } + MyReactome_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_8C, MyFileName="8C_50kb-NearestGenes_") + + my_folder_8D <- paste(myPath_t1, "8D_Reactome-within20kb", sep="/") + if( ! file.exists(my_folder_8D) ) { dir.create(my_folder_8D, recursive = TRUE) } + MyReactome_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_8D, MyFileName="8D_20kb-NearestGenes_") + + my_folder_8E <- paste(myPath_t1, "8E_Reactome-within5kb", sep="/") + if( ! file.exists(my_folder_8E) ) { dir.create(my_folder_8E, recursive = TRUE) } + MyReactome_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_8E, MyFileName="8E_5kb-NearestGenes_") + + + + + my_folder_9A <- paste(myPath_t1, "9A_KEGG-All", sep="/") + if( ! file.exists(my_folder_9A) ) { dir.create(my_folder_9A, recursive = TRUE) } + MyKEGG_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_9A, MyFileName="9A_all-NearestGenes_") + + my_folder_9B <- paste(myPath_t1, "9B_KEGG-within100kb", sep="/") + if( ! file.exists(my_folder_9B) ) { dir.create(my_folder_9B, recursive = TRUE) } + MyKEGG_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_9B, MyFileName="9B_100kb-NearestGenes_") + + my_folder_9C <- paste(myPath_t1, "9C_KEGG-within50kb", sep="/") + if( ! file.exists(my_folder_9C) ) { dir.create(my_folder_9C, recursive = TRUE) } + MyKEGG_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_9C, MyFileName="9C_50kb-NearestGenes_") + + my_folder_9D <- paste(myPath_t1, "9D_KEGG-within20kb", sep="/") + if( ! file.exists(my_folder_9D) ) { dir.create(my_folder_9D, recursive = TRUE) } + MyKEGG_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_9D, MyFileName="9D_20kb-NearestGenes_") + + my_folder_9E <- paste(myPath_t1, "9E_KEGG-within5kb", sep="/") + if( ! file.exists(my_folder_9E) ) { dir.create(my_folder_9E, recursive = TRUE) } + MyKEGG_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_9E, MyFileName="9E_5kb-NearestGenes_") + + + + my_folder_10A <- paste(myPath_t1, "10A_ModuleKEGG-All", sep="/") + if( ! file.exists(my_folder_10A) ) { dir.create(my_folder_10A, recursive = TRUE) } + MyModuleKEGG_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_10A, MyFileName="10A_all-NearestGenes_") + + my_folder_10B <- paste(myPath_t1, "10B_ModuleKEGG-within100kb", sep="/") + if( ! file.exists(my_folder_10B) ) { dir.create(my_folder_10B, recursive = TRUE) } + MyModuleKEGG_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_10B, MyFileName="10B_100kb-NearestGenes_") + + my_folder_10C <- paste(myPath_t1, "10C_ModuleKEGG-within50kb", sep="/") + if( ! file.exists(my_folder_10C) ) { dir.create(my_folder_10C, recursive = TRUE) } + MyModuleKEGG_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_10C, MyFileName="10C_50kb-NearestGenes_") + + my_folder_10D <- paste(myPath_t1, "10D_ModuleKEGG-within20kb", sep="/") + if( ! file.exists(my_folder_10D) ) { dir.create(my_folder_10D, recursive = TRUE) } + MyModuleKEGG_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_10D, MyFileName="10D_20kb-NearestGenes_") + + my_folder_10E <- paste(myPath_t1, "10E_ModuleKEGG-within5kb", sep="/") + if( ! file.exists(my_folder_10E) ) { dir.create(my_folder_10E, recursive = TRUE) } + MyModuleKEGG_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_10E, MyFileName="10E_5kb-NearestGenes_") + + + + my_folder_11A <- paste(myPath_t1, "11A_DiseaseEnrichmentOnDisGeNET-All", sep="/") + if( ! file.exists(my_folder_11A) ) { dir.create(my_folder_11A, recursive = TRUE) } + MyEnrichDGN_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_11A, MyFileName="11A_all-NearestGenes_") + + my_folder_11B <- paste(myPath_t1, "11B_DiseaseEnrichmentOnDisGeNET-within100kb", sep="/") + if( ! file.exists(my_folder_11B) ) { dir.create(my_folder_11B, recursive = TRUE) } + MyEnrichDGN_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_11B, MyFileName="11B_100kb-NearestGenes_") + + my_folder_11C <- paste(myPath_t1, "11C_DiseaseEnrichmentOnDisGeNET-within50kb", sep="/") + if( ! file.exists(my_folder_11C) ) { dir.create(my_folder_11C, recursive = TRUE) } + MyEnrichDGN_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_11C, MyFileName="11C_50kb-NearestGenes_") + + my_folder_11D <- paste(myPath_t1, "11D_DiseaseEnrichmentOnDisGeNET-within20kb", sep="/") + if( ! file.exists(my_folder_11D) ) { dir.create(my_folder_11D, recursive = TRUE) } + MyEnrichDGN_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_11D, MyFileName="11D_20kb-NearestGenes_") + + my_folder_11E <- paste(myPath_t1, "11E_DiseaseEnrichmentOnDisGeNET-within5kb", sep="/") + if( ! file.exists(my_folder_11E) ) { dir.create(my_folder_11E, recursive = TRUE) } + MyEnrichDGN_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_11E, MyFileName="11E_5kb-NearestGenes_") + + + + + my_folder_12A <- paste(myPath_t1, "12A_DiseaseEnrichmentOnDO-All", sep="/") + if( ! file.exists(my_folder_12A) ) { dir.create(my_folder_12A, recursive = TRUE) } + MyEnrichDO_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_12A, MyFileName="12A_all-NearestGenes_") + + my_folder_12B <- paste(myPath_t1, "12B_DiseaseEnrichmentOnDO-within100kb", sep="/") + if( ! file.exists(my_folder_12B) ) { dir.create(my_folder_12B, recursive = TRUE) } + MyEnrichDO_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_12B, MyFileName="12B_100kb-NearestGenes_") + + my_folder_12C <- paste(myPath_t1, "12C_DiseaseEnrichmentOnDO-within50kb", sep="/") + if( ! file.exists(my_folder_12C) ) { dir.create(my_folder_12C, recursive = TRUE) } + MyEnrichDO_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_12C, MyFileName="12C_50kb-NearestGenes_") + + my_folder_12D <- paste(myPath_t1, "12D_DiseaseEnrichmentOnDO-within20kb", sep="/") + if( ! file.exists(my_folder_12D) ) { dir.create(my_folder_12D, recursive = TRUE) } + MyEnrichDO_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_12D, MyFileName="12D_20kb-NearestGenes_") + + my_folder_12E <- paste(myPath_t1, "12E_DiseaseEnrichmentOnDO-within5kb", sep="/") + if( ! file.exists(my_folder_12E) ) { dir.create(my_folder_12E, recursive = TRUE) } + MyEnrichDO_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_12E, MyFileName="12E_5kb-NearestGenes_") + + + + my_folder_13A <- paste(myPath_t1, "13A_DiseaseEnrichmentOnTheNetworkOfCancerGenesDatabase-All", sep="/") + if( ! file.exists(my_folder_13A) ) { dir.create(my_folder_13A, recursive = TRUE) } + MyEnrichNCG_OneGroup_1_g(MyMaxDis=900000000, MyPeaksAnno=peakAnno, MyFolder=my_folder_13A, MyFileName="13A_all-NearestGenes_") + + my_folder_13B <- paste(myPath_t1, "13B_DiseaseEnrichmentOnTheNetworkOfCancerGenesDatabase-within100kb", sep="/") + if( ! file.exists(my_folder_13B) ) { dir.create(my_folder_13B, recursive = TRUE) } + MyEnrichNCG_OneGroup_1_g(MyMaxDis=100000, MyPeaksAnno=peakAnno, MyFolder=my_folder_13B, MyFileName="13B_100kb-NearestGenes_") + + my_folder_13C <- paste(myPath_t1, "13C_DiseaseEnrichmentOnTheNetworkOfCancerGenesDatabase-within50kb", sep="/") + if( ! file.exists(my_folder_13C) ) { dir.create(my_folder_13C, recursive = TRUE) } + MyEnrichNCG_OneGroup_1_g(MyMaxDis=50000, MyPeaksAnno=peakAnno, MyFolder=my_folder_13C, MyFileName="13C_50kb-NearestGenes_") + + my_folder_13D <- paste(myPath_t1, "13D_DiseaseEnrichmentOnTheNetworkOfCancerGenesDatabase-within20kb", sep="/") + if( ! file.exists(my_folder_13D) ) { dir.create(my_folder_13D, recursive = TRUE) } + MyEnrichNCG_OneGroup_1_g(MyMaxDis=20000, MyPeaksAnno=peakAnno, MyFolder=my_folder_13D, MyFileName="13D_20kb-NearestGenes_") + + my_folder_13E <- paste(myPath_t1, "13E_DiseaseEnrichmentOnTheNetworkOfCancerGenesDatabase-within5kb", sep="/") + if( ! file.exists(my_folder_13E) ) { dir.create(my_folder_13E, recursive = TRUE) } + MyEnrichNCG_OneGroup_1_g(MyMaxDis=5000, MyPeaksAnno=peakAnno, MyFolder=my_folder_13E, MyFileName="13E_5kb-NearestGenes_") + } + + myTempFun_1 <- function() { + dev.off() + dev.off() + dev.off() + dev.off() + dev.off() + } + tryCatch( + myTempFun_1(), + error = function(err){"myTempFun_1_YongPeng."} + ) + +} +###################################################################################################################################################### + + + + + +###################################################################################################################################################### +for(i in c(1:length(peakFiles)) ) { + myTempFunction3 <- function() { + ##MyAnnotation_OneGroup_1_g( myFile_t1 = peakFiles[i], myPath_t1 = paste(outDir_g, "/File.", i, "._", peakFiles_onlyName[i], sep="") ) + } + tryCatch( + myTempFunction3(), + error = function(err){"My_OneGroup_1_g:error_3."} + ) +} +###################################################################################################################################################### + + + + + + + + + + +########################################### +myPath_multi_1 = paste(outDir_g, "Multiple_Files",sep="/") +if( ! file.exists( myPath_multi_1 ) ) { dir.create(myPath_multi_1, recursive = TRUE) } + + +peakAnnoList_1 <- lapply(peakFiles, annotatePeak, TxDb=my_txdb_g, annoDb=my_orgdb_g, tssRegion=c(-upstream_g, downstream_g), verbose=FALSE, sameStrand = stranded_g , overlap="all") +names( peakAnnoList_1 ) <- peakFiles_onlyName + +pdf( file=paste(myPath_multi_1, "Genomic_Locations.pdf", sep="/"), width=10, height=length(peakFiles_onlyName)/4 + 1 ) + plotAnnoBar(peakAnnoList_1) + plotDistToTSS(peakAnnoList_1) +dev.off() + +sink( file = paste(myPath_multi_1, "Genomic_Locations.txt", sep="/") ) +print(peakAnnoList_1) +sink() + +########################################### + + + + + + + + + diff --git a/Call_Mutations/figures/Anno_v20210106.sh b/Call_Mutations/figures/Anno_v20210106.sh new file mode 100644 index 0000000..c853f8d --- /dev/null +++ b/Call_Mutations/figures/Anno_v20210106.sh @@ -0,0 +1,14 @@ + + + Rscript Anno_v20210106.R \ + --inputDir=2-getFraction \ + --outDir=4-Annotation \ + --refGenome=hg38 \ + --stranded=TRUE \ + --fast=TRUE \ + --upstream=1000 \ + --downstream=0 + + + + diff --git a/Call_Mutations/figures/figures.pl b/Call_Mutations/figures/figures.pl new file mode 100644 index 0000000..41730c8 --- /dev/null +++ b/Call_Mutations/figures/figures.pl @@ -0,0 +1,118 @@ +#!/usr/bin/env perl +use strict; +use warnings; +use v5.12; +## Perl5 version >= 5.12 +## Suffixes of all self-defined global variables must be "_g". +################################################################################################################################################################################################### + + + + + +################################################################################################################################################################################################### +my $input_g = ''; ## such as "1-finalSites", global variable. +my $output_g = ''; ## such as "2-getFraction", global variable. +{ +## Help Infromation +my $HELP = ' + ------------------------------------------------------------------------------------------------------------------------------------------------------ + ------------------------------------------------------------------------------------------------------------------------------------------------------ + + Usage: + perl figures.pl [-version] [-help] [-in inputFile] [-out outFile] + For instance: + perl figures.pl -in 1-finalSites -out 2-getFraction > figures.runLog.txt + + Optional arguments: + -version Show version number of this program and exit. + -help Show this help message and exit. + + Required arguments: + -in inputFile "inputFile" is the name of your input file path. (no default) + -out outFile "outFile" is the name of output path. (no default) + ------------------------------------------------------------------------------------------------------------------------------------------------------ + ------------------------------------------------------------------------------------------------------------------------------------------------------ +'; + +## Version Infromation +my $version = " Version 0.2, 2020-08-01."; + +## Keys and Values +if ($#ARGV == -1) { say "\n$HELP\n"; exit 0; } ## when there are no any command argumants. +if ($#ARGV%2 == 0) { @ARGV = (@ARGV, "-help") ; } ## when the number of command argumants is odd. +my %args = @ARGV; + +## Initialize Variables +$input_g = '1-finalSites'; ## This is only an initialization value or suggesting value, not default value. +$output_g = '2-getFraction'; ## This is only an initialization value or suggesting value, not default value. + +## Available Arguments +my $available = " -version -help -in -out "; +my $boole = 0; +while( my ($key, $value) = each %args ) { + if ( ($key =~ m/^\-/) and ($available !~ m/\s$key\s/) ) {say "\n\tCann't recognize $key"; $boole = 1; } +} +if($boole == 1) { + say "\tThe Command Line Arguments are wrong!"; + say "\tPlease see help message by using 'perl figures.pl -help' \n"; + exit 0; +} + +## Get Arguments +if ( exists $args{'-version' } ) { say "\n$version\n"; exit 0; } +if ( exists $args{'-help' } ) { say "\n$HELP\n"; exit 0; } +if ( exists $args{'-in' } ) { $input_g = $args{'-in' }; }else{say "\n -in is required.\n"; say "\n$HELP\n"; exit 0; } +if ( exists $args{'-out' } ) { $output_g = $args{'-out' }; }else{say "\n -out is required.\n"; say "\n$HELP\n"; exit 0; } + +## Conditions +$input_g =~ m/^\S+$/ || die "\n\n$HELP\n\n"; +$output_g =~ m/^\S+$/ || die "\n\n$HELP\n\n"; + +## Print Command Arguments to Standard Output +say "\n + ################ Arguments ############################### + Input Path: $input_g + Output Path: $output_g + ########################################################## +\n"; +} +################################################################################################################################################################################################### + + + + +opendir(my $FH_input1_g, $input_g) || die; +my @files = readdir($FH_input1_g); + + + + + + +################################################################################################################################################################################################### + +for(my $i=0; $i<=$#files; $i=$i+1) { + my $temp = $files[$i]; + next unless $temp =~ m/\.bed$/; + print "$temp\n"; + system("perl 2-getFraction.pl -f spikeins.fit.txt -in $input_g/$temp -out $output_g/$temp "); +} + +################################################################################################################################################################################################### + + + + + + + + + + + + + + + + diff --git a/Call_Mutations/figures/heatmap.R b/Call_Mutations/figures/heatmap.R new file mode 100644 index 0000000..99f9ed2 --- /dev/null +++ b/Call_Mutations/figures/heatmap.R @@ -0,0 +1,157 @@ +suppressPackageStartupMessages( library(ComplexHeatmap) ) +suppressPackageStartupMessages( library(corrplot) ) +library(gplots) + +reset_outliers2 <- function(x, na.rm = TRUE ) { + qnt <- quantile(x, probs=c(0.1, 0.9) , type=1, na.rm = na.rm ) + y <- x + y[x < qnt[1] ] <- qnt[1] + y[x > qnt[2] ] <- qnt[2] + y +} + + + +myScaleMatrix2 <- function( matrix_temp8, upper_temp8 = 1, lower_temp8 = -1 ) { + rawMatrix_2 = reset_outliers2(matrix_temp8) + rawMatrix_2 = lower_temp8 + (upper_temp8 - lower_temp8) * ( rawMatrix_2 - min(rawMatrix_2) )/( max(rawMatrix_2)- min(rawMatrix_2) ) + return(rawMatrix_2) +} + + + +MyHeatmaps_1_f <- function(matrix2, path2, fileName2, height2=30, width2=30, is.corr2=TRUE, my_col2 ) { + matrix2 = myScaleMatrix2( matrix2 ) + pdf( file = paste(path2, fileName2, sep="/"), width=width2, height=height2 ) + print( corrplot(matrix2, method = "circle", type = "full", title = "", is.corr = is.corr2, order = "hclust", hclust.method = "ward.D2", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "circle", type = "upper", title = "", is.corr = is.corr2, order = "hclust", hclust.method = "ward.D2", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "number", type = "upper", title = "", is.corr = is.corr2, order = "hclust", hclust.method = "ward.D2", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "pie", type = "upper", title = "", is.corr = is.corr2, order = "hclust", hclust.method = "ward.D2", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "color", type = "upper", title = "", is.corr = is.corr2, order = "hclust", hclust.method = "ward.D2", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "color", type = "full", title = "", is.corr = is.corr2, order = "hclust", hclust.method = "ward.D2", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + dev.off() + +} + + + +MyHeatmaps_1A_f <- function(matrix2, path2, fileName2, height2=30, width2=30, is.corr2=TRUE, my_col2 ) { + ## matrix2 = myScaleMatrix2( matrix2 ) + pdf( file = paste(path2, fileName2, sep="/"), width=width2, height=height2 ) + print( corrplot(matrix2, method = "circle", type = "full", title = "", is.corr = is.corr2, order = "original", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "circle", type = "upper", title = "", is.corr = is.corr2, order = "original", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "number", type = "full", title = "", is.corr = is.corr2, order = "original", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "pie", type = "full", title = "", is.corr = is.corr2, order = "original", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "color", type = "full", title = "", is.corr = is.corr2, order = "original", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + print( corrplot(matrix2, method = "color", type = "full", title = "", is.corr = is.corr2, order = "original", tl.col = "black", tl.srt = 45, col = my_col2 ) ) + dev.off() + +} + + + +MyHeatmaps_2_f <- function(matrix2, path2, fileName2, height2=30, width2=30 ) { + matrix2 = myScaleMatrix2( matrix2 ) + pdf( file = paste(path2, fileName2, sep="/"), width=width2, height=height2 ) + print( heatmap(x = matrix2, col = colorRampPalette(c("blue", "white", "red"))(20), symm = FALSE, scale = "none" ) ) + print( heatmap(x = matrix2, col = colorRampPalette(c("green4", "white", "purple"))(20), symm = FALSE, scale = "none" ) ) + print( heatmap(x = matrix2, col = colorRampPalette(c("cyan", "white", "red"))(20), symm = FALSE, scale = "none" ) ) + print( heatmap(x = matrix2, col = colorRampPalette(c("cyan", "white", "purple"))(20), symm = FALSE, scale = "none" ) ) + print( heatmap(x = matrix2, col = colorRampPalette(c("blue", "white", "purple"))(20), symm = FALSE, scale = "none" ) ) + dev.off() + +} + + +MyHeatmaps_3_f <- function(matrix2, path2, fileName2, height2=30, width2=30 ) { + matrix2 = myScaleMatrix2( matrix2 ) + pdf( file = paste(path2, fileName2, sep="/"), width=width2, height=height2 ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("blue", "white", "red"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("blue", "white", "red"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("green4", "white", "purple"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("green4", "white", "purple"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("cyan", "white", "red"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("cyan", "white", "red"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("cyan", "white", "purple"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("cyan", "white", "purple"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("blue", "white", "purple"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + print( ComplexHeatmap::Heatmap(matrix2, col = colorRampPalette(c("blue", "white", "purple"))(20) , heatmap_legend_param = list(legend_height = unit(10, "cm")) ) ) + dev.off() + +} + + + +################################################# +rawMatrix_1 <- read.table("Intervene_pairwise_frac_matrix.txt", header=TRUE, sep="\t", quote = "", comment.char = "") +rawMatrix_1 = rawMatrix_1[,-1] +rownames(rawMatrix_1) +colnames(rawMatrix_1) +dim(rawMatrix_1) + + +outPath = "figures" +if( ! file.exists(outPath) ) { dir.create(outPath, recursive = TRUE) } + + + +rawMatrix_2 <- rawMatrix_1 +dim(rawMatrix_2) + + +for(i in c(1:nrow(rawMatrix_2))) { + rawMatrix_2[i,i] = 0 +} + + + +my_col1=colorRampPalette( c( "black", "black", "black", "pink", "red", "red", "red4") ) + + +MyHeatmaps_1A_f(matrix2=as.matrix(rawMatrix_2), path2=outPath, fileName2="heatmap.pdf", + height2=9, width2=6, is.corr2=FALSE, my_col2=my_col1(10) ) + + + + + +rawMatrix_2 = rawMatrix_2*100 +rawMatrix_3 = round(rawMatrix_2, digits = 0) +rawMatrix_3 + + +my_col2=colorRampPalette( c( "grey80", "yellow3", "red" ), bias = 1.5 ) + + +pdf( file="heatmap-3.pdf", height=5, width=5 ) +heatmap.2 (x=as.matrix(rawMatrix_3), + + # dendrogram control + dendrogram = "none", + Rowv = FALSE, + Colv="Rowv" , + symm = FALSE, + + # data scaling + scale = "none" , + na.rm=TRUE, + + # colors + col=my_col2(40), + trace = "none", + + # cell labeling + cellnote= as.matrix(rawMatrix_3) , + notecex=1.0, + notecol="white", + na.color=par("bg") + +) +dev.off() + + + + + + + diff --git a/Call_Mutations/figures/histogram.R b/Call_Mutations/figures/histogram.R new file mode 100644 index 0000000..5e1550b --- /dev/null +++ b/Call_Mutations/figures/histogram.R @@ -0,0 +1,39 @@ + + +################################################# + +outDir_g = "histogram" +if( ! file.exists(outDir_g) ) { dir.create(outDir_g, recursive = TRUE) } + +options(digits=10) +continue_on_error <- function() { + print( "NOTE: THERE WAS AN ERROR HERE. We are continuing because we have set 'options(error=continue_on_error())' " ) +} +options( error=continue_on_error ) ## This option is very important. +############################################# + + + + + + +############################################# +matrix1 = read.table(file= "number.txt", header = F, sep = character(0), quote = "\"'", dec = ".") +dim(matrix1) + + +ratio1 = as.numeric( matrix1[,1] ) +col_names1 = matrix1[,2] + + +pdf(paste(outDir_g, "number_of_peaks.pdf", sep="/"), height=5, width=10) + +xx1 <- barplot( ratio1, xaxt = 'n', xlab = 'samples', ylab = "number", col="red", width = 1 ) +text(x = xx1, y = ratio1, label = ratio1, pos = 3, cex = 0.9, col = "black") +axis(1, at=xx1, labels=col_names1, tick=FALSE, las=2, line=0, cex.axis=0.6, cex=5 ) + +dev.off() + + + + \ No newline at end of file