diff [APliBio]Nebula tools suite/Nebula/CreateControlSubSet/createControlSubSet.pl @ 0:2ec3ba0e9e70 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/[APliBio]Nebula tools suite/Nebula/CreateControlSubSet/createControlSubSet.pl	Thu Oct 25 08:18:25 2012 -0400
@@ -0,0 +1,280 @@
+#:t:::::::::::::::::g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+#:t::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+#:::::::::::::z;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+#::::::::::::i@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+#::::::::::::@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@$@@@@
+#:::::::::::3@@@@@@@@@@@@@@@@@@@@@@@@@B@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+#::::::::::3@@@@@@@@@@@@@@@@@@@@@BEEESSE5EEEEBBM@@@@@@@@@@@@@@@@@@@@@@@@@@
+#::::::::::3@@@@@@@@@@@@@@@@@@@@BEEEEEE35EE55E2355E5SBMB@@@@@@@@@@@@@@@@@$
+#::::::::::@@@@@@@@@@@@@@@@@@@EEEE55533t3tttt::::::!!!!7755E755SBBMMM@@@MM
+#::::::::::3@@@@@@@@@@@@@@@@@@EEEE2t3ttttt:::::::::::::::::::::::!7?5225EE
+#::::::::::3@@@@@@@@@@@@@@@@@@EEEEE31t::::::::::::::::::::::::::::::::3E5@
+#::::::::::3@@@@@@@@@@@@@@@@@@EEEEEEtt:::::::::::::::::::::::::::::::::353
+#::::::::::3@@@@@@@@@@@@@@@@@@EEEEEE1ttz::::::::::::::::::::::::::::::::35
+#:::::::::::@@@@@@@@@@@@@@@@@@EEEEEEEtz1::::::::::::::::::::::::::::::::t:
+#:::::::::!3@@@@@@@@@@@@@@@@@@@EEEEEttt::::::::::::::::::::::::::::::::;zz
+#::::::::::@@@@@@@@@@@@@@@@@@@@EEEEEttt:::::z;z:::::::::::::::::::::::::13
+#::::::::::3B@@@@@@@@@@@@@@@@@@EEEEEEE3tt:czzztti;:::::::::::::::::::::::3
+#::::ttt::::3@@@@@@@@@@@@@@@@EEEEE5EE25Ezt1EEEz5Etzzz;;;;:::::::::::::::::
+#:::::::::::I9@@@@@@@@@@@@@@@@@@@@@@@@@@EEEEEE@@@@@@@@@@@@@@Ez;:::::::::::
+#:::::::::::::E@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@Ez::::::
+#::::::::::::::E@@@@@@@@@@@@@@@@@@@@@@@@@@@@@BE5EBB@@@@@@@@@@@@@@@EEE:::::
+#:::::::::::::::@@@@@@@@@@@@@@@@@@@@@@@@@@@@E1::35@@@@@@@@@@ME3MMME2::::::
+#:::::::::::::::?@@@@@@@@@@@@@@@@@@M@@@@@@@EE:::::3SB@@BBESEEt::::::::::::
+#::::::::::::::::J$@@@@@@@B@@@@@@@@@@@@@@@@EE:::::::!35E33t:::::::::::::::
+#:::::::::::::::::3@E@@@EE5EESE5EESE@@@@@@@Et::::::::::::tz:::::::::::::::
+#:::::::::::::::::J@E$@EEE5133555SE@@@@@@@@Et:::::::::::::::::::::::::::::
+#::::::::::::::::::E@E@EEEEtt3523EEE@@@@@@@E::::::::::::::::::::::::::::::
+#:t::::::::::::::::JEE3@@@EEEEEEEEEE@@@@@@@E:::::::::t;:::::::::::::::::::
+#:t:::::::::::::::::!5ES@EEEEEEEEES@@@@@@@@@E;:::;;;:3Ez::::::::::::::::::
+#:t::::::::::::::::::::JE@@EEEEEEE@@@@@@@@@@@@@@@@ME!:::;:::::::::::::::::
+#:tz::::::::::::::::::::JE@@@EEEE@@@@@@@@@@@@@@EE!:::::::t::::::::::::::::
+#:t::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@ESBE::::::::::::::::::::::::::
+#:::::::::::::::::::::::::Q@@@@@@@@@@@@@@@@EE3EE;:::::zzzz::::::::::::::::
+#:::::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@@@@@NN@@@@@@Ez:::::::::::::::
+#:zt:::::::::::::::::::::::3@@@@EE@@@@@@@@@@EEEEt::;z113E5t:::::::::::::::
+#::tt:::::::::::::::::::::::3@@@E@@@@@@@@@@@@@@@@BEt::::::::::::::::t:::::
+#:tt:t:::::::::::::::::::::::?S@@@@@@@@@@@BBEEE51!::::::::::::::zzzEt:::::
+#::::::::::::::::::::::::::::::3Q@@@@@@@BEEEEEt:::::::::::::;zz@@@EE::::::
+#::::::::::::::::::::::::::::::::75B@@@@@EEEtt;:::::::::;zz@@@@BEEEtz:::::
+#::::::::::::::::::::::::::::::::::::?9@@@@@@@@@@@E2Ezg@@@@@B@@@EEEE1t::::
+#:::::::::::::::::::::::::::::::::::::::3@@@@@@@@@@@@@@@@@@@E@EEEEEEEzzz::
+#::::::::::::::::::::::::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@EEEEEEE5ttttt
+#:::::::::::::::::::::::::::::::;g@@@@@@@@@@@@@@@@@@@@@@@@@@EEEEEEEEEEEtzt
+#::::::::::::::::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@E@@EEEEEEEEEEEE@@@
+#::::::::::::::::::::::::::g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEE3EEEE@@@@@@@
+#:::::::::::::::::::::;;g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEt33@@@@@@@@@@
+#:::::::::::::::::;g@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@E@@@@@@EEEtg@@@@@@@@@@@@
+#::::::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@EEEE@@@@@@@@@@@@@@@@@@@@@@@@
+#:::::::::::::@@@@@@@@@@@@@@@@@$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+#::::::::::;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+#
+# Copyleft ↄ⃝ 2012 Institut Curie
+# Author(s): Valentina Boeva, Alban Lermine (Institut Curie) 2012
+# Contact: [e-mail], [e-mail]
+# This software is distributed under the terms of the GNU General
+# Public License, either Version 2, June 1991 or Version 3, June 2007.
+
+#!/usr/bin/perl
+
+#filter out dulpicates from SAMPLE (optional) and create a control dataset w/o duplicates with the same number of reads as in the SAMPLE
+
+use strict;
+use warnings;
+use diagnostics;
+
+my $usage = qq{
+    $0
+
+    -----------------------------
+    mandatory parameters:
+    
+    -f CHiP_file               	
+    -c control_file
+    -t type [bam, sam, eland]
+    -o output file
+    -----------------------------
+    optional parameters:
+    
+    none			
+};
+
+if(scalar(@ARGV) == 0){
+    print $usage;
+    exit(0);
+}
+
+## mandatory arguments
+
+my $filename = "";
+my $output_fname = "";
+
+my $controlFilename = "";
+my $type = "";
+my $sampleOutput = "";
+my $samtools_bin_dir=;
+
+
+## optional arguments
+
+## parse command line arguments
+
+while(scalar(@ARGV) > 0){
+    my $this_arg = shift @ARGV;
+    if ( $this_arg eq '-h') {print "$usage\n"; exit; }
+
+    elsif ( $this_arg eq '-f') {$filename = shift @ARGV;}
+    elsif ( $this_arg eq '-c') {$controlFilename = shift @ARGV;}
+    elsif ( $this_arg eq '-t') {$type = shift @ARGV;}
+    elsif ( $this_arg eq '-o') {$output_fname = shift @ARGV;}
+    elsif ( $this_arg eq '-s') {$sampleOutput = shift @ARGV;}
+
+
+    elsif ( $this_arg =~ m/^-/ ) { print "unknown flag: $this_arg\n";}
+}
+
+if ( $filename eq ""){
+    die "you should specify chip file\n";
+}
+if( $controlFilename eq ""){
+    die "you should specify control file\n";
+}
+if( $type eq ""){
+    die "you should specify file type (bam, sam or eland)\n";
+}
+if( $output_fname eq ""){
+    die "you should specify output filename\n";
+}
+
+
+print "\n-----------------\n\n";
+
+my %hash;
+my $chipCount = 0;
+my @header;
+
+
+if ($type eq "eland") {
+    open FILE, "< $filename " || die "$filename : $!\n";    
+    while(<FILE>){
+	    my @fields = split(/\t/,$_);
+	    my $entry = $fields[6].":".$fields[7]."-".$fields[8];
+	    unless (exists($hash{$entry})) {
+		    $hash{$entry} = $_;
+		    $chipCount++;
+	    }
+    }
+} elsif ($type eq "sam") {
+    open FILE, "< $filename " || die "$filename : $!\n";    
+    while(<FILE>){
+	    if (m/^@/) {		
+		push(@header,$_);
+		next;
+	    }
+	    my @fields = split(/\t/,$_);
+	    next if (scalar(@fields)<10);
+	    my $entry = $fields[2].":".$fields[3]."-".$fields[1];
+	    unless (exists($hash{$entry})) {
+		    $hash{$entry} = $_;
+		    $chipCount++;
+	    }
+    }
+} elsif ($type eq "bam") {
+    open(FILE, "$samtools_bin_dir/samtools view -h $filename |") or die "$0: can't open ".$filename.":$!\n";
+    while(<FILE>){
+	    if (m/^@/) {		
+		push(@header,$_);
+		next;
+	    }
+	    my @fields = split(/\t/,$_);
+	    next if (scalar(@fields)<10);
+	    my $entry = $fields[2].":".$fields[3]."-".$fields[1];
+	    unless (exists($hash{$entry})) {
+		    $hash{$entry} = $_;
+		    $chipCount++;
+	    }
+    }
+}
+close FILE;
+print "ChIP: $chipCount\n";
+
+if ($sampleOutput ne "") {
+    
+    open OUT, "> $sampleOutput" || die "$sampleOutput: $!\n";
+
+    if ($type eq "bam" || $type eq "sam") { #print header
+	for my $headerLine (@header) {
+	    print OUT $headerLine;
+	}
+    }
+    for my $line (values %hash) {
+	print OUT $line;	  
+    }    
+    close OUT;
+}
+
+delete @hash{keys %hash};
+@header = ();
+
+my $controlCount = 0;
+if ($type eq "eland") {
+    open FILE, "< $controlFilename " || die "$controlFilename : $!\n";    
+    while(<FILE>){
+	    my @fields = split(/\t/,$_);
+	    my $entry = $fields[6].":".$fields[7]."-".$fields[8];
+	    unless (exists($hash{$entry})) {
+		    $hash{$entry} = $_;
+		    $controlCount++;
+	    }
+    }
+} elsif ($type eq "sam") {
+    open FILE, "< $controlFilename " || die "$controlFilename : $!\n";    
+    while(<FILE>){
+	    if (m/^@/) {		
+		push(@header,$_);
+		next;
+	    }
+	    my @fields = split(/\t/,$_);
+	    my $entry = $fields[2].":".$fields[3]."-".$fields[1];
+	    unless (exists($hash{$entry})) {
+		    $hash{$entry} = $_;
+		    $controlCount++;
+	    }
+    }
+} elsif ($type eq "bam") {
+    open(FILE, "$samtools_bin_dir/samtools view -h $controlFilename |") or die "$0: can't open ".$controlFilename.":$!\n";
+    while(<FILE>){
+	    if (m/^@/) {		
+		push(@header,$_);
+		next;
+	    }
+	    my @fields = split(/\t/,$_);
+	    my $entry = $fields[2].":".$fields[3]."-".$fields[1];
+	    unless (exists($hash{$entry})) {
+		    $hash{$entry} = $_;
+		    $controlCount++;
+	    }
+    }
+}
+close FILE;
+print "Control: $controlCount\n";
+my $prob = $chipCount/$controlCount;
+
+open OUT, "> $output_fname" || die "$output_fname: $!\n";
+
+if ($type eq "bam" || $type eq "sam") { #print header
+    for my $headerLine (@header) {
+	print OUT $headerLine;
+    }
+} 
+my $count = 0;
+
+for my $line (values %hash) {
+	my $rand = rand();
+	
+	if ($rand < $prob) {
+		print OUT $line;
+		$count ++;	
+	}
+	last if ($count == $chipCount);	   
+}
+
+
+if ($count < $chipCount) {
+
+	$prob = ($chipCount-$count)/$controlCount*1.1;
+
+	for my $line (values %hash) {
+		my $rand = rand();
+	
+		if ($rand < $prob) {
+			print OUT $line;
+			$count ++;	
+		}
+		last if ($count == $chipCount);	  
+	} 	
+}
+
+print "count = $count\n";
+close OUT;
+