Merge branch 'master' of https://github.com/friendsofstrandseq/pipeline

README.md

@@ -10,7 +10,12 @@ Preliminary SV calling using Strand-seq data - summarized in a [Snakemake](https
   
     * Install [mosaicatcher](https://github.com/friendsofstrandseq/mosaicatcher) (*currently you will need the `develop` branch*)
     * Get the R-scripts from [strandsequtils](https://github.com/friendsofstrandseq/strandsequtils)
-    * Install [Strand-Phaser](https://github.com/daewoooo/StrandPhaseR)
+    * Install BSgenome.Hsapiens.UCSC.hg38:
+      ```
+      source("https://bioconductor.org/biocLite.R")
+      biocLite('BSgenome.Hsapiens.UCSC.hg38')
+      ```
+    * [Strand-Phaser](https://github.com/daewoooo/StrandPhaseR) is installed automagically
 
   2. **Set up the configuration of the smakemake pipeline**
   

Snake.config.json

@@ -10,6 +10,7 @@
     "samtools"      : "samtools",
     "bcftools"      : "bcftools",
 
+    "snv_calls"     : "",
 
     "exclude_file"  : "/g/korbel/meiers/tools/mosaicatcher/mosaicatcher/data/exclude/hg19.exclude",
     "variable_bins" : {

Snakefile

@@ -289,13 +289,16 @@ rule prepare_strandphaser_config_per_chrom:
             print("bsGenome         = '", config["R_reference"], "'", sep = "", file = f)
 
 
-
+def locate_snv_vcf(wildcards):
+	if config["snv_calls"] == "":
+		return "snv_calls/{}.{}.vcf".format(config["sample"], wildcards.chrom)
+	else:
+		return "external_snv_calls/{}.{}.vcf".format(config["sample"], wildcards.chrom)
 
 rule run_strandphaser_per_chrom:
     input:
-        mergedbam    = "snv_calls/merged.bam",
         wcregions    = "strand_states/" + config["sample"] + ".strandphaser_input.txt",
-        snppositions = "snv_calls/" + config["sample"] + ".{chrom}.vcf",
+        snppositions = locate_snv_vcf,
         configfile   = "log/StrandPhaseR.{chrom}.config",
         strandphaser = "utils/R-packages/StrandPhaseR/R/StrandPhaseR",
         bamfolder    = "bam"
@@ -403,3 +406,15 @@ rule merge_SNV_calls:
     shell:
         config["bcftools"] + " concat -O v -o {output} {input}"
 
+rule split_external_snv_calls:
+    input:
+        vcf = config["snv_calls"]
+    output:
+        vcf = "external_snv_calls/" + config["sample"] + ".{chrom}.vcf"
+    log: "external_snv_calls/" + config["sample"] + ".{chrom}.vcf.log"
+    params:
+        bcftools = config["bcftools"],
+        #sample = config["sample"]
+    shell:
+        "({params.bcftools} view --samples " + config["sample"] + " --types snps {input.vcf} {wildcards.chrom} | bcftools view --genotype het - > {output.vcf}) > {log} 2>&1"
+