echo "MissSample . "\
"reads/chrom0.30X.1.fa.gz:reads/chrom0.30X.2.fa.gz,"\
"reads/chrom1.30X.1.fa.gz:reads/chrom1.30X.2.fa.gz" > samples.txt
../../scripts/make-pipeline.pl -r ../data/chr22/chr22_17M_18M.fa 31 proj samples.txt > task.k31.mk
kmers: 31
outdir: proj
sample_file: samples.txt
sample_names: MissSample
Usage: make -f <script> [options] [target]
  --always-run          Run/list all commands, inc. those already run
  --dry-run             List commands, don't run them
  CTXDIR=<mccortexdir>  Path to McCortex directory e.g. CTXDIR=~/mccortex
  MEM=<MEM>             Maximum memory to use e.g. MEM=80G
  NTHREADS=<N>          Maximum number of job threads to use

mkdir -p reads
mkdir -p diploid
../../libs/bioinf-perl/sim_mutations/sim_mutations.pl --snps 1000 --indels 100 --invs 0 diploid 2 ../data/chr22/chr22_17M_18M.fa
ref: 'chr22_17M_18M'
Genome size: 1,000,000
 snps: 999 / 1,000 (99.90%) generated
 insertions: 54 / 50 (108.00%) generated
 deletions: 45 / 50 (90.00%) generated
 inversions: 0 / 0 generated
cat diploid/genome0.fa | tr -d '-' | ../../libs/seq_file/bin/dnacat -u -F - > diploid/chrom0.fa
cat diploid/genome1.fa | tr -d '-' | ../../libs/seq_file/bin/dnacat -u -F - > diploid/chrom1.fa
../../libs/bioinf-perl/sim_mutations/sim_vcf.pl ../data/chr22/chr22_17M_18M.fa diploid/genome0.fa diploid/mask0.fa diploid/genome1.fa diploid/mask1.fa > truth.k31.vcf
ref: 'chr22_17M_18M'
2 Genome and mask pairs loaded
../../libs/bcftools/bcftools norm --remove-duplicates --fasta-ref ../data/chr22/chr22_17M_18M.fa --multiallelics +both truth.k31.vcf > truth.k31.norm.vcf
Lines total/modified/skipped:	1098/26/0
../../libs/htslib/bgzip -f truth.k31.norm.vcf
../../libs/bcftools/bcftools index -f truth.k31.norm.vcf.gz
../../libs/readsim/readsim -l 150 -r diploid/chrom0.fa -d 30 -e 0.01 reads/chrom0.30X
Sampling from diploid/chrom0.fa
 sequencing depth: 30.00
 read length: 150
 read pairs: yes
 insert length: 250
 insert stddev: 0.20 * insert = 50.00
 seq error rate: 1.00%
 Loaded contigs: genome0[999877]
 Genome size: 999877
Sampling 99987 paired-end reads...
Wrote 29996100 bases to: reads/chrom0.30X.1.fa.gz and reads/chrom0.30X.2.fa.gz
../../libs/readsim/readsim -l 150 -r diploid/chrom1.fa -d 30 -e 0.01 reads/chrom1.30X
Sampling from diploid/chrom1.fa
 sequencing depth: 30.00
 read length: 150
 read pairs: yes
 insert length: 250
 insert stddev: 0.20 * insert = 50.00
 seq error rate: 1.00%
 Loaded contigs: genome1[999900]
 Genome size: 999900
Sampling 99990 paired-end reads...
Wrote 29997000 bases to: reads/chrom1.30X.1.fa.gz and reads/chrom1.30X.2.fa.gz
make -f task.k31.mk CTXDIR=../.. MEM=1G bubblevcf
make[1]: Entering directory `/data1/users/turner/cortex_sims/ninja-cortex/results/bubble_calling'
mkdir -p proj/k31/graphs
mkdir -p proj/k31/links
mkdir -p proj/k31/contigs
mkdir -p proj/k31/bubbles
mkdir -p proj/k31/breakpoints
mkdir -p proj/k31/ref
mkdir -p proj/vcfs
../../bin/mccortex31 build  -m 1G -t 2 -k 31 --sample MissSample --seq2 reads/chrom0.30X.1.fa.gz:reads/chrom0.30X.2.fa.gz --seq2 reads/chrom1.30X.1.fa.gz:reads/chrom1.30X.2.fa.gz proj/k31/graphs/MissSample.raw.ctx >& proj/k31/graphs/MissSample.raw.ctx.log
../../bin/mccortex31 clean  -m 1G -t 2 --fallback 2 --covg-before proj/k31/graphs/MissSample.raw.covg.csv -o proj/k31/graphs/MissSample.clean.ctx proj/k31/graphs/MissSample.raw.ctx >& proj/k31/graphs/MissSample.clean.ctx.log
../../bin/mccortex31 inferedges  -m 1G -t 2 proj/k31/graphs/MissSample.clean.ctx >& proj/k31/graphs/MissSample.inferedges.ctx.log
../../bin/mccortex31 build  -m 1G -t 2 -k 31 --sample ref --seq ../data/chr22/chr22_17M_18M.fa proj/k31/ref/ref.ctx >& proj/k31/ref/ref.ctx.log
../../bin/mccortex31 thread  -m 1G -t 2 --seq reads/chrom0.30X.1.fa.gz --seq reads/chrom0.30X.2.fa.gz --seq reads/chrom1.30X.1.fa.gz --seq reads/chrom1.30X.2.fa.gz -o proj/k31/links/MissSample.se.raw.ctp.gz proj/k31/graphs/MissSample.clean.ctx >& proj/k31/links/MissSample.se.raw.ctp.gz.log
../../bin/mccortex31 links -L 5000 -T 0.001 proj/k31/links/MissSample.se.raw.ctp.gz > proj/k31/links/MissSample.se.thresh.txt 2> proj/k31/links/MissSample.se.thresh.txt.log
THRESH=`grep 'suggested_cutoff=' proj/k31/links/MissSample.se.thresh.txt | grep -oE '[0-9,]+$'`; \
	../../bin/mccortex31 links -c "$THRESH" -o proj/k31/links/MissSample.se.clean.ctp.gz proj/k31/links/MissSample.se.raw.ctp.gz >& proj/k31/links/MissSample.se.clean.ctp.gz.log
../../bin/mccortex31 thread  -m 1G -t 2 -p proj/k31/links/MissSample.se.clean.ctp.gz --seq2 reads/chrom0.30X.1.fa.gz:reads/chrom0.30X.2.fa.gz --seq2 reads/chrom1.30X.1.fa.gz:reads/chrom1.30X.2.fa.gz -o proj/k31/links/MissSample.pe.raw.ctp.gz proj/k31/graphs/MissSample.clean.ctx >& proj/k31/links/MissSample.pe.raw.ctp.gz.log
../../bin/mccortex31 links -L 5000 -T 0.001 proj/k31/links/MissSample.pe.raw.ctp.gz > proj/k31/links/MissSample.pe.thresh.txt 2> proj/k31/links/MissSample.pe.thresh.txt.log
THRESH=`grep 'suggested_cutoff=' proj/k31/links/MissSample.pe.thresh.txt | grep -oE '[0-9,]+$'`; \
	../../bin/mccortex31 links -c "$THRESH" -o proj/k31/links/MissSample.pe.clean.ctp.gz proj/k31/links/MissSample.pe.raw.ctp.gz >& proj/k31/links/MissSample.pe.clean.ctp.gz.log
../../bin/mccortex31 bubbles  -m 1G -t 2 --haploid 1 -o proj/k31/bubbles/bubbles.txt.gz -p 0:proj/k31/links/MissSample.pe.clean.ctp.gz proj/k31/graphs/MissSample.clean.ctx proj/k31/ref/ref.ctx >& proj/k31/bubbles/bubbles.txt.gz.log
../../scripts/cortex_print_flanks.sh proj/k31/bubbles/bubbles.txt.gz > proj/k31/bubbles/bubbles.flanks.fa.gz
../../libs/bwa/bwa index ../data/chr22/chr22_17M_18M.fa
[bwa_index] Pack FASTA... 0.01 sec
[bwa_index] Construct BWT for the packed sequence...
[bwa_index] 0.31 seconds elapse.
[bwa_index] Update BWT... 0.01 sec
[bwa_index] Pack forward-only FASTA... 0.01 sec
[bwa_index] Construct SA from BWT and Occ... 0.14 sec
[main] Version: 0.7.12-r1044
[main] CMD: ../../libs/bwa/bwa index ../data/chr22/chr22_17M_18M.fa
[main] Real time: 0.491 sec; CPU: 0.491 sec
../../libs/bwa/bwa mem ../data/chr22/chr22_17M_18M.fa proj/k31/bubbles/bubbles.flanks.fa.gz > proj/k31/bubbles/bubbles.flanks.sam
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1000 sequences (329123 bp)...
[M::mem_process_seqs] Processed 1000 reads in 0.144 CPU sec, 0.215 real sec
[main] Version: 0.7.12-r1044
[main] CMD: ../../libs/bwa/bwa mem ../data/chr22/chr22_17M_18M.fa proj/k31/bubbles/bubbles.flanks.fa.gz
[main] Real time: 0.222 sec; CPU: 0.150 sec
../../bin/mccortex31 calls2vcf -F proj/k31/bubbles/bubbles.flanks.sam -o proj/k31/bubbles/bubbles.raw.vcf proj/k31/bubbles/bubbles.txt.gz ../data/chr22/chr22_17M_18M.fa >& proj/k31/bubbles/bubbles.raw.vcf.log
../../scripts/bash/vcf-sort proj/k31/bubbles/bubbles.raw.vcf > proj/k31/bubbles/bubbles.sort.vcf
../../libs/bcftools/bcftools norm --remove-duplicates --fasta-ref ../data/chr22/chr22_17M_18M.fa --multiallelics +both proj/k31/bubbles/bubbles.sort.vcf | \
	../../scripts/bash/vcf-rename > proj/k31/bubbles/bubbles.norm.vcf
Lines total/modified/skipped:	1019/22/0
../../libs/htslib/bgzip -f proj/k31/bubbles/bubbles.norm.vcf
../../libs/bcftools/bcftools index -f proj/k31/bubbles/bubbles.norm.vcf.gz
../../libs/bcftools/bcftools concat --allow-overlaps --remove-duplicates proj/k31/bubbles/bubbles.norm.vcf.gz | \
	../../scripts/bash/vcf-rename | ../../libs/bcftools/bcftools view --output-type z --output-file proj/vcfs/bubbles.k31.vcf.gz -
../../libs/bcftools/bcftools index -f proj/vcfs/bubbles.k31.vcf.gz
make[1]: Leaving directory `/data1/users/turner/cortex_sims/ninja-cortex/results/bubble_calling'
../../libs/bcftools/bcftools isec truth.k31.norm.vcf.gz proj/vcfs/bubbles.k31.vcf.gz -p truthisec
McCortex Missed:  189 / 1098 (17.21%)
McCortex FP:       12 /  921 ( 1.30%)
McCortex Found:   909 / 1098 (82.79%)
cd cortex && make K=31
make[1]: Entering directory `/data1/users/turner/cortex_sims/ninja-cortex/results/bubble_calling/cortex'
mkdir -p ref
/apps/well/stampy/1.0.23-py2.6/stampy.py -G ref/chr22_17M_18M /data1/users/turner/cortex_sims/ninja-cortex/results/data/chr22/chr22_17M_18M.fa
stampy: Building genome...
stampy: Input files: ['/data1/users/turner/cortex_sims/ninja-cortex/results/data/chr22/chr22_17M_18M.fa']
stampy: Done
/apps/well/stampy/1.0.23-py2.6/stampy.py -g ref/chr22_17M_18M -H ref/chr22_17M_18M
stampy: Building hash table...
stampy: Initializing...
stampy: Counting...
stampy: Working... (0.0 %)         
stampy: Working... (10.0 %)         
stampy: Initializing hash...         
stampy: Flagging high counts...           
stampy: Working... (20.0 %)         
stampy: Creating hash...            
stampy: Working... (27.3 %)         
stampy: Working... (36.4 %)         
stampy: Working... (45.5 %)         
stampy: Working... (54.5 %)         
stampy: Working... (63.6 %)         
stampy: Working... (72.7 %)         
stampy: Working... (81.8 %)         
stampy: Working... (90.9 %)         
stampy: Writing...             
stampy: Finished building hash table
stampy: Done
../../..//bin/mccortex31 build -k 31 -s REF -1 /data1/users/turner/cortex_sims/ninja-cortex/results/data/chr22/chr22_17M_18M.fa ref/ref.k31.ctx >& ref/ref.k31.ctx.log
echo /data1/users/turner/cortex_sims/ninja-cortex/results/data/chr22/chr22_17M_18M.fa > ref/ref.falist
printf "MrSample\t.\t%s\t%s\n" reads.1.falist reads.2.falist > samples.txt
(echo `pwd`/../reads/chrom0.30X.1.fa.gz; \
	 echo `pwd`/../reads/chrom1.30X.1.fa.gz) > reads.1.falist
(echo `pwd`/../reads/chrom0.30X.2.fa.gz; \
	 echo `pwd`/../reads/chrom1.30X.2.fa.gz) > reads.2.falist
~/cortex/releases/CORTEX_release_v1.0.5.21/scripts/calling/run_calls.pl \
--first_kmer 31 \
--last_kmer 31 \
--kmer_step 2 \
--fastaq_index samples.txt \
--auto_cleaning yes \
--bc yes \
--pd no \
--outdir cortex_run \
--outvcf chr22_17M_18M \
--ploidy 2 \
--stampy_hash ref/chr22_17M_18M \
--stampy_bin /apps/well/stampy/1.0.23-py2.6/stampy.py \
--list_ref_fasta ref/ref.falist \
--refbindir ref/ \
--genome_size 1000000 \
--qthresh 5 \
--mem_height 20 --mem_width 100 \
--vcftools_dir ~/bioinf/vcftools_0.1.12b/ \
--do_union yes \
--ref CoordinatesAndInCalling \
--workflow independent \
--logfile runcalls.k31.log
Warning message:
In xy.coords(x, y, xlabel, ylabel, log) :
  9276 y values <= 0 omitted from logarithmic plot
sort -k1,1d -k2,2n
sort -k1,1d -k2,2n
( ../../..//scripts/bash/vcf-header cortex_run/vcfs/chr22_17M_18M_union_BC_calls_k31.decomp.vcf | \
	  grep -v '^##contig' | \
	  grep -v '^#CHROM' | \
	  sed 's/, Description=/,Description=/g'; \
	  echo '##INFO=<ID=KMER,Number=1,Type=Integer,Description="Kmer used for calling">'; \
	  echo '##contig=<ID=chr22_17M_18M,length=1000000,assembly=hg19>'; \
	  ../../..//scripts/bash/vcf-header cortex_run/vcfs/chr22_17M_18M_union_BC_calls_k31.decomp.vcf | grep '^#CHROM' ) > new_header.k31.txt
( cat new_header.k31.txt; \
	  ~/c/vcf-hack/bin/vcfref -s cortex_run/vcfs/chr22_17M_18M_union_BC_calls_k31.decomp.vcf /data1/users/turner/cortex_sims/ninja-cortex/results/data/chr22/chr22_17M_18M.fa | grep -v '^#' | sort -k1,1d -k2,2n ) > cortex.k31.sort.vcf
Loading /data1/users/turner/cortex_sims/ninja-cortex/results/data/chr22/chr22_17M_18M.fa
Loaded: 'chr22_17M_18M'
 Done.
../../..//libs/bcftools/bcftools norm --remove-duplicates --fasta-ref /data1/users/turner/cortex_sims/ninja-cortex/results/data/chr22/chr22_17M_18M.fa --multiallelics +both cortex.k31.sort.vcf > cortex.k31.norm.vcf
Lines total/modified/skipped:	906/0/0
../../..//libs/htslib/bgzip cortex.k31.norm.vcf
../../..//libs/bcftools/bcftools index cortex.k31.norm.vcf.gz
make[1]: Leaving directory `/data1/users/turner/cortex_sims/ninja-cortex/results/bubble_calling/cortex'
../../libs/bcftools/bcftools isec truth.k31.norm.vcf.gz cortex/cortex.k31.norm.vcf.gz -p truthisec2
Cortex Missed:  201 / 1098 (18.31%)
Cortex FP:        9 /  906 ( 0.99%)
Cortex Found:   897 / 1098 (81.69%)
make -f task.k31.mk CTXDIR=../.. MEM=1G breakpointvcf
make[1]: Entering directory `/data1/users/turner/cortex_sims/ninja-cortex/results/bubble_calling'
../../bin/mccortex31 breakpoints  -m 1G -t 2 -s ../data/chr22/chr22_17M_18M.fa -o proj/k31/breakpoints/breakpoints.txt.gz -p 0:proj/k31/links/MissSample.pe.clean.ctp.gz proj/k31/graphs/MissSample.clean.ctx proj/k31/ref/ref.ctx >& proj/k31/breakpoints/breakpoints.txt.gz.log
../../bin/mccortex31 calls2vcf -o proj/k31/breakpoints/breakpoints.raw.vcf proj/k31/breakpoints/breakpoints.txt.gz ../data/chr22/chr22_17M_18M.fa >& proj/k31/breakpoints/breakpoints.raw.vcf.log
../../scripts/bash/vcf-sort proj/k31/breakpoints/breakpoints.raw.vcf > proj/k31/breakpoints/breakpoints.sort.vcf
../../libs/bcftools/bcftools norm --remove-duplicates --fasta-ref ../data/chr22/chr22_17M_18M.fa --multiallelics +both proj/k31/breakpoints/breakpoints.sort.vcf | \
	../../scripts/bash/vcf-rename > proj/k31/breakpoints/breakpoints.norm.vcf
Lines total/modified/skipped:	1946/44/0
../../libs/htslib/bgzip -f proj/k31/breakpoints/breakpoints.norm.vcf
../../libs/bcftools/bcftools index -f proj/k31/breakpoints/breakpoints.norm.vcf.gz
../../libs/bcftools/bcftools concat --allow-overlaps --remove-duplicates proj/k31/breakpoints/breakpoints.norm.vcf.gz | \
	../../scripts/bash/vcf-rename | ../../libs/bcftools/bcftools view --output-type z --output-file proj/vcfs/breakpoints.k31.vcf.gz -
../../libs/bcftools/bcftools index -f proj/vcfs/breakpoints.k31.vcf.gz
make[1]: Leaving directory `/data1/users/turner/cortex_sims/ninja-cortex/results/bubble_calling'
../../libs/bcftools/bcftools isec truth.k31.norm.vcf.gz proj/vcfs/breakpoints.k31.vcf.gz -p truthisec
McCortex-brkpt Missed:  131 / 1098 (11.93%)
McCortex-brkpt FP:       21 /  988 ( 2.13%)
McCortex-brkpt Found:   967 / 1098 (88.07%)
