completed sequenza rule for WES and WGS, and completed second-pass ru…

…les for 2-pass freec
CCBR · Jan 12, 2019 · d01d82c · d01d82c
1 parent a5135e4
commit d01d82c
Show file tree

Hide file tree

Showing 5 changed files with 221 additions and 0 deletions.
diff --git a/Results-template/Scripts/make_freec_pass2_exome_tn_config.pl b/Results-template/Scripts/make_freec_pass2_exome_tn_config.pl
@@ -0,0 +1,68 @@
+#!/usr/bin/perl -w
+use strict;
+use List::Util 'shuffle';
+
+#INPUT
+
+#my $mergedmaf = $ARGV[1] . '_out/oncotator_out/' . $ARGV[1] . '_merged.maf'; #to fix...
+#open C, ">$mergedmaf";
+
+my $outfile = $ARGV[0] . '/freec_exome_config.txt';
+my $chrLenFile = $ARGV[1];
+my $chrFiles = $ARGV[2];
+my $tumormateFile = $ARGV[3];
+my $controlmateFile = $ARGV[4];
+my $makePileup = $ARGV[5];
+my $fastaFile = $ARGV[6];
+my $SNPfile = $ARGV[7];
+my $targets = $ARGV[8];
+my @line=();
+my $contamination='';
+my $ploidy='';
+my $rep=0;
+
+my $infile=$ARGV[9];
+open G, "<$infile";
+while (<G>){
+	chomp;
+  	last if m/^$/;
+  	@line = split;
+  	next if ($line[0] =~ m'cellularity');
+  	if ($rep==0) {
+		$contamination=(1-$line[0]);
+		$ploidy=$line[1];
+		$rep++;
+	}
+}
+
+open C, ">$outfile";
+
+print C '[general]' . "\n\n";
+
+print C "BedGraphOutput = TRUE\ndegree = 1\nforceGCcontentNormalization = 1\nminCNAlength = 3\nnoisyData = TRUE\nreadCountThreshold = 50\n";
+print C "chrLenFile = $chrLenFile\n"
+print C "ploidy = $ploidy\ncontamination=$contamination\nbreakPointThreshold = 0.8\nwindow = 0\n";
+print C "chrFiles = $chrFiles\n";
+print C "minimalSubclonePresence = 30\nprintNA = FALSE\ncontaminationAdjustment = TRUE\nmaxThreads = 24\nnumberOfProcesses = 24\n";
+print C "outputDir = $ARGV[0]\n";
+
+print C '[sample]' . "\n\n";
+
+print C "mateFile = $tumormateFile\n";
+print C "inputFormat = BAM\nmateOrientation = FR\n\n";
+
+print C '[control]' . "\n\n";
+
+print C "mateFile = $controlmateFile\n";
+print C "inputFormat = BAM\nmateOrientation = FR\n";
+
+print C '[target]' . "\n\n";
+
+print C "captureRegions = $targets\n";
+
+print C '[BAF]' . "\n\n";
+
+print C "makePileup = $makePileup\n";
+print C "fastaFile = $fastaFile\n";
+print C "minimalCoveragePerPosition = 20\nminimalQualityPerPosition = 20\n";
+print C "SNPfile = $SNPfile";
diff --git a/Results-template/Scripts/make_freec_pass2_wgs_tn_config.pl b/Results-template/Scripts/make_freec_pass2_wgs_tn_config.pl
@@ -0,0 +1,64 @@
+#!/usr/bin/perl -w
+use strict;
+use List::Util 'shuffle';
+
+#INPUT
+
+#my $mergedmaf = $ARGV[1] . '_out/oncotator_out/' . $ARGV[1] . '_merged.maf'; #to fix...
+#open C, ">$mergedmaf";
+
+my $outfile = $ARGV[0] . '/freec_wgs_config.txt';
+my $chrLenFile = $ARGV[1];
+my $chrFiles = $ARGV[2];
+my $tumormateFile = $ARGV[3];
+my $controlmateFile = $ARGV[4];
+my $makePileup = $ARGV[5];
+my $fastaFile = $ARGV[6];
+my $SNPfile = $ARGV[7];
+my @line=();
+my $contamination='';
+my $ploidy='';
+my $rep=0;
+
+my $infile=$ARGV[8];
+open G, "<$infile";
+while (<G>){
+	chomp;
+  	last if m/^$/;
+  	@line = split;
+  	next if ($line[0] =~ m'cellularity');
+  	if ($rep==0) {
+		$contamination=(1-$line[0]);
+		$ploidy=$line[1];
+		$rep++;
+	}
+}
+
+open C, ">$outfile";
+
+print C '[general]' . "\n\n";
+
+print C "chrLenFile = $chrLenFile\n"
+print C "ploidy = 2,3,4,5,6\nbreakPointThreshold = 0.8\nwindow = 1000\n";
+print C "chrFiles = $chrFiles\n";
+print C "minimalSubclonePresence = 20\ncontaminationAdjustment = TRUE\nmaxThreads = 24\nnumberOfProcesses = 24\n";
+print C "outputDir = $ARGV[0]\n";
+
+print C '[sample]' . "\n\n";
+
+print C "mateFile = $tumormateFile\n";
+print C "inputFormat = BAM\nmateOrientation = FR\n\n";
+
+print C '[control]' . "\n\n";
+
+print C "mateFile = $controlmateFile\n";
+print C "inputFormat = BAM\nmateOrientation = FR\n";
+
+print C '[target]' . "\n\n";
+
+print C '[BAF]' . "\n\n";
+
+print C "makePileup = $makePileup\n";
+print C "fastaFile = $fastaFile\n";
+print C "minimalCoveragePerPosition = 20\nminimalQualityPerPosition = 20\n";
+print C "SNPfile = $SNPfile";
diff --git a/Rules/freec_exome_somatic_pass2.rl b/Rules/freec_exome_somatic_pass2.rl
@@ -0,0 +1,27 @@
+if config['project']['annotation'] == "hg19":
+ rule freec_exome_somatic_pass2:
+     input: newnormbam=temp("freec_out/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam"),
+            newtumorbam=temp("freec_out/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam"),
+            index1=lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bai",
+            index2=lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bai",
+            fit="sequenza_out/{x}/{x}"+"_alternative_solutions.txt",
+            freectargets=config['project']['workpath']+"/freec_targets.bed",
+     output: cnvs="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_CNVs",
+             normalpileup=temp("freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam_minipileup.pileup"),
+             tumorpileup=temp("freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_minipileup.pileup"),
+     params: dir=config['project']['workpath'],fasta=config['references'][pfamily]['FREECFASTA'],tumorsample=lambda wildcards: config['project']['pairs'][wildcards.x][1],normalsample=lambda wildcards: config['project']['pairs'][wildcards.x][0],lengths=config['references'][pfamily]['FREECLENGTHS'],chroms=config['references'][pfamily]['FREECCHROMS'],pile=config['references'][pfamily]['FREECPILEUP'],snps=config['references'][pfamily]['FREECSNPS'],rname="pl:freec"
+     shell: "mkdir -p freec_out/pass2; mkdir -p freec_out/pass2/{params.normalsample}+{params.tumorsample}; perl Scripts/make_freec_pass2_exome_tn_config.pl {params.dir}/freec_out/pass2/{params.normalsample}+{params.tumorsample} {params.lengths} {params.chroms} {params.dir}/freec_out/{input.newtumorbam} {params.dir}/freec_out/{input.newnormalbam} {params.pile} {params.fasta} {params.snps} {input.freectargets} {input.fit}; module load samtools/1.9; module load freec/11.5; module load bedtools/2.27.1; freec -conf {params.dir}/freec_out/pass2/{params.normalsample}+{params.tumorsample}/freec_wgs_config.txt"
+
+else:
+ rule freec_exome_somatic_pass2:
+     input: normal=lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam",
+            tumor=lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam",
+            index1=lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bai",
+            index2=lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bai",
+            fit="sequenza_out/{x}/{x}"+"_alternative_solutions.txt",
+            freectargets=config['project']['workpath']+"/freec_targets.bed",
+     output: cnvs="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_CNVs",
+             normalpileup="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam_minipileup.pileup",
+             tumorpileup="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_minipileup.pileup",
+     params: dir=config['project']['workpath'],fasta=config['references'][pfamily]['FREECFASTA'],tumorsample=lambda wildcards: config['project']['pairs'][wildcards.x][1],normalsample=lambda wildcards: config['project']['pairs'][wildcards.x][0],lengths=config['references'][pfamily]['FREECLENGTHS'],chroms=config['references'][pfamily]['FREECCHROMS'],pile=config['references'][pfamily]['FREECPILEUP'],snps=config['references'][pfamily]['FREECSNPS'],rname="pl:freec"
+     shell: "mkdir -p freec_out; mkdir -p freec_out/pass2; mkdir -p freec_out/pass2/{params.normalsample}+{params.tumorsample}; perl Scripts/make_freec_pass2_exome_tn_config.pl {params.dir}/freec_out/pass2/{params.normalsample}+{params.tumorsample} {params.lengths} {params.chroms} {params.dir}/{input.tumor} {params.dir}/{input.normal} {params.pile} {params.fasta} {params.snps} {input.freectargets} {input.fit}; module load freec/11.5; module load samtools/1.9; module load bedtools/2.27.1; freec -conf {params.dir}/freec_out/pass2/{params.normalsample}+{params.tumorsample}/freec_wgs_config.txt"
diff --git a/Rules/freec_wgs_somatic_pass2.rl b/Rules/freec_wgs_somatic_pass2.rl
@@ -0,0 +1,41 @@
+if config['project']['annotation'] == "hg19":
+ rule freec_wgs_somatic_pass2:
+     input: newnormbam=temp("freec_out/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam"),
+            newtumorbam=temp("freec_out/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam"),
+            index1=lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bai",
+            index2=lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bai",
+            fit="sequenza_out/{x}/{x}"+"_alternative_solutions.txt",
+     output: cnvs="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_CNVs",
+             newnormbam=temp("freec_out/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam"),
+             newtumorbam=temp("freec_out/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam"),
+             normalpileup="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam_minipileup.pileup",
+             tumorpileup="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_minipileup.pileup",
+     params: dir=config['project']['workpath'],fasta=config['references'][pfamily]['FREECFASTA'],tumorsample=lambda wildcards: config['project']['pairs'][wildcards.x][1],normalsample=lambda wildcards: config['project']['pairs'][wildcards.x][0],lengths=config['references'][pfamily]['FREECLENGTHS'],chroms=config['references'][pfamily]['FREECCHROMS'],pile=config['references'][pfamily]['FREECPILEUP'],snps=config['references'][pfamily]['FREECSNPS'],rname="pl:freec"
+     shell: "module load samtools/1.9
+             mkdir -p freec_out
+             mkdir -p freec_out/pass2
+             mkdir -p freec_out/pass2/{params.normalsample}+{params.tumorsample}
+             perl Scripts/make_freec_pass2_wgs_tn_config.pl {params.dir}/freec_out/pass2/{params.normalsample}+{params.tumorsample} {params.lengths} {params.chroms} {params.dir}/freec_out/{input.newtumorbam} {params.dir}/freec_out/{input.newnormalbam} {params.pile} {params.fasta} {params.snps} {input.fit}
+             module load freec/11.5
+             module load bedtools/2.27.1
+             freec -conf {params.dir}/freec_out/pass2/{params.normalsample}+{params.tumorsample}/freec_wgs_config.txt"
+
+else:
+ rule freec_wgs_somatic_pass2:
+     input: normal=lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam",
+            tumor=lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam",
+            index1=lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bai",
+            index2=lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bai",
+            fit="sequenza_out/{x}/{x}"+"_alternative_solutions.txt",
+     output: cnvs="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_CNVs",
+             normalpileup="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam_minipileup.pileup",
+             tumorpileup="freec_out/pass2/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_minipileup.pileup",
+     params: dir=config['project']['workpath'],fasta=config['references'][pfamily]['FREECFASTA'],tumorsample=lambda wildcards: config['project']['pairs'][wildcards.x][1],normalsample=lambda wildcards: config['project']['pairs'][wildcards.x][0],lengths=config['references'][pfamily]['FREECLENGTHS'],chroms=config['references'][pfamily]['FREECCHROMS'],pile=config['references'][pfamily]['FREECPILEUP'],snps=config['references'][pfamily]['FREECSNPS'],rname="pl:freec"
+     shell: "mkdir -p freec_out
+             mkdir -p freec_out/pass2
+             mkdir -p freec_out/pass2/{params.normalsample}+{params.tumorsample}
+             perl Scripts/make_freec_pass2_wgs_tn_config.pl {params.dir}/freec_out/pass2/{params.normalsample}+{params.tumorsample} {params.lengths} {params.chroms} {params.dir}/{input.tumor} {params.dir}/{input.normal} {params.pile} {params.fasta} {params.snps} {input.fit}
+             module load freec/11.5
+             module load samtools/1.9
+             module load bedtools/2.27.1
+             freec -conf {params.dir}/freec_out/pass2/{params.normalsample}+{params.tumorsample}/freec_wgs_config.txt"
diff --git a/Rules/sequenza.rl b/Rules/sequenza.rl
@@ -0,0 +1,21 @@
+rule sequenza:
+    input: normalpileup="freec_out/pass1/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam_minipileup.pileup",
+           tumorpileup="freec_out/pass1/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_minipileup.pileup",
+    output: normalpileupgz=temp("sequenza_out/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][0]+".recal.bam_minipileup.pileup.gz"),
+            tumorpileupgz=temp("sequenza_out/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".recal.bam_minipileup.pileup.gz"),
+            segments="sequenza_out/{x}/{x}"+"_segments.txt",
+            fit="sequenza_out/{x}/{x}"+"_alternative_solutions.txt",
+            seqz=temp("sequenza_out/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".seqz.gz"),
+            binseqz=temp("sequenza_out/{x}/"+lambda wildcards: config['project']['pairs'][wildcards.x][1]+".bin100.seqz.gz"),
+    params: pair={x},dir=config['project']['workpath']/sequenza_out/{x},tumorsample=lambda wildcards: config['project']['pairs'][wildcards.x][1],normalsample=lambda wildcards: config['project']['pairs'][wildcards.x][0],rname="pl:sequenza"
+    threads: 8
+    shell: "mkdir -p sequenza_out
+            mkdir -p sequenza_out/{params.pair}
+            module load sequenza-utils/2.2.0
+            module load samtools/1.9
+            gzip -c {input.normalpileup} > {output.normalpileupgz}
+            gzip -c {input.tumorpileup} > {output.tumorpileupgz}
+            sequenza-utils bam2seqz -p -gc {params.gc} -n {output.normalpileupgz} -t {output.tumorpileupgz} | gzip > {output.seqz}
+            sequenza-utils seqz_binning -w 100 -s {output.seqz} | gzip > {output.binseqz}
+            module load R/3.5
+            Rscript Scripts/run_sequenza.R {output.binseqz} {params.dir} {threads} {params.pair}"