售后

脚本

perl peak_overlap_venn.pl \
  -peakfile peakfile.txt \
  -outdir . \
  -mode peak

具体脚本

#!/usr/bin/perl
use Getopt::Long;
use File::Basename qw(basename dirname);

my $ver    = "1.0";
my $Writer = "chenzixi\@novogene.com";
my $Data   = "2019/12/25";
#######################################################################################
# ------------------------------------------------------------------
# GetOptions
# ------------------------------------------------------------------
my ($peakfile,$outdir,$mode);
GetOptions(
                        "h|?" =>\&help,
                        "peakfile:s"=>\$peakfile,
                        "outdir:s"=>\$outdir,
						"mode:s"=>\$mode,
                        ) || &help;
&help unless ($peakfile && $outdir );
sub help
{
        print <<"    Usage End.";
        Description:
        Writer  : $Writer
        Data    : $Data
        Version : $ver
        function: ......
        Usage:
                -peakfile    peakfile, eg /TJPROJ1/RNA/WORK/Pipline/aftersale/peak_overlap_venn/testdata/peakfile.txt   Required

                -outdir      output result dir                                                                          Required

				-mode		 peak or dml[peak]																		Optional

                -h           Help                                                                                       document

#################
for DML mode,use both start/end as the position(start=end)
#################

    Usage End.
        exit;
}
# ------------------------------------------------------------------
# GetOptions End
# ------------------------------------------------------------------
my @pallete=qw/#E43B2D #377EB7 #4DAF4A #984EA3 #F27E33 #FEF840 #A65629 #F180BF #999999/;
$mode||="peak";
my $min=1;

if (($mode ne "peak")&&($mode ne "dml")){
	die "only peak/dml mode suppported, please check!\n";
}


my %samplelist;
my @samples;
my $one_dml_file;
open (IN,"< $peakfile") || die $!;
while(<IN>){
    chomp;
    my @line=split /\t/;
    $samplelist{$line[0]}=$line[1];
	$one_dml_file=$line[1];
    push @samples,$line[0];
}
close IN;

### Test start=end if using DML mode
if ($mode eq "dml"){
	open (IN,"$one_dml_file");
	my $line1=<IN>;
	close IN;
	chomp $line1;
	my @line=split(/\t/,$line1);
	die "If using dml mode, the start and end position should be the same. Please check!\n" if($line[1] != $line[2]);

}


my $remove=scalar @pallete - scalar @samples;#print $remove;
for(my $i=0;$i<$remove;$i++){
    pop @pallete;
}
my $colors="\"".join("\",\"",@pallete)."\"";

chdir($outdir);
my $sh="$outdir\/Peak_Overlap.r";
open (SH,"> $sh") || die $!;
print SH <<"  END.";
library(ChIPpeakAnno)

  END.

my $i=0;
foreach my $sample (@samples){
    $i++;
    print SH <<"  END.";
###  read file for $sample
data$i<-read.delim("$samplelist{$sample}",header=F,sep = "\\t")
$sample<-unique(GRanges(seqnames=data$i\[,1],
               ranges=IRanges(start=data$i\[,2],
                              end=data$i\[,3],
                              names=paste(data$i\[,1],data$i\[,2],data$i\[,3],sep=":")),
               strand=rep("+",length(data$i\[,1]))))
  END.
}

my $allsample=join(",",@samples);
print SH <<"  END.";
### overlap peaks
ol <- findOverlapsOfPeaks($allsample,minoverlap=$min,
                          connectedPeaks="merge")
### make venn
pdf("$outdir\/Peak_Overlap.pdf")
makeVennDiagram(ol, col="#00000000",euler.d=TRUE,scaled=TRUE,
                fill=c($colors),
                margin=c(0.05,0.05,0.05,0.05))
dev.off()
system("convert -density 1200 -resize 600 $outdir\/Peak_Overlap.pdf $outdir\/Peak_Overlap.png")
### output overlap peaklist
for (i in names(ol\$peaklist)){
  change<-gsub("///","_",i)
  filename<-paste(change,"all.id.xls",sep=".")
  temp<-as.data.frame(ol\$peaklist[[i]])
  for (j in c(1:length(temp\$peakNames))){
    temp\$samplepeaks[j]<-paste(as.vector(temp\$peakNames[[j]]),collapse = ";")
  }
  temp <- subset(temp, select = -peakNames)
  write.table(temp,filename,sep = "\t",quote=F,row.names = F,col.names = F)
  rm(temp)
}
  END.
close SH;
#`/PUBLIC/software/RNA/R/R-3.5.1/bin/Rscript $sh`;
`source /TJPROJ1/RNA/WORK/software/ChIPpeakAnno_dev/source_before_use;/TJPROJ1/RNA/WORK/software/ChIPpeakAnno_dev/R-4.0.2/bin/Rscript $sh`;
#`echo source /TJPROJ1/RNA/WORK/software/ChIPpeakAnno_dev/source_before_use > $outdir\/Peak_Overlap.sh`;
#`echo /TJPROJ1/RNA/WORK/software/ChIPpeakAnno_dev/R-4.0.2/bin/Rscript $sh >> $outdir\/Peak_Overlap.sh`;
#`sh $outdir\/Peak_Overlap.sh`;
`ls $outdir\/*.xls|xargs sed -i '1ichr\tstart\tend\tlength\tstrand\tsamplepeak'`;