minor bug fixed; switch to the setting for release

snakemake/Snakefile

@@ -1,10 +1,10 @@
 import os
 import pandas as pd
 
-configfile: 'test_pe.yaml'
+# configfile: 'test_pe.yaml'
 # configfile: 'test_se.yaml'
 # configfile: "config_local.yaml"
-# configfile: "config.yaml"
+configfile: "config.yaml"
 # configfile: "config_mouse.yaml"
 
 ''' Load from config '''
@@ -47,9 +47,6 @@ THREADS = config['THREADS']
 RAND_SEED = config['RAND_SEED']
 ''''''
 
-# Prefixes for filtered VCFs
-# PREFIX_VCF_F = os.path.join(DIR, '{CHROM}_filtered')
-
 # Bowtie 2 index extensions
 IDX_ITEMS = ['1', '2', '3', '4', 'rev.1', 'rev.2']
 
@@ -61,7 +58,7 @@ DIR_MAJOR = os.path.join(DIR, 'major')
 # Prefixes and directory paths for population reference contruction and indexing
 DIR_POP_GENOME = os.path.join(DIR, 'pop_genome/')
 POP_DIRNAME = 'thrds{0}_S{1}_b{2}_ld{3}'.format(POP_THRSD, POP_STOCHASTIC, POP_BLOCK_SIZE, POP_USE_LD)
-WG_POP_GENOME_SUFFIX = EXP_LABEL + '-{POP_LEVEL}_{GROUP}_' + POP_DIRNAME
+WG_POP_GENOME_SUFFIX = EXP_LABEL + '-' + POP_LEVEL + '_{GROUP}_' + POP_DIRNAME
 DIR_POP_GENOME_BLOCK = os.path.join(DIR_POP_GENOME, POP_DIRNAME + '/')
 DIR_POP_GENOME_BLOCK_IDX = os.path.join(DIR_POP_GENOME_BLOCK, 'indexes/')
 
@@ -115,14 +112,12 @@ rule filter_vcf:
         vcf = os.path.join(DIR_VCF, VCF_PREFIX + '{CHROM}' + VCF_SUFFIX),
         chrom_map = CHROM_MAP
     output:
-        # vcf = PREFIX_VCF_F + '.vcf'
         vcf = temp(os.path.join(DIR, '{CHROM}_filtered.vcf'))
     shell:
         # Take PASS variants
         # Does not remove mnps, since they will be needed for constructing personalized reference genome, 
         # and will be removed when building major and refflow references.
         '{BCFTOOLS} view -r {wildcards.CHROM} -c 1 -f PASS {input.vcf} | {BCFTOOLS} annotate --rename-chrs {input.chrom_map} -o {output.vcf}'
-        # '{BCFTOOLS} view -r {wildcards.CHROM} -c 1 -f PASS -V mnps,other {input.vcf} | {BCFTOOLS} annotate --rename-chrs {input.chrom_map} -o {output.vcf}'
 
 rule aggregate_vcf:
     input:
@@ -131,11 +126,3 @@ rule aggregate_vcf:
         vcf = os.path.join(DIR, EXP_LABEL + '_filtered.vcf.gz')
     shell:
         '{BCFTOOLS} concat -O z -o {output.vcf} {input.vcf}'
-
-# rule prepare_chrom_genome:
-#     input:
-#         genome = GENOME
-#     output:
-#         os.path.join(DIR, 'chr{CHROM}.fa')
-#     shell:
-#         '{SAMTOOLS} faidx {input.genome} {CHR_PREFIX}{wildcards.CHROM} > {output};'

snakemake/config.yaml

@@ -1,9 +1,15 @@
 # Genomic reads to process
+ALN_MODE : 'single-end'
 READS1 : '../test/SRR622457_1-1k.fastq'
+READS2 : ''
 
 # Name of the tested sample
 INDIV : 'test'
 
+# Experiment label that will become prefixes for many files
+# For example, 'wg' (stands for whole-genome), 'chr1'
+EXP_LABEL : 'wg'
+
 # Directory where the outputs will be put
 DIR: 'run'
 
@@ -34,6 +40,10 @@ VCF_SUFFIX : '.shapeit2_integrated_snvindels_v2a_27022019.GRCh38.phased.vcf.gz'
 
 # 1KG super populations used to build second pass population genomes
 GROUP : ['EUR', 'AMR', 'EAS', 'SAS', 'AFR'] 
+POP_LEVEL : 'superpop'
+# Replace with the below if using a second-pass reference set based on 1KG populations
+# GROUP : ['CHB', 'JPT', 'CHS', 'CDX', 'KHV', 'CEU', 'TSI', 'FIN', 'GBR', 'IBS', 'YRI', 'LWK', 'GWD', 'MSL', 'ESN', 'ASW', 'ACB', 'MXL', 'PUR', 'CLM', 'PEL', 'GIH', 'PJL', 'BEB', 'STU', 'ITU']
+# POP_LEVEL : 'pop'
 
 # mapping quality cutoff to split read into committed and deferred groups
 ALN_MAPQ_THRSD : '10' 
@@ -73,9 +83,9 @@ LENGTH_MAP : '../resources/GRCh38.length_map'
 CHROM_MAP : '../resources/GRCh38.chrom_map'
 
 # Paths of software
-BCFTOOLS : 'bcftools' # there is a bug with `bcftools consensus` in v1.9, so we recommend using the latest dev version
+BCFTOOLS : 'bcftools'
 SAMTOOLS : 'samtools'
-LEVIOSAM : '../levioSAM/leviosam'
+LEVIOSAM : 'leviosam'
 PYTHON : 'python'
 DIR_SCRIPTS : '../src'
 

snakemake/shared/alignment_single_end.Snakefile

@@ -3,7 +3,6 @@
 rule align_to_major:
     input:
         reads1 = READS1,
-        reads2 = READS2,
         idx = expand(
             os.path.join(DIR, 'major/indexes/' + EXP_LABEL + '-maj.{idx}.bt2'),
             idx = IDX_ITEMS)