Replacing impute2 with Beagle; Support for hg19 and hg38 reference genomes

.gitignore

@@ -4,3 +4,4 @@ build/
 *.swp
 *.DS_Store
 #.+
+hg19_GRCh37_1000genomes

README.ACEseqWorkflow.txt

@@ -15,6 +15,10 @@ isNoControlWorkflow	        false   Run analysis with matching control and estim
 
 == Changelist
 
+* Version update 6.0.0
+- Changed the phasing routine. The program "impute2" was replaced by "Beagle". Files and tools were renamed accordingly.
+- Generally renamed all tools and files from "imputeGenotype" to "phaseGenotype" (and so on) as the subroutine does not actually perform imputation but rather phasing.
+
 * Version update to 5.0.1
 - fixed density(NA) bug and index bug for frequencies (as.character) in clustering step
 

buildversion.txt

@@ -1,2 +1,2 @@
-5.0
+6.0
 0

documentation/source/methods.rst

@@ -2,7 +2,7 @@ Methods - Theory
 ================
 
 
-| ACEseq can be used to estimate copy-numbers from WGS data using a tumor vs. control approach. Thus a pre-requesite is WGS data from healthy tissue and tumor tissue of the same patient with at least 30x coverage. Samtools [] mpileup is used to determine the coverage for tumor and control sample - position specific for each single nucleotide polymorphism (SNP) position recorded in dbSNP and per 1 kb window. To get chromosome specific allele frequencies, the genotypes of SNP positions are phased with Impute2 [] and A and B allele are assigned accordingly. Haploblocks are defined as regions with consecutively phased SNPs. Subsequently, B-allele frequencies (BAFs) are estimated for all SNP positions in tumor and control with sufficient coverage in the control: 
+| ACEseq can be used to estimate copy-numbers from WGS data using a tumor vs. control approach. Thus a pre-requesite is WGS data from healthy tissue and tumor tissue of the same patient with at least 30x coverage. Samtools [] mpileup is used to determine the coverage for tumor and control sample - position specific for each single nucleotide polymorphism (SNP) position recorded in dbSNP and per 1 kb window. To get chromosome specific allele frequencies, the genotypes of SNP positions are phased with Beagle [] and A and B allele are assigned accordingly. Haploblocks are defined as regions with consecutively phased SNPs. Subsequently, B-allele frequencies (BAFs) are estimated for all SNP positions in tumor and control with sufficient coverage in the control: 
 
   .. math::
      \label{eq:BAF}

documentation/source/parameters.rst

@@ -11,7 +11,7 @@ Multiple parameters can be set with ACEseq though not all are necessary to chang
  svOutputDirectory,${outputAnalysisBaseDirectory}/SV_calls,path,
  crestOutputDirectory,${outputAnalysisBaseDirectory}/crest,path,
  cnvSnpOutputDirectory,${aceseqOutputDirectory}/cnv_snp,path,
- imputeOutputDirectory,${aceseqOutputDirectory}/phasing,path,
+ phasingOutputDirectory,${aceseqOutputDirectory}/phasing,path,
  plotOutputDirectory,${aceseqOutputDirectory}/plots,path,
  runWithoutControl,false,boolean,use control for analysis (false|true)
  minHT,5,integer,minimum number of consecutive SNPs to be considered for haploblocks
@@ -41,7 +41,7 @@ Multiple parameters can be set with ACEseq though not all are necessary to chang
  min_distance,0.05,float,min_distance
  haplogroupFilePrefix,haploblocks_chr,string,prefix for file with haplogroups per chromosome
  haplogroupFileSuffix,txt,string,suffix for file with haplogroups per chromosome
- haplogroupFilePath,${imputeOutputDirectory}/${haplogroupFilePrefix},path,
+ haplogroupFilePath,${phasingOutputDirectory}/${haplogroupFilePrefix},path,
  min_length_purity,1000000,integer,minimum length of segments to be considered for tumor cell content and ploidy estimation
  min_hetSNPs_purity,500,integer,minimum number of control heterozygous SNPs in segments to be considered for tumor cell content and ploidy estimation
  dh_stop,max,string,
@@ -72,18 +72,9 @@ Multiple parameters can be set with ACEseq though not all are necessary to chang
  PATIENTSEX,male,string,patient sex used in case of no control workflow (male|female|klinefelter)
  CNV_ANNO_SUFFIX,cnv.anno.tab.gz,string,suffix for mappability annotated chromosome-wise 1kb coverage files
  CNV_SUFFIX,cnv.tab.gz,string,suffix chromosome-wise 1kb coverage files
- FILE_UNPHASED_PRE,${imputeOutputDirectory}/${unphasedGenotypesFilePrefix},path,
- FILE_UNPHASED_GENOTYPE,${imputeOutputDirectory}/unphased_genotype_chr,path,
- FILE_PHASED_PRE,${imputeOutputDirectory}/${phasedGenotypesFilePrefix},path,
- FILE_PHASED_GENOTYPE,${imputeOutputDirectory}/phased_genotype_chr,path,
- FILE_INFO,info,string,
- FILE_INFO_SAMPLE,info_by_sample,string,
- FILE_HAPS,haps,string,
- FILE_HAPS_CONF,haps_confidence,string,
- FILE_SUMMARY,summary,string,
- FILE_WARNINGS,warnings,string,
- FILE_PART,part,string,
- FILE_SAMPLE_G,${imputeOutputDirectory}/sample_g.txt,path,sample_g file used by imputation on X chromosome for females
+ FILE_UNPHASED_PRE,${phasingOutputDirectory}/${unphasedGenotypesFilePrefix},path,
+ FILE_PHASED_GENOTYPE,${phasingOutputDirectory}/phased_genotype_chr,path,
+ FILE_SAMPLE_G,${phasingOutputDirectory}/sample_g.txt,path,sample_g file used by imputation on X chromosome for females
  MALE_FAKE_CONTROL_PRE,${pathToACEseqResults}/cnv_snp/${pid}.chr,path,path and prefix to chromosome-wise 1kb coverage file used for fake control workflow for male patients
  FEMALE_FAKE_CONTROL_PRE,${pathToACEseqResults}/cnv_snp/${pid}.chr,path,path and prefix to chromosome-wise 1kb coverage file used for fake control workflow for female patients
  PLOT_PRE,${aceseqOutputDirectory}/${pid}_plot,path,
@@ -104,16 +95,13 @@ Multiple parameters can be set with ACEseq though not all are necessary to chang
  CHROMOSOME_LENGTH_FILE,${path}/chrlengths.txt,path,
  REPLICATION_TIME_FILE,${path}/ReplicationTime_10cellines_mean_10KB.Rda,path,"replication timing file"
  GC_CONTENT_FILE,${path}/hg19_GRch37_100genomes_gc_content_10kb.txt,path,
- GENETIC_MAP_FILE,${path}/genetic_map_chr${CHR_NAME}_combined_b37.txt,path,"impute files"
- KNOWN_HAPLOTYPES_FILE,${path}/ALL.chr${CHR_NAME}.integrated_phase1_v3. 20101123.snps_indels_svs.genotypes.nomono.haplotypes.gz,path,"impute files"
- KNOWN_HAPLOTYPES_LEGEND_FILE,${path}ALL.chr${CHR_NAME}.integrated_phase1_v3. 20101123.snps_indels_svs.genotypes.nomono.legend.gz,path,"impute files"
- GENETIC_MAP_FILE_X,${path}/genetic_map_chrX_nonPAR_combined_b37.txt,path,"impute files"
- KNOWN_HAPLOTYPES_FILE_X,${path}/ALL_1000G_phase1integrated_v3_chrX_nonPAR_impute.hap.gz,path,"impute files"
- KNOWN_HAPLOTYPES_LEGEND_FILE_X,${path}/ALL_1000G_phase1integrated_v3_chrX_nonPAR_impute.legend.gz,path,"impute files"
  outputExecutionDirectory,${path}/exec_${executionTimeString},,"path to log files"
- imputeBaseDirectory,${path}/,path,"directory for impute files"
  mergedBamSuffix,merged.mdup.bam,string,"A list of all known suffixes for merged bam files. I.e. merged.dupmark.bam, merged.mdup.bam..." 
  mergedBamSuffixList,${mergedBamSuffix},string,"A list of all known suffixes for merged bam files. I.e. merged.dupmark.bam, merged.mdup.bam..."
  defaultMergedBamSuffix,${mergedBamSuffix},string,The default suffix for merged bam files when they are created by Roddy.
  libloc_PSCBS,,string,path to PSCBS library in R
  libloc_flexclust,,string,path to felxclust library in R
+ BEAGLE_REFERENCE_FILE,${baseDirectoryReference}/tools_data/Beagle/chr${CHR_NAME}.1kg.phase3.v5a.b37.bref3,path
+ BEAGLE_REFERENCE_FILE_X,${baseDirectoryReference}/tools_data/Beagle/chrX.1kg.phase3.v5a.b37.bref3,path
+ BEAGLE_GENETIC_MAP,${baseDirectoryReference}/tools_data/genetic_maps/plink.chr${CHR_NAME}.GRCh37.map,path
+ BEAGLE_GENETIC_MAP_X,${baseDirectoryReference}/tools_data/genetic_maps/plink.chrX.GRCh37.map,path

installation/downloadReferences.sh

@@ -8,16 +8,16 @@ trap 'echo "Download incomplete. Please restart script."' ERR
 
 # convenience function to create a directory and change into it
 mkdir_cd() {
-	mkdir -p "$1"
-	cd "$1"
+    mkdir -p "$1"
+    cd "$1"
 }
 
 # compute an MD5 sum over all files found recursively in the current directory
 # and check it against the MD5 sum given in the variable EXPECTED_MD5SUM
 check_md5sum() {
-	local FILES=$(find -type f | sort)
-	local MD5SUM=$([ -n "$FILES" ] && cat $FILES | md5sum | cut -f1 -d' ')
-	[ "$EXPECTED_MD5SUM" = "$MD5SUM" ]
+    local FILES=$(find -type f | sort)
+    local MD5SUM=$([ -n "$FILES" ] && cat $FILES | md5sum | cut -f1 -d' ')
+    [ "$EXPECTED_MD5SUM" = "$MD5SUM" ]
 }
 
 
@@ -31,118 +31,128 @@ mkdir_cd hg19_GRCh37_1000genomes
 # download reference genome
 ###############################################################################
 (
-	mkdir_cd sequence/1KGRef
+    mkdir_cd sequence/1KGRef
 
-	EXPECTED_MD5SUM=12a0bed94078e2d9e8c00da793bbc84e
-	check_md5sum && exit 0 || echo downloading reference genome....
+    EXPECTED_MD5SUM=12a0bed94078e2d9e8c00da793bbc84e
+    check_md5sum && exit 0 || echo downloading reference genome....
 
-	wget -c ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz
-	gunzip hs37d5.fa.gz
+    wget -c ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz
+    gunzip hs37d5.fa.gz
 
-	check_md5sum
+    check_md5sum
 )
 
-
 ###############################################################################
-# download dbSNP database
+# download Beagle database
 ###############################################################################
 (
-	DBSNP_VERSION=135
-	mkdir_cd databases/dbSNP/dbSNP_$DBSNP_VERSION
+    mkdir_cd databases/Beagle_references
 
-	EXPECTED_MD5SUM=fed2a31b5a5d8fe12e072576c0c17199
-	check_md5sum && exit 0 || echo downloading dbSNP file....
+    EXPECTED_MD5SUM=529947b78dbedfb7a947b530da510887
+    check_md5sum && exit 0 || echo downloading Beagle reference files....
 
-	# CITATION
-	# As a NCBI Resource: "Sherry ST, Ward MH, Kholodov M, Baker J, Phan L, Smigielski EM, Sirotkin K. dbSNP: the NCBI database of genetic variation. Nucleic Acids Res. 2001 Jan 1;29(1):308-11."
-	# As a whole for a specific build (use this!) : "Database of Single Nucleotide Polymorphisms (dbSNP). Bethesda (MD): National Center for Biotechnology Information, National Library of Medicine. (dbSNP Build ID: 141 ). Available from: http://www.ncbi.nlm.nih.gov/SNP/"
-	# A single or a range of Submitted SNP (ss) or Reference SNP (rs) entries: "Database of Single Nucleotide Polymorphisms (dbSNP). Bethesda (MD): National Center for Biotechnology Information, National Library of Medicine. dbSNP accession:{ss1 or ss1 – ss100}, (dbSNP Build ID: 141). Available from: http://www.ncbi.nlm.nih.gov/SNP/"
+    for chr in {1..22}; do wget -c "http://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase3_v5a/b37.bref3/chr${chr}.1kg.phase3.v5a.b37.bref3"; done
+    wget -c "http://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase3_v5a/b37.bref3/chrX.1kg.phase3.v5a.b37.bref3"
 
-	# DOWNLOAD
-	DBSNP_BASE_URL="ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b150_GRCh37p13/VCF"
-	wget -c "$DBSNP_BASE_URL/README.txt"
-	wget -c "$DBSNP_BASE_URL/00-All.vcf.gz"
+    check_md5sum
+)
 
-	# POST PROCESSING
-	# extract SNPs from dbSNP version 135 and older
-	zcat 00-All.vcf.gz |
-	awk '/^#/{print} /VC=SNV/{ v=$8; sub(/.*dbSNPBuildID=/, "", v); sub(/;.*/, "", v); if (v~/^[0-9]+$/ && int(v)<='$DBSNP_VERSION') print }' |
-	bgzip > 00-All.SNV.vcf.gz
-	tabix -p vcf 00-All.SNV.vcf.gz
+###############################################################################
+# download genetic maps
+###############################################################################
+(
+    mkdir_cd databases/genetic_maps
+
+    EXPECTED_MD5SUM=37a6b0f29695cd87a3249065a1eae420
+    check_md5sum && exit 0 || echo downloading genetic maps....
 
-	# CLEANUP
-	rm -f 00-All.vcf.gz
+    wget -c "http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/plink.GRCh37.map.zip"
+    unzip plink.GRCh37.map.zip
+    rm plink.GRCh37.map.zip
 
-	check_md5sum
+    check_md5sum
 )
 
 ###############################################################################
-# download mappability file
+# download dbSNP database
 ###############################################################################
 (
-	mkdir_cd databases/UCSC
+    DBSNP_VERSION=135
+    mkdir_cd databases/dbSNP/dbSNP_$DBSNP_VERSION
 
-	EXPECTED_MD5SUM=3d12d0a4d7afdb52cfd10f886d48b5f0
-	check_md5sum && exit 0 || echo downloading mappability file....
+    EXPECTED_MD5SUM=fed2a31b5a5d8fe12e072576c0c17199
+    check_md5sum && exit 0 || echo downloading dbSNP file....
 
-	wget -c http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeCrgMapabilityAlign100mer.bigWig
-	bigWigToBedGraph wgEncodeCrgMapabilityAlign100mer.bigWig /dev/stdout | bgzip > wgEncodeCrgMapabilityAlign100mer_chr.bedGraph.gz
-	tabix -p bed wgEncodeCrgMapabilityAlign100mer_chr.bedGraph.gz
-	rm -f wgEncodeCrgMapabilityAlign100mer.bigWig
+    # CITATION
+    # As a NCBI Resource: "Sherry ST, Ward MH, Kholodov M, Baker J, Phan L, Smigielski EM, Sirotkin K. dbSNP: the NCBI database of genetic variation. Nucleic Acids Res. 2001 Jan 1;29(1):308-11."
+    # As a whole for a specific build (use this!) : "Database of Single Nucleotide Polymorphisms (dbSNP). Bethesda (MD): National Center for Biotechnology Information, National Library of Medicine. (dbSNP Build ID: 141 ). Available from: http://www.ncbi.nlm.nih.gov/SNP/"
+    # A single or a range of Submitted SNP (ss) or Reference SNP (rs) entries: "Database of Single Nucleotide Polymorphisms (dbSNP). Bethesda (MD): National Center for Biotechnology Information, National Library of Medicine. dbSNP accession:{ss1 or ss1 – ss100}, (dbSNP Build ID: 141). Available from: http://www.ncbi.nlm.nih.gov/SNP/"
 
-	check_md5sum
+    # DOWNLOAD
+    DBSNP_BASE_URL="ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b150_GRCh37p13/VCF"
+    wget -c "$DBSNP_BASE_URL/README.txt"
+    wget -c "$DBSNP_BASE_URL/00-All.vcf.gz"
+
+    # POST PROCESSING
+    # extract SNPs from dbSNP version 135 and older
+    zcat 00-All.vcf.gz |
+    awk '/^#/{print} /VC=SNV/{ v=$8; sub(/.*dbSNPBuildID=/, "", v); sub(/;.*/, "", v); if (v~/^[0-9]+$/ && int(v)<='$DBSNP_VERSION') print }' |
+    bgzip > 00-All.SNV.vcf.gz
+    tabix -p vcf 00-All.SNV.vcf.gz
+
+    # CLEANUP
+    rm -f 00-All.vcf.gz
+
+    check_md5sum
 )
 
 ###############################################################################
-# download replication timings
+# download mappability file
 ###############################################################################
 (
-	mkdir_cd databases/ENCODE
+    mkdir_cd databases/UCSC
 
-	EXPECTED_MD5SUM=2a63b34a737383af2a3f7eb32801a5fa
-	check_md5sum && exit 0 || echo downloading replication timing file....
+    EXPECTED_MD5SUM=3d12d0a4d7afdb52cfd10f886d48b5f0
+    check_md5sum && exit 0 || echo downloading mappability file....
 
-	wget -c https://raw.githubusercontent.com/eilslabs/ACEseqWorkflow/github/installation/ReplicationTime_10cellines_mean_10KB.Rda
+    wget -c http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeCrgMapabilityAlign100mer.bigWig
+    bigWigToBedGraph wgEncodeCrgMapabilityAlign100mer.bigWig /dev/stdout | bgzip > wgEncodeCrgMapabilityAlign100mer_chr.bedGraph.gz
+    tabix -p bed wgEncodeCrgMapabilityAlign100mer_chr.bedGraph.gz
+    rm -f wgEncodeCrgMapabilityAlign100mer.bigWig
 
-	check_md5sum
+    check_md5sum
 )
 
 ###############################################################################
-# download chromosome statistics
+# download replication timings
 ###############################################################################
 (
-	mkdir_cd stats
+    mkdir_cd databases/ENCODE
 
-	EXPECTED_MD5SUM=801bdaa8c3b0d5c18a0637b0b29fd337
-	check_md5sum && exit 0 || echo downloading stats files....
+    EXPECTED_MD5SUM=2a63b34a737383af2a3f7eb32801a5fa
+    check_md5sum && exit 0 || echo downloading replication timing file....
 
-	wget -c http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz
-	zcat chromInfo.txt.gz | grep -Pv "(_)|(chrM)" | sed -e '1i\#chrom\tsize\tfileName' > chrlengths.txt
-	rm -f chromInfo.txt.gz
+    wget -c https://raw.githubusercontent.com/eilslabs/ACEseqWorkflow/github/installation/ReplicationTime_10cellines_mean_10KB.Rda
 
-	wget -c https://raw.githubusercontent.com/eilslabs/ACEseqWorkflow/github/installation/hg19_GRch37_100genomes_gc_content_10kb.txt
-
-	check_md5sum
+    check_md5sum
 )
 
 ###############################################################################
-# download IMPUTE database
+# download chromosome statistics
 ###############################################################################
 (
-	mkdir_cd databases/1000genomes/IMPUTE
+    mkdir_cd stats
 
-	EXPECTED_MD5SUM=261a28d6b6917340cd82ada2d7185e17
-	check_md5sum && exit 0 || echo downloading impute files....
+    EXPECTED_MD5SUM=801bdaa8c3b0d5c18a0637b0b29fd337
+    check_md5sum && exit 0 || echo downloading stats files....
 
-	wget -c https://mathgen.stats.ox.ac.uk/impute/ALL.integrated_phase1_SHAPEIT_16-06-14.nomono.tgz
-	tar -xzvf ALL.integrated_phase1_SHAPEIT_16-06-14.nomono.tgz
-	rm -f ALL.integrated_phase1_SHAPEIT_16-06-14.nomono.tgz
+    wget -c http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz
+    zcat chromInfo.txt.gz | grep -Pv "(_)|(chrM)" | sed -e '1i\#chrom\tsize\tfileName' > chrlengths.txt
+    rm -f chromInfo.txt.gz
 
-	wget -c https://mathgen.stats.ox.ac.uk/impute/ALL_1000G_phase1integrated_v3_impute.tgz
-	tar -xzvf ALL_1000G_phase1integrated_v3_impute.tgz
-	rm -f ALL_1000G_phase1integrated_v3_impute.tgz
+    wget -c https://raw.githubusercontent.com/eilslabs/ACEseqWorkflow/github/installation/hg19_GRch37_100genomes_gc_content_10kb.txt
 
-	check_md5sum
+    check_md5sum
 )
 
 echo "All files downloaded successfully"

resources/analysisTools/copyNumberEstimationWorkflow/add_haplotypes.py

@@ -16,7 +16,7 @@
 
 from python_modules import Tabfile
 
-parser=argparse.ArgumentParser(description = "Add imputed haplotypes")
+parser=argparse.ArgumentParser(description = "Add phased haplotypes")
 parser.add_argument('--inputsuffix', '-i', type = str, help = "suffix of vcf file with phased genotypes")
 parser.add_argument('--inputpath', '-p', type = str, help = "path and prefix to vcf file with phased genotypes")
 parser.add_argument('--out', '-o', default=sys.stdout, help = "Output file. Default: STDOUT")

resources/analysisTools/copyNumberEstimationWorkflow/artifact.homoDels.potentialArtifacts.hg38_liftover.txt

@@ -0,0 +1,10 @@
+#chromosome	start	end	rec
+2	95574244	96074272	2906
+2	131432422	132252426	2910
+11	48878445	50210921	2909
+14	21896742	22531023	159
+15	24084780	24584944	2910
+18	14119958	14910027	2910
+19	42705844	43295853	2910
+23	90094938	90885199	2909
+23	92274954	93125000	2862

resources/analysisTools/copyNumberEstimationWorkflow/bashLib.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+
+# Copyright (c) 2017 The ACEseq workflow developers.
+# Distributed under the MIT License (license terms are at https://www.github.com/eilslabs/ACEseqWorkflow/LICENSE.txt).
+
+dieWith() {
+  local ec="${2:-$?}"
+  local msg="${1:?No error message}"
+  echo "$msg: exit code $ec" >> /dev/stderr
+  exit "$ec"
+}

resources/analysisTools/copyNumberEstimationWorkflow/beagle_create_fake_samples.py

@@ -0,0 +1,29 @@
+#!/usr/bin/python
+
+import argparse
+
+def is_comment_line(line):
+    return line.startswith("#")
+
+parser = argparse.ArgumentParser()
+parser.add_argument('--in_file', help='Input .vcf file')
+parser.add_argument('--out_file', help='Out .vcf file')
+args = parser.parse_args()
+
+vcf_infile = open( args.in_file, "r" )
+outfile    = open( args.out_file, "w" )
+
+for vcf_line in vcf_infile:
+
+    if is_comment_line(vcf_line):
+        if vcf_line.startswith("#CHROM"):
+            vcf_line = vcf_line.rstrip().split("\t")
+            vcf_line.append('sample0')
+            vcf_line = "\t".join( vcf_line )+"\n"
+
+    else:
+        vcf_line = vcf_line.rstrip().split("\t")
+        vcf_line.append(vcf_line[9])
+        vcf_line = "\t".join( vcf_line )+"\n"
+
+    outfile.write( vcf_line )