From d21597aaa01ac546821568515975fb77b3aa6dc1 Mon Sep 17 00:00:00 2001 From: arangrhie Date: Fri, 7 Feb 2020 10:25:48 -0500 Subject: [PATCH] clean up and move merqury --- README.md | 10 +- scripts/README.md | 232 ------------------------ scripts/_submit_cn.sh | 33 ---- scripts/_submit_cn_hap.sh | 104 ----------- scripts/_submit_meryl2_build.sh | 161 ---------------- scripts/_submit_meryl2_build_10x.sh | 105 ----------- scripts/_submit_meryl2_build_hapmers.sh | 135 -------------- scripts/_submit_split.sh | 33 ---- scripts/build/concat_splits.sh | 23 --- scripts/build/kmerHistToPloidyDepth.jar | Bin 20371 -> 0 bytes scripts/build/meryl2_count.sh | 58 ------ scripts/build/meryl2_count_10x.sh | 55 ------ scripts/build/meryl2_diff.sh | 35 ---- scripts/build/meryl2_filt.sh | 30 --- scripts/build/meryl2_hapmers.sh | 79 -------- scripts/build/meryl2_union_sum.sh | 119 ------------ scripts/build/split.sh | 56 ------ scripts/build/split_10x.sh | 56 ------ scripts/cn_hap.sh | 62 ------- scripts/eval/bedCalcN50.jar | Bin 35204 -> 0 bytes scripts/eval/spectra-cn.sh | 225 ----------------------- scripts/plot/plot_blob.R | 47 ----- scripts/plot/plot_block.sh | 58 ------ scripts/plot/plot_spectra_asm.R | 53 ------ scripts/plot/plot_spectra_cn.R | 66 ------- scripts/trio/bedMerToPhaseBlock.jar | Bin 31354 -> 0 bytes scripts/trio/exclude_reads.sh | 21 --- scripts/trio/fastaGetGaps.jar | Bin 18792 -> 0 bytes scripts/trio/hap_blob.sh | 74 -------- scripts/trio/phase_block.sh | 68 ------- scripts/trio/spectra-hap.sh | 107 ----------- 31 files changed, 2 insertions(+), 2103 deletions(-) delete mode 100644 scripts/README.md delete mode 100755 scripts/_submit_cn.sh delete mode 100755 scripts/_submit_cn_hap.sh delete mode 100755 scripts/_submit_meryl2_build.sh delete mode 100755 scripts/_submit_meryl2_build_10x.sh delete mode 100755 scripts/_submit_meryl2_build_hapmers.sh delete mode 100755 scripts/_submit_split.sh delete mode 100755 scripts/build/concat_splits.sh delete mode 100755 scripts/build/kmerHistToPloidyDepth.jar delete mode 100755 scripts/build/meryl2_count.sh delete mode 100755 scripts/build/meryl2_count_10x.sh delete mode 100755 scripts/build/meryl2_diff.sh delete mode 100755 scripts/build/meryl2_filt.sh delete mode 100755 scripts/build/meryl2_hapmers.sh delete mode 100755 scripts/build/meryl2_union_sum.sh delete mode 100755 scripts/build/split.sh delete mode 100755 scripts/build/split_10x.sh delete mode 100755 scripts/cn_hap.sh delete mode 100755 scripts/eval/bedCalcN50.jar delete mode 100755 scripts/eval/spectra-cn.sh delete mode 100755 scripts/plot/plot_blob.R delete mode 100755 scripts/plot/plot_block.sh delete mode 100755 scripts/plot/plot_spectra_asm.R delete mode 100755 scripts/plot/plot_spectra_cn.R delete mode 100755 scripts/trio/bedMerToPhaseBlock.jar delete mode 100755 scripts/trio/exclude_reads.sh delete mode 100755 scripts/trio/fastaGetGaps.jar delete mode 100755 scripts/trio/hap_blob.sh delete mode 100755 scripts/trio/phase_block.sh delete mode 100755 scripts/trio/spectra-hap.sh diff --git a/README.md b/README.md index 2499f90..32cab25 100644 --- a/README.md +++ b/README.md @@ -24,13 +24,7 @@ export PATH=/path/to/meryl/bin:$PATH This is 'sequence', a utility for working with sequence files. -# Scripts +# Evaluate assemblies with k-mers and more -Here, bash scripts are provided to give a guideline for the following use-cases: +See [Merqury](https://github.com/marbl/merqury). -1. [Build k-mer dbs in batches](https://github.com/marbl/meryl/tree/master/scripts#1-build-k-mer-dbs) -2. [Copy-number comparison (spectra-cn and spectra-hap)](https://github.com/marbl/meryl/tree/master/scripts#2-copy-number-comparison-spectra-cn) -3. [Haplotype specific mers for trio binning and making blob-plots](https://github.com/marbl/meryl/tree/master/scripts#3-haplotype-specific-mers-for-trio-binning-and-making-blob-plots) -4. [Include or exclude read pairs having specific kmers in one of the reads](https://github.com/marbl/meryl/blob/master/scripts/README.md#4-include-or-exclude-read-pairs-having-specific-kmers-in-one-of-the-reads) - -See [scripts](https://github.com/marbl/meryl/tree/master/scripts) for more details. diff --git a/scripts/README.md b/scripts/README.md deleted file mode 100644 index 0279114..0000000 --- a/scripts/README.md +++ /dev/null @@ -1,232 +0,0 @@ -# Scripts - -Here, bash scripts are provided to provide guidelines for the following purposes: - -1. Build k-mer dbs in batches -2. Copy-number comparison (spectra-cn and spectra-hap) -3. Haplotype specific mers for trio binning and making blob-plots - -All scripts are tested on slurm environment for diploid genomes up to 3.3Gb. -Modify the `_submit_*` scripts to match your system. -`cpus` and `mem` are suggestive. - - -# 1. Build k-mer dbs - -`_submit_meryl2_build.sh` is doing the following 1. Build and 2. Merge steps automatically. -`_submit_meryl2_build_10x.sh` trims the 10X barcodes from illumina reads, and runs the 1. Build 2. Merge steps. -The barcode trimming step uses `scaff_reads` from [scaff10x](https://github.com/wtsi-hpag/Scaff10X). - -## 1. Build - -Count k-mers in , could be fastq/fasta file, gzipped or not. - -``` -meryl count k=$k_size output $output.meryl -``` - -Will generate canonical k-mer counts. - -`Threads`: By default, meryl uses the maximum number of cores available, up to 64. -This can be limited by defining `threads=$cores`. -If you want to see results faster, split the file to multiple pieces and count each kmers and merge at the end. - -`Memory`: Meryl automatically detects available memory and uses that. -To prevent overhead for allocating too much memory, use something like `memory=48` which will use up to 48G. -When the k-mer space is too big and requires more memory, meryl automatically splits countig jobs in batches and merges at the end. - -Run `meryl` with no args to see more details. - -## 2. Merge - -If counted from multiple fastq files, merge them at the end. -``` -ulimit -Sn 32000 -meryl memory=$mem union-sum output $output_meryl $meryl_dbs -``` -meryl will try to make 64 x 2 files per meryl dbs to merge. - -## 3. Histogram -What’s in my meryl db? -This will give you a short summary of the k-mers, unique / distinct / total. - -``` -meryl statistics $output_meryl | head -``` - -Without the head, the histogram is printed after the summary, which could go very long. - -Get the histogram: This gives a traditional 2-column histogram, ready for genomescope or any other tool. - -``` -meryl histogram $output_meryl > $output_meryl.hist -``` - -Ploidy estimation including erroneous k-mer boundaries can be made by this script: -``` -Java -jar -Xmx256m /path/to/this/meryl-accessories/scripts/kmerHistToPloidyDepth.jar $output_meryl.hist > $output_meryl.hist.ploidy -``` - -The output is a tdf file, with -``` -Ploidy Depth Boundary -``` - -Depth is the depth where we see a ‘peak’ for ploidy > 0. -Boundaries are the dips in the histogram. - -# 2. Copy-number comparison (spectra-cn) -This script is the meryl version spectra-cn plot, introduced and implemented by [KAT](https://github.com/TGAC/KAT). -The scripts here provide a way to draw unstacked, stacked histograms by copy numbers or specific to assemblies. No dependency to KAT. - -## 1. Get kmers from the read set and assembly (assemblies) -Use the above tricks for building the meryl db for the reads. -For each assembly, -``` -meryl count memory=8 output $meryl_db -``` -Should work fine with most of the genomes. This is using 8G of memory (or lower), taking a few seconds to build the $meryl_db. - -## 2. Get assembly-specific and shared kmers -When working with 2 genome assemblies, -such as one primary and alternate haplotype assemblies from diploid genomes, -or two haplotype assemblies generated with trio-binning, -the interest will be to see how many kmers are over / under represented compared to the expected copy numbers. - -Here, we will get assembly-specific (private to one of the assemblies) kmers and the shared kmers in separate temporary meryl dbs for convenience. - -hap1 and hap2 can be set as - ->hap1=paternal - ->hap2=maternal - -Or - ->hap1=pri - ->hap2=alt - - -``` -meryl difference output $name.${hap2}_only.meryl $name.${hap2}.meryl $name.${hap1}.meryl -meryl difference output $name.${hap1}_only.meryl $name.${hap1}.meryl $name.${hap2}.meryl -meryl intersect output $name.${hap1}_shrd.meryl $name.${hap1}.meryl $name.${hap2}.meryl -meryl intersect output $name.${hap2}_shrd.meryl $name.${hap2}.meryl $name.${hap1}.meryl -``` - -This generates 4 meryl dbs; and we need one more that sums up all the kmer counts: -``` -meryl union-sum output $name.${hap1}_${hap2}_shrd.meryl $name.${hap1}_only.meryl $name.${hap1}_shrd.meryl $name.${hap2}_only.meryl $name.${hap2}_shrd.meryl -``` - -## 3. Intersect with the read meryl db -Now, intersect the above kmer set of interest with the kmers from the read set, and get the kmer counts. - -* 0-counts in the read set; i.e. only seen in the assembly. -Most likely assembly errors or kmers difficult to get from illumina reads - -``` -meryl difference output $name.only.meryl $name.${hap1}_${hap2}_shrd.meryl $name.k21.meryl -``` - -This kmer set is a useful indicator of bp errors. -The number of k-mers are associate-able with bp, so roughly the number of distinct kmers in this set are representing number of erroneous bases in the assembly. -It is difficult to distinguish kmers not sequence-able with illumina reads though, so it might be useful for comparing overal bp error rate in different assemblies from the same individual. - -``` -meryl statistics $name.only.meryl | head -``` - -Switch the two inputs now to get: - -* 0-counts in the asm; i.e. only seen in the illumina reads -``` -meryl difference output $name.k21.0.meryl $name.k21.meryl $name.${hap1}_${hap2}_shrd.meryl -``` - -* 1~4 copy kmers -``` -for i in $(seq 1 4) -do - meryl intersect output $name.k21.$i.meryl $name.k21.meryl [ equal-to $i $name.${hap1}_${hap2}_shrd.meryl ] -done -``` - -* \>4 copy kmers -``` -meryl intersect output $name.k21.gt$i.meryl $name.k21.meryl [ greater-than $i $name.${hap1}_${hap2}_shrd.meryl ] -``` - -## 4. Spectra-cn plots - -The final histogram file can be obtained with: -``` -echo -e "Copies\tkmer_multiplicity\tCount" > $name.spectra-cn.hist -meryl histogram $name.k21.0.meryl | awk '{print "illumina-only\t"$0}' >> $name.spectra-cn.hist -for i in $(seq 1 4) -do - meryl histogram $name.k21.$i.meryl | awk -v cn=$i '{print cn"\t"$0}' >> $name.spectra-cn.hist -done -meryl histogram $name.k21.gt$i.meryl | awk -v cn=">$i" '{print cn"\t"$0}' >> $name.spectra-cn.hist -``` - -## 5. Spectra-hap plots -Since we are interested in comparing the portion of kmers each haplotype (or pseudo haplotype) assemblies have, -let’s make some plots with kmers grouped by specific / shared kmers. - -``` -meryl intersect output $name.k21.$hap1.meryl $name.k21.meryl $name.${hap1}_only.meryl -meryl intersect output $name.k21.$hap2.meryl $name.k21.meryl $name.${hap2}_only.meryl -meryl intersect output $name.k21.shrd.meryl $name.k21.meryl $name.${hap1}_shrd.meryl -``` - -Again, get the histograms by groups. -``` -echo -e "kmer\tkmer_multiplicity\tCount" > $name.spectra-hap.hist -meryl histogram $name.k21.0.meryl | awk '{print "illumina-only\t"$0}' >> $name.spectra-hap.hist -meryl histogram $name.k21.$hap1.meryl | awk -v hap="$hap1-only" '{print hap"\t"$0}' >> $name.spectra-hap.hist -meryl histogram $name.k21.$hap2.meryl | awk -v hap="$hap2-only" '{print hap"\t"$0}' >> $name.spectra-hap.hist -meryl histogram $name.k21.shrd.meryl | awk -v hap="shared" '{print hap"\t"$0}' >> $name.spectra-hap.hist -``` - -`spectracn.R` is an example to generate plots for `$name.spectra-hap.hist`. More instructions will be posted later. - - -# 3. Haplotype specific mers for trio binning and making blob-plots - -`_submit_meryl2_build_hapmers.sh` does the binning. -For making blob-plots, use the temporary scripts uploaded to [here]( https://github.com/arangrhie/TrioBinning/blob/master/hap_kmer_plot.R). -(Update will be followed) - -# 4. Include or exclude read pairs having specific kmers in one of the reads -`meryl-lookup` has a utility to include / exclude reads (or read pairs) containing specific kmers. - - -One use-case is to exclude all read pairs containing unwanted haplotype reads. -In this case, -we can exclude i.e. all the paternal reads from Hi-C reads for aligning or scaffolding with maternal reads only. - -For single-end reads: -``` -# Include all reads having kmers in $meryl -meryl-lookup -include -mers $meryl -sequence $read | pigz -c > $out.fastq.gz - -# Exclude all reads having kmers in $meryl -meryl-lookup -exclude -mers $meryl -sequence $read | pigz -c > $out.fastq.gz -``` - -For paired end reads: -``` -# Include all reads having kmers in $meryl -meryl-lookup -include -mers $meryl -sequence $read1 -sequence2 $read2 -r2 $out.R2.fastq.gz | pigz -c > $out.R1.fastq.gz - -# Exclude all reads having kmers in $meryl -meryl-lookup -exclude -mers $meryl -sequence $read1 -sequence2 $read2 -r2 $out.R2.fastq.gz | pigz -c > $out.R1.fastq.gz -``` - -`exclude_reads.sh` is a wrapper script to make it easy to filter out reads in the above use-case. -``` -./exclude_reads.sh haplotype-Paternal.filt.meryl R1.fastq.gz R2.fastq.gz mygenome.mat > mygenome.mat.log -``` -will generate `mygenome.mat.R1.fastq.gz` and `mygenome.mat.R2.fastq.gz`. diff --git a/scripts/_submit_cn.sh b/scripts/_submit_cn.sh deleted file mode 100755 index 9cbc357..0000000 --- a/scripts/_submit_cn.sh +++ /dev/null @@ -1,33 +0,0 @@ -#! /bin/bash - -if [[ "$#" -lt 3 ]]; then - echo "Usage: ./_submit_cn.sh [asm2.fasta] " - exit -1 -fi - -readdb=$1 -k=$2 -asm1=$3 -asm2=$4 -out=$5 - -args="$readdb $k $asm1 $asm2 $out" -if [ -z $out ]; then - out="$asm2" - args="$readdb $k $asm1 $out" -fi - -cpus=32 -mem=48g -name=$out.spectra-cn -script="$tools/meryl/scripts/eval/spectra-cn.sh" -partition=quick -walltime=4:00:00 -path=`pwd` - -mkdir -p logs -log=logs/$name.%A.log - -echo "\ -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args diff --git a/scripts/_submit_cn_hap.sh b/scripts/_submit_cn_hap.sh deleted file mode 100755 index 67beb69..0000000 --- a/scripts/_submit_cn_hap.sh +++ /dev/null @@ -1,104 +0,0 @@ -#! /bin/bash - -if [[ "$#" -lt 5 ]]; then - echo "Usage: ./_submit_cn_hap.sh [asm2.fasta] " - exit -1 -fi - -readdb=$1 -hap1=$2 -hap2=$3 -k=$4 -asm1=$5 -asm2=$6 -out=$7 - -if [ -z $out ]; then - out=$6 - asm2="" -fi - -mkdir -p logs -# All jobs are expected to finish within 4 hours -partition=quick -walltime=4:00:00 -path=`pwd` -extra="" - -#### Get spectra-cn plots and QV stats -cpus=32 -mem=48g -name=$out.spectra-cn -script="$tools/meryl/scripts/eval/spectra-cn.sh" -args="$readdb $k $asm1 $asm2 $out" -log=logs/$name.%A.log - -echo "\ -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args > cn.jid -jid=`cat cn.jid` - -#### Get blob plots -cpus=8 -mem=10g - -script="$tools/meryl/scripts/trio/hap_blob.sh" -# ./hap_blob.sh [asm2.fasta] -args="$hap1 $hap2 $asm1 $asm2 $out" -name=$out.blob -log=logs/$name.%A.log - -echo "\ -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args > blob.jid - - -#### Get haplotype specfic spectra-cn plots -cpus=24 -mem=10g -extra="--dependency=afterok:$jid" # Re-uses asm.meryl dbs in spectra-cn.sh. - -name=$out.spectra-hap -script="$tools/meryl/scripts/trio/spectra-hap.sh" -args="$readdb $hap1 $hap2 $k $asm1 $asm2 $out" -log=logs/$name.%A.log - -echo "\ -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args > hap.jid - -#### Get phase blocks -cpus=12 -mem=72g -extra="" - -script="$tools/meryl/scripts/trio/phase_block.sh" -# ./phase_block.sh - - -if [[ "$asm2" == "" ]] ; then - # Only one assembly given. - args="$asm1 $hap1 $hap2 $out" - name=$out.phase-block -else - args="$asm1 $hap1 $hap2 $out.${asm1/.fasta/}" - name=$out.phase-block.${asm1/.fasta/} -fi -log=logs/$name.%A.log - -echo "\ -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args > block1.jid - -if [[ "$asm2" == "" ]] ; then - exit 0 -fi - -args="$asm2 $hap1 $hap2 $out.${asm2/.fasta/}" -name=$out.phase-block.${asm2/.fasta/} -log=logs/$name.%A.log - -echo "\ -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args > block2.jid - diff --git a/scripts/_submit_meryl2_build.sh b/scripts/_submit_meryl2_build.sh deleted file mode 100755 index 9091f2e..0000000 --- a/scripts/_submit_meryl2_build.sh +++ /dev/null @@ -1,161 +0,0 @@ -#!/bin/bash - -# meryl_count=/path/to/meryl/scripts/meryl_count/ -meryl_count=$tools/meryl/scripts/meryl_count - -if [ -z $1 ]; then - echo "Usage: ./_submit_meryl2_build.sh [mem=T]" - echo -e "\t: kmer size k" - echo -e "\t: ls *.fastq.gz > input.fofn. accepts fasta, fastq, gzipped or not." - echo -e "\t: Final merged meryl db will be named as .meryl" - echo -e "\t[mem=T]: Submit memory option on sbatch [DEFAULT=TRUE]. Set it to F to turn it off." - exit -1 -fi - -k=$1 -input_fofn=$2 -out_prefix=$3 -mem_opt=$4 - - -LEN=`wc -l $input_fofn | awk '{print $1}'` - -mkdir -p logs - -# Split files when >12GB -cpus=20 -if [[ "$mem_opt" = "F" ]]; then - mem="" -else - mem="--mem=4g" -fi -name=$out_prefix.split -script=$meryl_count/split.sh -partition=quick -walltime=4:00:00 -path=`pwd` -log=logs/$name.%A_%a.log - -wait_for="" -split=0 - -split_arrs="" - -for i in $(seq 1 $LEN) -do - fq=`sed -n ${i}p $input_fofn` - GB=`du -k $fq | awk '{printf "%.0f", $1/1024/1024}'` - if [[ $GB -lt 12 ]]; then - echo "$fq is $GB, less than 12GB. Skip splitting." - echo $fq >> $input_fofn.$i - else - echo "$fq is $GB, over 12GB. Will split and run meryl in parallel." - echo "Split files will be in $input_fofn.$i" - args="$input_fofn" - split_arrs="$split_arrs$i," - wait_for="${wait_for}afterok:$split_jid," - split=1 - echo - fi -done - -if [ $split -eq 1 ]; then - split_arrs=${split_arrs%,} - echo " - sbatch -D $path -J $name --array=$split_arrs --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime --error=$log --output=$log $script $args" - sbatch -D $path -J $name --array=$split_arrs --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime --error=$log --output=$log $script $args | awk '{print $NF}' > split_jid - split_jid=`cat split_jid` - wait_for="afterok:$split_jid" - cpus=2 - if [[ "$mem_opt" = "F" ]]; then - mem="" - else - mem="--mem=1g" - fi - name=$out_prefix.concat - script=$meryl_count/concat_splits.sh - args="$k $input_fofn $out_prefix" - partition=quick - walltime=10:00 - path=`pwd` - log=logs/$name.%A.log - wait_for="--dependency=${wait_for%,}" - echo " - sbatch -D $path -J $name --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime $wait_for --error=$log --output=$log $script $args" - sbatch -D $path -J $name --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime $wait_for --error=$log --output=$log $script $args - exit 0 -else - rm $input_fofn.* - wait_for="" -fi - -offset=$((LEN/1000)) -leftovers=$((LEN%1000)) - -cpus=32 # Max: 64 per each .meryl/ file writer -if [[ "$mem_opt" = "F" ]]; then - mem="" -else - mem="--mem=40g" -fi -name=$out_prefix.count -script=$meryl_count/meryl2_count.sh -partition=quick -walltime=4:00:00 -path=`pwd` -log=logs/$name.%A_%a.log - -if [ -e meryl_count.jid ]; then - echo "Removing meryl_count.jid" - cat meryl_count.jid - rm meryl_count.jid -fi - -for i in $(seq 0 $offset) -do - args="$k $input_fofn $i" - if [[ $i -eq $offset ]]; then - arr_max=$leftovers - else - arr_max=1000 - fi - echo "\ - sbatch -J $name $mem --partition=$partition --cpus-per-task=$cpus -D $path --array=1-$arr_max --time=$walltime --error=$log --output=$log $script $args" - sbatch -J $name $mem --partition=$partition --cpus-per-task=$cpus -D $path --array=1-$arr_max --time=$walltime --error=$log --output=$log $script $args | awk '{print $NF}' >> meryl_count.jid -done - -# Wait for these jobs -WAIT="afterok:"`cat meryl_count.jid | tr '\n' ',afterok:'` -WAIT=${WAIT%,} - -## Collect .meryl list -if [ -e meryl_count.meryl.list ]; then - echo "Removing meryl_count.meryl.list" - cat meryl_count.meryl.list - rm meryl_count.meryl.list -fi - -for line_num in $(seq 1 $LEN) -do - input=`sed -n ${line_num}p $input_fofn` - name=`echo $input | sed 's/.fastq.gz$//g' | sed 's/.fq.gz$//g' | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` - name=`basename $name` - echo "$name.$k.$line_num.meryl" >> meryl_count.meryl.list -done - -cpus=48 # Max: 64 per each .meryl/ file writer -if [[ "$mem_opt" = "F" ]]; then - mem="" -else - mem="--mem=32g" -fi -walltime=4:00:00 -partition=norm -name=$out_prefix.union_sum -script=$meryl_count/meryl2_union_sum.sh -log=logs/$name.%A.log -args="$k meryl_count.meryl.list $out_prefix" -echo "\ -sbatch -J $name $mem --partition=$partition --cpus-per-task=$cpus -D $path --dependency=$WAIT --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name $mem --partition=$partition --cpus-per-task=$cpus -D $path --dependency=$WAIT --time=$walltime --error=$log --output=$log $script $args | awk '{print $NF}' > meryl_union_sum.jid - diff --git a/scripts/_submit_meryl2_build_10x.sh b/scripts/_submit_meryl2_build_10x.sh deleted file mode 100755 index 05a6df8..0000000 --- a/scripts/_submit_meryl2_build_10x.sh +++ /dev/null @@ -1,105 +0,0 @@ -#!/bin/bash - -# meryl_count=/path/to/meryl/scripts/meryl_count/ -meryl_count=$tools/meryl/scripts/meryl_count - -if [ -z $1 ]; then - echo "Usage: ./_submit_meryl2_build.sh [mem=T]" - echo -e "\t: kmer size k" - echo -e "\t: Read 1. The first 23 bases will get stripped off." - echo -e "\t: Read 2. Will be processed as normal." - echo -e "\t: Final merged meryl db will be named as .meryl" - echo -e "\t[mem=T]: Submit memory option on sbatch [DEFAULT=TRUE]. Set it to F to turn it off." - exit -1 -fi - -k=$1 -R1=$2 -R2=$3 -out_prefix=$4 -mem_opt=$5 - -mkdir -p logs - -# Split files >10GB -cpus=20 -if [[ "$mem_opt" = "F" ]]; then - mem="" -else - mem="--mem=4g" -fi -name=$out_prefix.split -partition=quick -walltime=4:00:00 -path=`pwd` -log=logs/$name.%A_%a.log - -wait_for="" -split=0 - -script=$meryl_count/split_10x.sh -LEN=`wc -l $R1 | awk '{print $1}'` - -echo "R1 will be split to trim off the barcodes." -split_arrs="1-$LEN" -args="$R1" -echo " -sbatch -D $path -J $name --array=$split_arrs --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime --error=$log --output=$log $script $args" -sbatch -D $path -J $name --array=$split_arrs --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime --error=$log --output=$log $script $args | awk '{print $NF}' > split_jid -split_jid=`cat split_jid` -wait_for="${wait_for}afterok:$split_jid," - - -##### -echo "$R2 will be split if >12G" - -script=$meryl_count/split.sh -LEN2=`wc -l $R2 | awk '{print $1}'` - -split_arrs="" - -for i in $(seq 1 $LEN2) -do - fq=`sed -n ${i}p $R2` - GB=`du -k $fq | awk '{printf "%.0f", $1/1024/1024}'` - if [[ $GB -lt 12 ]]; then - echo "$fq is $GB, less than 12GB. Skip splitting." - echo $fq >> $R2.$i - else - echo "$fq is $GB, over 12GB. Will split and run meryl in parallel. Split files will be in $R2.$i" - split_arrs="$split_arrs$i," # keep the line nums $i to split - split=1 - echo - fi -done - -if [[ $split -eq 1 ]]; then - split_arrs=${split_arrs%,} - args="$R2" - echo " - sbatch -D $path -J $name --array=$split_arrs --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime --error=$log --output=$log $script $args" - sbatch -D $path -J $name --array=$split_arrs --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime --error=$log --output=$log $script $args | awk '{print $NF}' >> split_jid - split_jid=`cat split_jid | tail -n1` - wait_for="${wait_for}afterok:$split_jid," -fi - - -cpus=2 -if [[ "$mem_opt" = "F" ]]; then - mem="" -else - mem="--mem=1g" -fi -name=$out_prefix.concat -script=$meryl_count/concat_splits.sh -args="$k $R1 $out_prefix $R2" -partition=quick -walltime=10:00 -path=`pwd` -log=logs/$name.%A.log -wait_for="--dependency=${wait_for%,}" -echo "$wait_for" -echo " -sbatch -D $path -J $name --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime $wait_for --error=$log --output=$log $script $args" -sbatch -D $path -J $name --partition=$partition $mem --cpus-per-task=$cpus --time=$walltime $wait_for --error=$log --output=$log $script $args - diff --git a/scripts/_submit_meryl2_build_hapmers.sh b/scripts/_submit_meryl2_build_hapmers.sh deleted file mode 100755 index 4ca3a22..0000000 --- a/scripts/_submit_meryl2_build_hapmers.sh +++ /dev/null @@ -1,135 +0,0 @@ -#!/bin/bash - -# build=/path/to/codes/and/scripts/ -build=$tools/meryl/scripts/build - -if [ -z $1 ]; then - echo "Usage: ./_submit_meryl2_build_hapmers.sh " - echo -e "\t: kmer size k" - echo -e "\t: ls *.fastq.gz > hap1.fofn. accepts fasta, fastq, gzipped or not." - echo -e "\t: ls *.fastq.gz > hap2.fofn. accepts fasta, fastq, gzipped or not." - echo -e "\tIntermediate kmer dbs will be named after hap1 and hap2" - exit -1 -fi - -k=$1 -hap1_fofn=$2 -hap2_fofn=$3 - -hap1=${hap1_fofn/.fofn/} -hap2=${hap2_fofn/.fofn/} - -if [ -e build.jid ]; then - echo "Removing build.jid" - cat build.jid - rm build.jid -fi - -# Build for hap1 -for hap in "$hap1" "$hap2" -do - input_fofn=$hap.fofn - echo $input_fofn - - LEN=`wc -l $input_fofn | awk '{print $1}'` - offset=$((LEN/1000)) - leftovers=$((LEN%1000)) - - mkdir -p logs - - cpus=32 # Max: 64 per each .meryl/ file writer - mem=48g - name=build - script=$build/meryl2_count.sh - partition=norm - walltime=1-0 - path=`pwd` - log=logs/$name.%A_%a.log - - for i in $(seq 0 $offset) - do - args="$k $input_fofn $i" - if [[ $i -eq $offset ]]; then - arr_max=$leftovers - else - arr_max=1000 - fi - echo "\ - sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --array=1-$arr_max --time=$walltime --error=$log --output=$log $script $args" - sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --array=1-$arr_max --time=$walltime --error=$log --output=$log $script $args >> build.jid - done - - - # Wait for these jobs - WAIT="afterok:" - - job_nums=`wc -l build.jid | awk '{print $1}'` - - if [[ $job_nums -eq 1 ]]; then - jid=`cat build.jid` - WAIT=$WAIT$jid - else - - for jid in $(cat build.jid) - do - WAIT=$WAIT$jid"," - done - fi - - ## Collect .meryl list - if [ -e build.$hap.meryl.list ]; then - echo "Removing build.meryl.list" - cat build.$hap.meryl.list - rm build.$hap.meryl.list - fi - - for line_num in $(seq 1 $LEN) - do - input=`sed -n ${line_num}p $input_fofn` - name=`echo $input | sed 's/.fastq.gz$//g' | sed 's/.fq.gz$//g' | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` - name=`basename $name` - echo "$name.$k.$line_num.meryl" >> build.$hap.meryl.list - done - - cpus=48 # Max: 64 per each .meryl/ file writer - mem=72g - name=meryl_union_sum - script=$build/meryl2_union_sum.sh - args="$k build.$hap.meryl.list $hap" - log=logs/$name.%A_%a.log - - echo "\ - sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path --dependency=$WAIT --time=$walltime --error=$log --output=$log $script $args" - sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path --dependency=$WAIT --time=$walltime --error=$log --output=$log $script $args > meryl_union_sum.jid -done - -# Wait for these jobs -WAIT="afterok:" - -job_id=meryl_union_sum.jid -job_nums=`wc -l $job_id | awk '{print $1}'` - -if [[ $job_nums -eq 1 ]]; then - jid=`cat $job_jid` - WAIT=$WAIT$jid -else - -for jid in $(cat build.jid) -do - WAIT=$WAIT$jid"," -done -fi - -cpus=48 # Max: 64 per each .meryl/ file writer -mem=20g -name=meryl_diff -script=$build/meryl2_diff.sh -args="$hap1.meryl $hap2.meryl" -log=logs/$name.%A_%a.log -partition=quick - -echo "\ -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path --dependency=$WAIT --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path --dependency=$WAIT --time=$walltime --error=$log --output=$log $script $args > meryl_diff.jid - - diff --git a/scripts/_submit_split.sh b/scripts/_submit_split.sh deleted file mode 100755 index b65612c..0000000 --- a/scripts/_submit_split.sh +++ /dev/null @@ -1,33 +0,0 @@ -#! /bin/bash - -echo "Usage: ./_submit_split.sh " -echo "Split the input.fofn if >8GB" - -if [ -z $1 ]; then - echo "No input.fofn provided. Exit." - exit -1 -fi - -LEN=`wc -l $1 | awk '{print $1}'` - -cpus=6 -mem=4g -name=split -script=$tools/meryl/scripts/meryl_count/split.sh -args=$1 -partition=norm -walltime=1-0 -path=`pwd` -extra="--array=1-$LEN" - -mkdir -p logs -if [ -z $extra ]; then - log=logs/$name.%A.log -else - log=logs/$name.%A_%a.log -fi - -echo "\ -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args" -sbatch -J $name --mem=$mem --partition=$partition --cpus-per-task=$cpus -D $path $extra --time=$walltime --error=$log --output=$log $script $args - diff --git a/scripts/build/concat_splits.sh b/scripts/build/concat_splits.sh deleted file mode 100755 index 2d61c13..0000000 --- a/scripts/build/concat_splits.sh +++ /dev/null @@ -1,23 +0,0 @@ -#!/bin/bash - -echo "Usage: ./concat_splits.sh [input2.fofn]" -echo "[input2.fofn]: Only needed for 10X data" - -k=$1 -input_fofn=$2 -out_prefix=$3 -input2_fofn=$4 - -if [[ -z $input_fofn ]] ; then - echo "No given. Exit." - exit -1 -fi - -cat $input_fofn.* > $input_fofn.splits -if [[ ! -z $input2_fofn ]]; then - cat $input2_fofn.* >> $input_fofn.splits -fi - -echo " -$tools/meryl/scripts/_submit_meryl2_build.sh $k $input_fofn.splits $out_prefix" -$tools/meryl/scripts/_submit_meryl2_build.sh $k $input_fofn.splits $out_prefix diff --git a/scripts/build/kmerHistToPloidyDepth.jar b/scripts/build/kmerHistToPloidyDepth.jar deleted file mode 100755 index 6cacb89ebccaaabc24d766b75fc4f947bb9d9420..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 20371 zcmb5VW3XnymMvPlwr#7nZQHhO+qP}nwr$(CZPwf1Qy#OuC&*NTkzBJ#(V8JS~@ zImet!Rss+R=?54X*bf%yEy*AM@F4vF`XMEx$VVk6EleZDCoLu{q@YMGCHyu1jlU>cr$ff9fITI|n|{m%u_{>y?$|IC`}IqSL7 z=-KO8o6uNT8QIf_nmITs+Q?ekm>IeW8reFUQX5$6IXGmgS~w~$A%9(Ox;cR3rQrdh zwfgcS5EF7n9?Btsi3k2b;037DG^S3#d)S|bf?E_%DpzV;S9S{x4;MzN!iQ2MSI6~Q zRe$$(ajL4WUeeKelKDo`c+28&Ax<#-sRd2JHq$k;z4P_VeR6;2v)%nem;39`U#kBR zd!5=!Cc)S=OjJQflO;V_WFgIIjH%9$adW4vguaV&I*M_(obLmqDT#Nsnd%3JFE6O=?#>FJ0*$)JgR0shJ1ZXDRvG?ChgHux!n=gJcIn7 zpqbQ?$+B9*hBU-^!j~EOy}&J}We~%?QVjH6bB%!&D%$AIdj^6tm1IJ7#ns+OSOG^W z!U8s)at`3WU1Nqhp{Gku6d^($8OhPX>DH0Mk4^0;#nErj+U`ray3GDSdZk$?j*#D= zOh#RIx-jWlrtz|Drd~2W>;)splS7=tjA>IyY{*i~Rq%N?V`%LtHZBGLDfgGRfk9l^ zR>?8ijkCE>x8Y6I70}E@+$r~FM^2h?AQDuNiao_67_RvT=oHx4J;n% zT-uEXAbqLKa~k=&WTR3gqUnOF{Yp++PXvt`;_E?!Q4rSs_VjEClQiPDOHoLB+ZlE9 zLC6Aq^eae?2FTd!MML#(=U4rC^dShfOp+K2yG{L*IN(u7oIR0k3ANgNAe-Z(u@D3> z>Rnc#6oF^U!gl?|dL*hBDLi}5z#6uob5bHN=f>&{qpQscH4D|F4070#{GzF}j1%t) zL6W!|33L2jTk7K&?~*u_43Y7iVvy$6DsXgSrI?(XA|WYJmA+fwB)=y}Dj%|aH0gaZ zZ@fVsMj{BE+hgNtpCC?1X(e&u_Vwf*Hg$auHx?Rz9h}Klqqq>T&J-!ff2$q+>}kWg zc;2T*y*4C5*Q@dGxFfeo94bMO&~CbWAENoMS6IKIXPNMY*2l62b2Uh-~<^Jz^7;pr_8p}UAr zP&So7z!ej=iIp%X33Kq{74fQmTx0?}MrFyN4FZvDH%uDZpmuf_Ayv^LSk_m{Z;o6| z)4Fk^DhBD+azm59hxTfvy{w0~5NtC>>WOE$yL^*~=8=t;wd#JV+R~W2>@}LP7j$UG z6m7nC2QELEQSoyQIXBsb7?b5L2<#)FR>$zzvZ&iz~i(pIMQDY#}I0soQb}ah>@CO5srvb zoj(`movW9sOo^;V$&3oelw1~F4VMUx6I@y=NC7QM)=Uo?&LOK89ODV;Y`&%iwh!fF zlv{Wduv1xQchuIoc^`Znvy~g8Wp>~ptmN0Ra*b8*Fy!|*JNVcovc836KI1r zi^@a!lY#|xN_Hnlww~>V=prM5&H+Pn#G40znYFJJ_kAz)kMHZvCO3TRgfjIykR_a5>JP7 zMF%gGdKsjKbt4c(=f|g|DA7B$#`CHXHOQK{eXyBaSuQ{~IF#OS2i=LlaRm*R|7O#i zc;;MZQTE!StS-6dg^6%b9R_~S2WZr&XgZl4dMaOH+eD~r3UN?P8L&E*MN`ZilekVP zTGxfih|zi5m0}p~Up$M%qzxjjV06Ne%8nN9XNu_-Ud3QkZJmo{S+W;+><+DR91Lq6 zgTt&-g@`gh-=J)Hu<&QVn~n_nu$gqwi+9}WVWq;;jQo5Ti<3*-0>_%hr=2=gqVk1Pg z0^c1OxhlzCB56~es^yhj{nK6|s)~sA7mpx4kusBVk;7>78p-;R+sW+CeyV9pgM&=P zf{qhVw=ctPbR#){Nri_`663F$85!5HzJXNhx``W70M^~9JqkB8z-?ayvnnJhouRyf zX0?2mTXfIX9^w0y)fhd!-mP<+TaHUTc4smGBU zj>sb5G(fjBb)VkwDmW+{ zc5*L8IBxp(T3Y*{vn%>MkyKr|iNj5g)r7qJ(ikv>sqn7rTH|&u+ zd0JN1E}9S%Kx=jOHr+Ew6v?7}p+PRXX-j_LyhJ9}V zM^t_venSVD%XZHwe6~DJSCpsjMvPkc(O1d*BXo?Y{$!|ks01E%jZ5$px4;Y!P~_1h zp`H+~X(~xo$dUsIBHS}5)~F5Z;u5vNpD4=dFv|lMazU7~0YjC%5^@pJ%ezLEG%R&L zo$_!fR1w?NB$f35loj-$Rc!W5soY;}5SQzAD@~}@d*$2Vy)xIrraUif`uM6jM$J6Yr*HU2bn38u47|QMB9=`mS-uBh%rM>x*Fwb};8LoENR_K zOzM1(KB3F*P(QK+8hkr=b{SlU-wQ)iPjYwd(hJrP#E3CTH%GV>N<6jeF;&ka%I}C3 zpJW?-EgMMMnY|N%pFm38{HEGOzh&;N9x>8?Z!x(FO*i>p_RZWuPWS>YYt(h^{d8y^ zo)$TF6YV7QS^G6Jv8Mu#OhaBIXtNwURZ!}R0xp83cgqifE3*yq3B zK>9bhU+2%yeXndf#y}4vTeMlHueZ0ghY&^9u;G#1!R5um zewnLXA*{DEn>Ta9Q@mf1`GHJKZBx(fh z-bN@zD0~)5>YONQdDFG)Xmn%-emC;K^ZoYuF_2ORZ%8T@aPgD%ZQyv{WE!Orc;cs8 zx}O#FFRgUFC)tYxP8qW4+06O<`u+DhU+ zQDq#zTjpdlvelo`Q2MhqqwDBUJg^L4Wd2_$0Y7maq> zqob~b{Gd@B{oE3$>Bq>?($(wn!|(SxVxr=j^vQ6@NeKuVt*ASvTd4)2r&}`mPfN0p zVEJd=nK|8BD+}gaj0&Bip|FkSZi#ZE5GTh%)Ob!rN)7pi?nGq1#nV@`Br|PZ=2=Bu z!AkcgMEpgoCKb)KgmNWe>blT4r@cyITSf~^8H$3Z;k51Ah`1?Giv;|NctgUDXuPP- zL~~claQ=lzWgLunB6Gu*-C>u#VOLuMT1h=3ZQku41CjP`D8L;JkFyj_M@Sdwhf!Z5 zdAzby%?vw@&Fu|1;7eLV*o)h<1^)YN__>rh@@i&k-I>$+j}#KeVEk27=FFeGT+f^Y z6~k-dR#zd8&$Vn?hZGowrn450UZ*+dV%kXle{w)!qqXYJA1+$wkIsS2=9&Jfa4~f# z-_Z8pyt$<3Oh1;`Smt>3i#0eVQi#MKfJ#jLR&+g2z(4Qa*4|o(P~3kbG$sBjO+Ok} zLnOx`ohqM(HVH_vWhJo9Hy>K>Le6IT03TkLRU_<}^0wr1jMVPXy^15w7^m?@z!+sf zo9$(IzxW`iC2gYkX3C@RO(zAmJr}ZcDcex5CX#He-Dt8ZyF$5gPIifi)k5CJAFfl3 zb^uol8gh_fm@vPUT22p?8qL)eTY;8!>s-Qkz7PjH7h?j6H{k&TlAZmrORy5g=vEfi z?eeU2EA}a=a=xgcaE0Wr=u>yGT6fY$4zaSh7?BMH^NGHMBGY7VeyVxPpT-5_l2X19 z_jZ@@DjXa#-m4;CvQ_moHD18PM^#sS49S&w)?*F!19mm6n9h%>IJOfE8jpi=XK59% z37qPyDS|&&!lvnrk7iaamqhWZIlIp4q>xD2hlI9crs9sC6KG>BP0_|0N(uIOqT$NE zorHWXBx$n`x9| zk8~$rkE|zV{)j)jntl}l_eiY%Fw+>EDMa@M*&VvXB1+z^bd#*alTr_SP?4z#I_eM^ zaA<7O)f#m_YxShMppS6tPT50aiaoH8l1fQeoZuwNInu!mH@w;;(jOOOAm=TB<&z}< zO>bK7YRl124TI!2Kb~c$#EiTyz+Q7s{oBTf#TYR*RN6R1?|A#Q2JVAdS%Uyds%L~k zy3|pNuo554vQkf{>gssBHBI(0#DKS&yH_iK)k)0>z?I?M{MOczBXgDyz+YWbU%`P8 z0QvFb_dn^9XGO<98DWI~qDx{jH2QiDW(G9EW|l@$dKO0Z|HpHRQn9eXIzsl=_HZ%j z%3d`-u(%9lO4-$@rz%9+XY^yGqpQ5~MFhn(MZbx@5?k&tc3wuMCfb`vu%_~}r9fIp zZf{g1FIDS^Qf%5+sPGmb#|({PSkoeLe2% ze8Kcey%q&S_+gVou^Hf*tqNO{v(aW|(n2V+uh|$)@f{?ZS%85rie(rot3%0TGA=~( z4h?qlI*5~1dfw~W+0QplO@eROOLnN_(i^P(yjy&mf*)zNS5OopIrAJ%r31!)ugpBU zL~YGGu%DP9$}0Go*X=nRnTNo{u3u_Y7(G2l4yIWoAfTAtd|t}0C=lS=4I)d7B*2w^ zSY_-NOb@^zQZ&SjU1dtpoTf0Dm-{o)Vij+;@%M0Fe$rH`GZeUHL7taDo+c#r%04@Q z)$}OJ)Sj4MR$r?%od*X^2*Q9fqmW&Iujsf#Ni?-1Vj>-f#WDTQj`6)7$wL6eYv;4o zD8h-GL#-lSFOjv`r`PgmOeN?^;EG+O@>)<9i`VsGAXLa`wt|C;)*vQ&tNE?nmjStU z`_oxQWA}np3b5u&eRcB%+cSEBAHpiBm?D( zM>VVJP0n=d%*hH@%^3_;e5fkv5Xvsk2grng$V%2!l6u&_jOpm|sG3HLPj`418@vWR z#M&9lY=TXZv!FL&R+Z(X#yDRRT?Zwom{D(kENf@n-x!0As1`<-B!jsoH_&|8K?23L zfD%F%_jd0!^eI17ccglBC1kU5xL7@Bf86PbSh|*4umhEGx_6m|>EunnXeGx6h0%c1 z9Zj*Wdr+T?o*4)1<)3Cv{?fW>tDlLKd{=z!4LW1`TFq;*Mh^vr>rz=6-Qb$6xlT-z z?r@y<eq2tPI*r%%K%lUlmmcEn zv!pUL5gn&F5$!HDa$)a3AdWGVA<03EY;)$S@hs*_uR$p(u|5S%5=^wFb^|hEze-{- zKnn}X7GIE$@ydlBU4Z`vq&6bEK zJE=4ldqyGXS`5|EqhaiF}XQ{2*lM zAsXSDy(PTkB{+F*1A$!65&7U;t84M-j~)ARqvYv;ZdGk)2Oej{2wX!xw2|Yp=>(4` zn2Ns9MpIWjl8qd{h+tArsJKx`?%P%RAk^31$B z1X~%z%sE&yZ?NjOu-Msab81m}qtOa-<6z_Mdr>F8b@O}c9cb$niR+buM;C_cMkVPR zN@x<`oN;Cy(n)={%%9FXK!{V4GbS%a;~QT8cOp;BeHX$`@ln()x5;iTNFH9^*OJ9% z+iIT5x5zn$?AgU31Frle_ch(S%ssW7<{Q89*1i^ci-*9*&LFy(OfYHyl*HP12`eQT zOpOATG7TmfOiIZr6SVv$!g{y=f(7BoCMRe5cW{qQ{%>P(8@#EO_~S0Z<`lx~0ql?A zm{?LwP|*P^O_CM|GVH;MLW3}=TO61Jg2#iO(YvsVpYQ>XT!B8fi0?F;Ri@ni2o1g5 zQz@S5G;WL&3aHtH*~*@kV6b%Y60WJ;{o%K<_?|I@Pn7uFW%>*UQeagS$Vc!(zZIyD zSjA=RUoG0(RqhylvMn_FB5=$1Ql3J|1V?ItN07e)Ev^n~!8+xi?L8lkB-bWRJ}I;FnkWr9XiPdJT$Q*FY+OSS zynBThR}Y*>U8|OEX06NxO;ROcROHqp)ta0U7ONQ+uSvF}amxocI+f(|*h9Xh+GjGp zA4;EO*R-uP#&gA--c8S2K~`oN%VqjFO^p2Rg}%r>XyLgqa&w0Sv!~SD-IUD#}{Eu+|=R@%4 z*FRQ|{(n^b@AQ@Omdye$ln>1O7idBV)zx^EiXt=|W)gAKlishu!3iR#4UGlDcFhds z(r=FL7|&{C0nA@;@{5O--Fd&P4EsKj#HPEAHN<0@D!8^nV&x(4xQ#(#<%`_QHDlA#|H% zrHIuXL*z*T_Uo_8Zj+P6qe}*E5@f}XV0JGOa=A~xeR6l@!m>lz7KHGZvX!!#uvA-W$!siZmXp!LwwDz24MQLM}i zPvUPUu1ej-e5G$K^F?(i7c7Jg6D7_0mG1V6*5^RDfMC(guaP|tw|uC!5@^eg5`XIiNMD#0Xoi4D3S7)%et^#3;%=50Yn)YjD-h#q3PUHtxXp z^j8Rnhb+$6h|6nUH|r`evmO(wn=R8-9+B3iGCm>9&0fT)=v)MitEUxtV znAjzFxAYK@z2adl9R)_58WX^E0Jg629*p3g3jJ^RT>KSbCK?_NhrDu=1~x3lRD-zO zYf>(#iwhNu?YZd2;S9gbaEkt2J~_y}n2_UB`ju$?x>nCQFQ(d$(J| zkVb(!xSE~8y^|JG*rZ|PH#@_}%{s+hSo^SjMi(lL2jB2exfGmU!60R($}$J(8JE9! zVja?4{Qmjl$05`|;^|+b;(v8({uxg;hGxcQ{}C7mDQc-;Dj{vL;6gzGjy-Av&H z2=BL|s9Ij=zYOUy4-VAd-9<&4xq!1G!BsmUWaBetXV_H!Tmr1#(pJ6pE}u*Cf`Etn zYV8*Cl^m(z902V1HR$%TSj=mO3u*(4ZR%;mhFX31B= z)cBG{LN@!%<{%J=dK}VeXpI}6=MV=eYQk7yy-vd|N_7LKa}oZ8kY5F=(?Mn#pRz!Y z)gJby?f_uTHPWesrs#oIr^#+H?+tLds_1~B9kLO5k`(zzcM$OA{;yGi6b4i=!xCm_ z$~y3=Dt)2$T#k_FN}y^Ll+*Dxqc-PSo!GBnEz&daN34!Uv-fIu8#IExY7L_chi?xk z)@pxD>XR=9x&rolzbnrUa}^9S_RM${t^hf}m$;xzjN{U0G*yduf{ty)z`Y6aF=sIP z*imi9XoxGnpw8uN1Qewub5HKN+OUwzvIdKH3FxC_rPt)3nZ*3Lk#5Y>RT@gquo^Ut z_@gXc72NvDFeOieK4s1h@+X|SN@$%0ctVnTac(zJ+dbwPHwEy-0xuB|3X>cFh8yhH zXhjgeYh$kjbB^c z4km}luk1bpP50nAg2gCvSz3SlKvkr#di=_sW9ip8Zf>;LJ zxiP=r&d?9v;bI<$qK|2<^^hhWhcDXh0Ip4gBhm|WHK8b8V=r&!*CS@h`~*M(`O@4p zBI!sYriJ0%s>@|9-pGV4T3rG~rQzj?C_XO&*`4ZHtqwMs^1xNu$IandD_jDv(KiJ~ z*soa-K8CB{tJ7}#;n^?LxV>M!R0oCehe%3%NF8G6XNJvJOZOCx(yR6!NJ@He^PHm8 zy|hReLoiJOiOV2k^?j=I0xo`4E2OXywoc!1(9TdqWu%|_U+W!9mO(+-GU|f}nWOA4 z!HQ|5V!*OwACI{MZ+OZVMdcpJ2nMZfyNRpyCq4vaSD2GJi&kk{hgbMV*B4gfJ|dBI zQ$=~w1Vu2vf?b}IzS#Vw6wJ%>xG%sTKal^(^nW#Yyrxn7lSC8$j}rYKl=4@ly%grG z_i5ogr3*g>!c!n+GcOpp%0#K;i#i0@-vEc{!fL{>*(CDVpLT4luzmq3g=39mXJut= z877D!3;=D$6WI+6VbTzO9Nm2MwjQF1Ew#>5C?l6SGKBW$u2e7|PBcN28XlHnxBM6`M#*%|K zDv+Rhp+$;s7wBQ5+vyv<7BP^y1NSc=#p^Poa-`3H5{alco>yOvm9BgYgDs~02?2~d zfd;u)k{ZzO=!|vTwm52 z>sti$2`y2YSWE**3hNrk$;_C?bh^EA#*P^Z8%VbH^9_cd0Q;jnmTD?_gv@;9g{=#R zcjnrRA3$;o)?)Xywoq2k@~hk+@EgFjH6e-W9TZmm*5oh73gb-hbV3k^foQ>Zu(lDj_Nw5 z&E-1NV`h!v;o;>4_lLb}UhhaO_uDAMya>r4on_amN6hYK4j*kau9gDPOD}FN z|02)#A$|c`AE`j8CX1DIY9YB-lHQfChR|Id*<}5w$<8R>8R;-x->k+Fbce3;;Mg;c zKz&BF|KLvOK5;Qxf0S|A9u{pOeGtFOQj6%d2K{O12Izbz449e@9c6HL{?r66^Y1tl zv&;Dk=mPR)S1D0*CHWaKyRFsJnio#mfO?usCu@b|QO7=VRb4cP@n+V#rhGY$a)`BB zxfb(H7`?Es+r(~{jEp34AdKq?mb+aH6K5n&!&4WSbjOW!EaX1={`w~|B)7_0hYg2d zhDle3^;v5d>7SaNRHM0epYBaSQAih36_>3-DreknZfkyIu_LiMzVz=ooyRhalOlFxI15Xzok4@BgGxOK(M_~vK0~1x zInL*R_1srWC)1@9H?&dk5t-yJa{HKhPQj@tp+fZ~njqy{&hCtp^24?xXZN4{>U{Kk zX2nz&41I+6!FI0wN%cn(^spsYPu!d2{AJDZ{+_M&KY`fWKLO!iS@YNb1qgrnGD}{{qMsLzOXQXZimj^Aa--pq*dO$Q z&mD0*RdU|yfcR3nQvPkH0~fq0Uo~Su+SO)zVse5xwaY;4Hdtd*jImULO=FZvJt-V za*b5WA5X-pJa3P zw_td6LB=U;8x|RiGiw^uX72#gnWn*w6X+Kx`V)1EKY7?^6QUQX$O=-pZ0ZN@!Xv;4 zH2zk{i4phA1EcgD>@1{_C!#hs57Qg^c-GF}GY~5??PmT%Mv{Nj>%Wxie@#ZhHuhF} zj{izq@Jgs_k%sidF6U%6}>K7QpZfH4;@K>)(6zrBjnjp5S8(3XO| znh?oy2mU$>Owlo0lXyUc>DT4jcb#q>fK4T1>(j-&oYx)#GsROm*V7UV;L%eh@*#NH z{n?t^N-w76xhsc3&oykDHH^k8@##!I_PDYOAEmJx$aY>~0c!yq>bM@N2WQarD&C%h z(cA>8hqvW2!U!jybu{nOsIY^X3~rhtw->E#paHfKzqh!#d*0>PFY~j&TBZtGF?)D( zgve$dD6BQvxxdP{^ZgTO~(R z5zA>2$_w8FU>KCD$D~^0IE-_NrYpk zxe*Df#`>Fc!U6o*3@DiD>WgK}AI|vtGyZFH|9Yc|{ySFrTlA1O(lh+8lMoawU9#WM-}mhor=QQ;Kzm^8r8>8- z&0a+7=DVRi!0=HeQ%l3k%hmhWMGnsw5N&YIM=4@!n++C*h~5Bm<(91@UMB6mF+Hb6 z@liaZ{QUC($?UNB9%BV1;j_^^e5u6h&n!A3QdVR3&P+y1iFJ4YooU-=Aa;GF#d?E` zW$Ll>2bewXR(yBm)it@@VUc1Zc7>lSAvOB}cKSKf#f8D{jpwGqm~t?n2i}fgjvfgh zvqV3juTG5vw}I?=MGzUk+Vcite1&yz&reObaXEg?VlYs#>5LK%fR zN`1M&dQE4ip`t~8%HR+Hc;f%MHfxRI@P^V7H6PQ$9srXEIppgTJ)tqmO?(lx9f*PV zd!}%+>E{akMC%&$$aFCj9w}*+L1v)Zr<6-yvt1lV#-gW;k`xpGGsw)GA!QoKChuc4 zxyXwKVFbAaF7k4G2@LoxF9H+;^mAVPXhpKR)8uUQuFO6vv_layjof7o(N7z;!|0B| z@5_9W#F+X+pP*f};B>^nFPk3lq)WI|4`T}>TVHa(z*`c3e9v2oN)vj8N|AQ@EVRxu zn$*HT51QH|ki(jfx2CB?FFPu*c5v2+Ts|rNDBKmD9&f&ezO@+N7$kAPHN_WL})&jyw-rahAQC96; zb-+dRRS)dZ?b+5*FUh`|APi~9Py#?A6 zsl;65qX9$+SKlw2eGu+^k_jZ}PiS%juW`yL&Kky*)g8m+CK!R9yq1FlhfQ`NRCMU7 z@9+*GPS%0fnN`L9nHf5}@OSZ7oCHUn;&-yGg^NnrK}q-qp-$5M3jA)lQ9c!q?Il&p z{UQV)eS#BJHxcb0=;ia|MIi_1`H4STMk_XuM&SHg1CZ7IEFP2kz+h}s$eGj;ck(cI z0FHFAF$Y`A5H2oqT#rV|V_>?6k}<7Qw=L`awtg4I4+%T5R;xb6`%jSHw*U0r^#_Fd z!rn-^Lu7j81*Z9i>~TQnISXUqx;r>o?TL(W-)G2v&HmHtc;}G)3b;z`Jv`KsczGf- zx9g$tt=#nNYi>|IBB8HZu4itS#+;PBJ1`*r&ee%F|FSep9XC{Fw8FW?US=!VhSfe{ zXUD6)xM>AHW078Wg+blb{}eF)b*J-}v*5~z1oi$zLBQbu$XR)<&VRoe_|NBmEY$ux zXZ>6DsI+0dK?~*O?E3+@M?!L%2ulj$IIxVzyMkg70wADcm=V3i7EdyEW;EIp6MwmC zN9{&yDbLAC2%9wI$(_}UwP`??waEECpaM4IIiQ2SqFKFqA6>K=dTGy|COz1$yzb6} z817{q{*&pLC-$V22KHkmzyDyls1oyCm+Ag#&(H6H)ZRd72#22f$Zn;hbggA);{B6w zVJyFD*$;!3$v7V_Re&$B^B|n1SrLP~Z9J>!LLxyKKcI4d*tY=Jw2vnQqe=3YlVSdl zeh)=enkiVDfBMF}S9^ei2nl;%2ylwTalkKB42FtFOj&kOcO6(|deM4-!Y zUSa5@#5m?#xd$kL2CrS`%9Ig=JkxS=FahLUVdh1nJ*5ILI46oiim2wf`H||uV8?@i zpviJ#7l_rf3;aNd?%_IX^}5PgZG|_FpQpv!`7Z8fg7h%&^F-)FR5#x3Lq#-HwgI1X zJA)3P)#n=d?l`tC=8>qr{+hk3TeP6ZgwZA>%F8U}I~qM<0d%F0Y2#?nz4JfF@!|$` z>Dh%Zs+ZVP+v>Jm;Q3ASiZB~U-%8^#e<>Lg0jmQyBJ=%V?PQXkjM~4){3&KUpj?zL z=GoL|SG~j{NRzp4!Q>RBh8%tMR*2 zqx8TG_FBsQ3-m8sI=$J~xBkQ~f2$(@y~I=f0nT4I{Kxk+|LNfJI~f}r*&7-DtJFs+ zTPPr!AbD%A*4n87x0&&3l;?{f={LYmplG5<#X-VFG z52Nvx6+tYSe6%q6Ao_%swmI=?00%E*pl`UmcAH`Q^R`*yd_Qd5{BS-%MMj3&pvGtz z&YY=@Tuq$8eUw5IV>skUg$ZdJBJr34M}{~#0{6R zhfXS;lZ#TfKYNC%C^ax8zgY|e0WnpLAlA#P`3anBfXsNeSSb*sP*4KR&J5ub4(6Lj ztUK)L!_MAUl`t$A7!kOSM?i0>TI|{qByPSgN$>hoM$EhgF#}z-BbeDzp{^-#Pji)h z^A-;ALw-G!FoiCTg1^N?ap>7|oFY6oKdUc^=V1il45{r{n9-#Zf}68O--PTr$uL~R z>9XINybhgUlNTL!3&UAA=2v}!sbh-L9hAlGK2cs{N}jjaSaB}{W$MlH9&ZxjdubrH zf(kGmsAFfO3R)&7a~h{rI5qhONa6=xaMb;;EP8_O3gdnk$3-4lq@ZvDsEhm2{2xLg z7%&B51ALcryCW6k(F&OotRlz!oN@FAbl3oM2u$GpL<$6q4)gTapz>&=Taj0fcBIao7wgOc?t7+3;FBARTLP~{_bx`0upP@_n;{S;w$tq=K7;Xb^PbO!@t&j*)t@ozkaZ zdhzbdqJ*W+P#oFhKRnKZLTaPQpb*5#vpERFi=zy6Hurxrkezrc=1!l6_8QuTx`vjg zmg;co@XUqx@iU31v#M$qp34$uH=4-1N}Y>_$#WEDPA6o}JPub)1Ph%pwtarpf0zO1 zi0vJPM_0H;HS;gl8lN6NM`*w7wMVbqfgJTxRj?NOY3huX+3izdsE3m>PLXJY4K2I2 zF74}A8sKf`?SX8|*)@>P;agze@SDgQ6$tB-Z%?rJfNvQK7HdDp@9Cdv2X2Agjr@(i{DAA~7 zpS1ZQ)1!8P@#4Ok8}U|d62Gq}LT#h&P|GmLdqQu}veOM<>p8Pgtx>iw>K?Y*V07I3 z>F9nCrM72iHE?DUI~V$8PDo2`Luh~_HV%A<+M}^um--s|!f8iR~p0Fd&BC zL?bMPk-A10dafUbL5X`55yujPe#yEtvg<63|fQ3tZA3+HqPK_T3?F$3BA3d=R7 z>i(Ko-5&wY%iQgdq9LKD)_#-ZWR?lqPnrH#Nqd5Qt|?ZRNW(F-jhovFhBuU}F&{af zr}pp9_0I;t@0f#eqpK7i1Q8GU3^=0|6ajZokqd76L9Bq@D$gnmNnA5F)mG2(jTwhr zgc~7Jtg?e8E5hEzfiiEd#2BIT5@2oQy}&y)nw$acZWCQj@~f1U|Vm zxzPlMn-WqUm0(TV!dg^vh!Uj$$aH@0WwULY!c|jw>dw=1nha&L$$Zaol{pA7^MOz( ze9P;xOA~l88?k1S@%u&y;Vr2@OfbiqQV;Li-@N^H_FH57G>uB?eWg2Kv|v&WU@O~V zESzkNZ=ZiM)9f#cz7urj42@8WfaDn7x`F6C?I@}i)bFkM;GRJBX!|Y~wIYcZeO7B? z3!qUS&XYOd=CwFHG71M~mB^u$KvBycpg+Eo1_}_u!tn@z9YA@l4VTK}?-(Q|snY!R zqK&qpqglfkad}k4>++V7QW2|s+y#v%JM_ZzCQJ+PWFr&^H!MY0@+r_3!hDCkIC}=F zCb4}DW1JgqfpiLAcQUWNAbeIat136M84|~MKoByI;}mNWmy2!?dNQv)s8Q3d@$S&( z8FPF4{ckRp?z*Rl{GTQ~w|^=(QYfMXg9l%4UNVj!U9PKNm5{a4lcO z#&<}c1rl(pSX7>}mPH(6&C*a_LSA~l|GR4dRxSAeY4uQlo^U5)ByB{F=G>!>|CPLX zKMs;$F`TBxji|B`)s32^i<$b;zU=pnz(#)GF^Eurv0bMZ5r<~Ew{6RG?keYoj`sI& zRfC&C<(U`aO|s{{v1p$JuJ4PTqZ75rPTX?oWm%Vgbu?smvN+K0lQm5iMZs)%?Q-PD z{;D&bY@+lNpvzsuvo-pbzv8~=#=uuy#`CrRkc7efAdS8`6{(Ol^tz# zqDq`1?n>+=PC?3vw)mA{4?hn1lk9GkYHh%c4hY&U7&awd_}EJ?6TzrFMjYhvhwIBr z%9&lL>t(R3S*N_nT#ehJAi#%4#vpA(R1BP7J?M@Basch|8jA1~N?XuTFJc7H%hDMo z5uwWZ0XV;{Z|V$7k=rT4Kht){aKhN#%w2RR!0$&yc}6H}dR|HRzas~!q)YghwGAIs zJ_UdJh&S^v9)uItrR#2~ZVwnKBOU1{+vX@^KInT?8d22wf!d$w%{x6M4bCiFk#1!v ziZZn@WtrzPxg(|oIqBWt|7rv{k|5szZ$v>E6?jUwl?LPrQ?-PW*Dp+9!k2xBLsLCl zc<84JXnB^&fU%o@& z0S>Ql9GNOO`=b3_IM_O-eJCSEJ$H|L zYAR<0h&Z(AkZ|PFCCzes;^ac;nJOHt%x52IDkU^u8I~D0_qcyU$=)?0biqH~XW&0& zUjOPA{kO#!c_R}cSKEL0KVv0i`gxIpW^OBif$qXW$EQ^;`-_GP8PP<;(Zp>BjMYh4 z3sDwO2Ksr9#s`CVdU1_9sfeAA^rt6XjyET}9-cn$z;}O6amqtIk@@UFM!j!lm~M3& zKk%}z^GdEqoN&@ZC6VT~%{lF4Ly#EC-f(Q`ZL}{$F$*kxK&s}^bd3=SjHSBqjds}% zZneqQ_Q&FYmJ|wf73Gfhj*eq|YRjm|5_1wXN^YH#mLixVS)Bm_XF!$RQ3za%yVG|P zzn{7e#?l$X@(W>f2wCc79fSf1-N z89VI2v(T55QC>Zi0!hh2KKiIFnjVGT3zCrzn)^cz76H{Ds33F-cV_QGv)7p-&12g) zrk-UbxGZk(i)JyRW|&>Zp`-vX*2-`U@MW@p{9A#{#6v+v=#M?Q_hb6vE^mWsQ@agIGbk@%fXTjzPAtGjJ=rKSWGa zfHi=NQd9MqQ%l3Td7 z{hK)%Fn#J4^=IKo|FpG#y>9sb-PTm>^=xhb!GlMqjTj+cM+<(b0=$@Q3_MzzLd7&*y1zA)qq${5p}Wail*{ z&W^OBVZKONq7SZ-T`=KVYw`|7Q1o3#L0k9*SJ%}Z>|Lk?rKzXseFD0Gz)#nqGy zw+R_mvk^_80B7LNTTW(4(=P|-U6&K*uOA;P;)FID%Y*8#_lwYY9oj(wBpI-eswK9U zbzrq3`2;=^UKYh*j5VE|`8l%S%wT{i+Bh0IXN7TdSPR(&jh!Yn8ADr31ceRcsn=fJ z(;y{}rX!g`>xs;?H3S^_1sy~{q8Ap>A?v+c>7)XBPu(wCEOjwy>vpgFlN%#s1sfEOodt99`djD!^%S4`2>D;>1-k*#hB+9Kn$9jkt7nzmk3KvbFKKv{Vr( zy&Ud@%_SX5z3or((QQ6qHd_!ZBYx`|JWM8RrPpX$EIKS!!bXsmc!gilX!d)iF+L}^ zr$hmhUCcQ+4nnB(up&1?`Yny6eGXCoDVtzQL^SM!?5%p&v{9lA-wB++F52z^>%b!w ze3dx#L*%p*SKvBf&Ks(O0-X#sQBoam=$PXE`O5zdfS}9ZN(X);tSCugll*3oVXwGQ zsVxWJ$Diz;z3&6Kh)$dgca<;jEp!x44!vaL(9D&E&KF`HZ;5TRZ9BTuSr;0%EOb$q za{%9rF5EMaXdWp%i`S1a$s2G8E+!^We9&`My$Xha|Qng(Jv`;01DD%Y~8YFe~-E#FRB(FB(rf!nSDN;r?;!l*~pvF zK(V1FMe%mZ80RY`35j`4)ndj*2I9r&>9y}=S}vj)J&=cd(K8NA>>;k2wfs>TQ)F3e z@v?lqcK>tdIR%bBzygvUwtPEW0&x6foIFWqX7A^?|8IRce=fZ|czpXE=in_13yt%g zcRX9y>w8Rih1vr3kMqiZuG{h7&yerqqj3MlH6o|f(iHxI3z_&=vN@{20>%Lt?uhae zy^xW`4D+CJMIvAM(O2xRiK4E1zR4?_Ks!&ElZIT{E4-w-Y(Uy-V9p%)2Y3b*LkI zhSL-KY5uTv!(P2|?@*}Q6DhmdOgvtwC?Yy7=V&&XVH zw0%WJ#N-~s^L%r(%nH04&rO&sFt_cthoMGnrjhQ<55Bz7JG_<7<@}zT)B5yL##guC zx6`${e@n^zSei9!+wbLEJpCpiwY$#?O_}h4*IVhtG0EgCs{g#`l z0sU|0T^3=?`BdJY>oLKjv86?EKgZe??U^AN%P*hFn5V@r@Y|_dAVg~MzCgz&H(d*d z;z?<3-X6UA^=by3yF%Y7*Y>^fOl`Qcdi@OB$zT5+@v_-(&GgLj_J7q| zDt}DRzIeP#IIKy+WQ~9`Yv!bB3%AetYWjqE-a@vWw~i$TO+2Di_JKvQ>Qac=iRFK% zT=AZJzNYOic%-}BGhj&|FkOnFln39@R!Yg>Nt>XepmxG}uOkKm$IeH;u;!`G68yw; zboGm)Q}42#zg!aac$w#nl>Z$+-l=Vr5}jf&edhb$wZ?@S{1$>Xk2&?3B+sh6GmDMA zI(63ErysJFs=j!2NiWUxwh$=)&n+)G_w*{@64L8lOyBGbPB8>6{=UJszN0~}@yvrn z!@bw8W);p^BbBh}tlkDAVU5XR8H$==b_(r#XFXn!x%-M|s9DRzhf+Lq7k-)5I^pSV zarq3fBRVlh*mGxnWsdf|w|s*kS0ZotL}fLuljZ%!rn#Keo4yGwRGuCt;+^|9Jx=mT z$NSg&j%rP>`|E30$@t$wZT3H-lwZ1;({GlFK70=!ok)bdml<>-5dusEGC_e48&Ab|J`uX9$UD(MD{B#;5ZEmy z!c>qFWJ5uv8S=t-(A)?Da0;Uu3Lp9f4f!H%X+t&q=-uYl?}9ZejjAumipcO3Ex0aV8+YU4K!o8yq@ z?m(*w5a1t9)6mu!AUhB=>4!YOfts*48Q^jt!boh6L>`0!&5IzwLt`8!;+iEvb}VQV z3%Lgk8pT3@J}csl#^z|`ejr-3+u<-7xj&BVP*CR)x$mWba3zp+z#hMG*c^x4um$Ei zBrw|nhiNFiH)Pj=GBt9e4b`(f&IHWE<~-zj7Sz&10 [offset line_num]" - echo -e "\t: k-size mers will be collected. REQUIRED." - echo -e "\t: list of fastq / fasta file. REQUIRED." - echo -e "\t[offset]: OPTIONAL. DEFAULT=0. For array job limit only." - echo -e "\t[line_num]: OPTIONAL. (offset * 1000 + line_num)'th line of input.fofn will be the input." - echo -e "\t\t\$SLURM_ARRAY_TASK_ID will be used if not specified." - exit -1 -fi - -k=$1 -input_fofn=$2 - - -if [ -z $3 ]; then - offset=0 -else - offset=$3 -fi - -if [ ! -z $4 ]; then - line_num=$4 -else - line_num=$SLURM_ARRAY_TASK_ID -fi - -# Note: Provide memory in Gb unit. SLURM provides $SLURM_MEM_PER_NODE in Mb. -# Give extra 4Gb to avoid 'Bus Error' form running out of memory. -if [[ -z $SLURM_MEM_PER_NODE ]]; then - mem=32 -else - mem=$(((SLURM_MEM_PER_NODE/1024)-4)) -fi -line_num=$(((offset * 1000) + $line_num)) - -# Read in the input path -input=`sed -n ${line_num}p $input_fofn` - -# Name it accordingly -name=`echo $input | sed 's/.fastq.gz$//g' | sed 's/.fq.gz$//g' | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` -name=`basename $name` - -output=$name.$k.$line_num.meryl - -if [ ! -d $output ]; then - # Run meryl count: Collect k-mer frequencies - echo " - meryl k=$k memory=$mem count $input output $output - " - meryl k=$k memory=$mem count $input output $output -else - echo "$output dir already exist. Nothing to do with $name." -fi - - diff --git a/scripts/build/meryl2_count_10x.sh b/scripts/build/meryl2_count_10x.sh deleted file mode 100755 index 94d0f4a..0000000 --- a/scripts/build/meryl2_count_10x.sh +++ /dev/null @@ -1,55 +0,0 @@ -#!/bin/bash - - -if [[ -z $1 ]] || [[ -z $2 ]]; then - echo "Usage: ./meryl2_count.sh [offset line_num]" - echo -e "\t: k-size mers will be collected. REQUIRED." - echo -e "\t: list of fastq.gz file. Read pair 1 from 10X reads. REQUIRED." - echo -e "\t[offset]: OPTIONAL. DEFAULT=0. For array job limit only." - echo -e "\t[line_num]: OPTIONAL. (offset * 1000 + line_num)'th line of input.fofn will be the input." - echo -e "\t\t\$SLURM_ARRAY_TASK_ID will be used if not specified." - echo -e "\t*NOTE* This script is trimming off the first 23 bases before kmer counting. Only useful for 10X barcode trimming." - exit -1 -fi - -k=$1 -input_fofn=$2 - - -if [ -z $3 ]; then - offset=0 -else - offset=$3 -fi - -if [ ! -z $4 ]; then - line_num=$4 -else - line_num=$SLURM_ARRAY_TASK_ID -fi - -# Note: Provide memory in Gb unit. SLURM provides $SLURM_MEM_PER_NODE in Mb. -# Give extra 4Gb to avoid 'Bus Error' form running out of memory. -mem=$(((SLURM_MEM_PER_NODE/1024)-4)) -line_num=$(((offset * 1000) + $line_num)) - -# Read in the input path -input=`sed -n ${line_num}p $input_fofn` - -# Name it accordingly -name=`echo $input | sed 's/.fastq.gz$//g' | sed 's/.fq.gz$//g' | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` -name=`basename $name` - -output=$name.$k.$line_num.meryl - -if [ ! -d $output ]; then - # Run meryl count: Collect k-mer frequencies - echo " - zcat $input | awk '{if (NR%2==1) {print $1} else {print substr($1,24)}}' | meryl k=$k threads=$SLURM_CPUS_PER_TASK memory=$mem count output $output - - " - zcat $input | awk '{if (NR%2==1) {print $1} else {print substr($1,24)}}' | meryl k=$k threads=$SLURM_CPUS_PER_TASK memory=$mem count output $output - -else - echo "$output dir already exist. Nothing to do with $name." -fi - - diff --git a/scripts/build/meryl2_diff.sh b/scripts/build/meryl2_diff.sh deleted file mode 100755 index 9d0f33a..0000000 --- a/scripts/build/meryl2_diff.sh +++ /dev/null @@ -1,35 +0,0 @@ -#!/bin/bash - -if [[ -z $1 ]] || [[ -z $2 ]]; then - echo "Usage: meryl2_diff.sh " - exit -1 -fi - -db1=$1 -db2=$2 - -db1=${db1/.meryl/} -db2=${db2/.meryl/} - -if [ ! -d ${db1}_only.meryl ]; then - echo "\ - meryl difference output ${db1}_only.meryl $db1.meryl $db2.meryl" - meryl difference output ${db1}_only.meryl $db1.meryl $db2.meryl - echo -fi - -if [ ! -d ${db2}_only.meryl ]; then - echo "\ - meryl difference output ${db2}_only.meryl $db2.meryl $db1.meryl" - meryl difference output ${db2}_only.meryl $db2.meryl $db1.meryl - echo -fi - -echo " -sh $tools/TrioBinning/meryl_count/meryl2_filt.sh ${db1}_only.meryl" -sh $tools/TrioBinning/meryl_count/meryl2_filt.sh ${db1}_only.meryl - -echo " -sh $tools/TrioBinning/meryl_count/meryl2_filt.sh ${db2}_only.meryl" -sh $tools/TrioBinning/meryl_count/meryl2_filt.sh ${db2}_only.meryl - diff --git a/scripts/build/meryl2_filt.sh b/scripts/build/meryl2_filt.sh deleted file mode 100755 index b8ef081..0000000 --- a/scripts/build/meryl2_filt.sh +++ /dev/null @@ -1,30 +0,0 @@ -#!/bin/bash - -if [ -z $1 ]; then - echo "Usage: ./meryl2_filt.sh " - exit -1 -fi - -db=$1 -db=${db/.meryl/} - -echo "Generate $db.hist" -meryl histogram $db.meryl > $db.hist - -echo " -java -jar -Xmx1g $tools/TrioBinning/meryl_count/kmerHistToPloidyDepth.jar $db.hist -" -java -jar -Xmx1g $tools/TrioBinning/meryl_count/kmerHistToPloidyDepth.jar $db.hist > $db.hist.ploidy - -cat $db.hist.ploidy - -filt=`sed -n 2p $db.hist.ploidy | awk '{print $NF}'` - -echo "Filter out kmers <= $filt" - -echo " -meryl greater-than $filt output $db.gt$filt.meryl $db.meryl -" -meryl greater-than $filt output $db.gt$filt.meryl $db.meryl - - diff --git a/scripts/build/meryl2_hapmers.sh b/scripts/build/meryl2_hapmers.sh deleted file mode 100755 index 728e84a..0000000 --- a/scripts/build/meryl2_hapmers.sh +++ /dev/null @@ -1,79 +0,0 @@ -#!/bin/sh - -if [[ -z $1 ]] || [[ -z $2 ]] || [[ -z $3 ]] || [[ -z $4 ]]; then - echo "Usage: ./meryl_count_asm.sh " - exit -1 -fi - -hap_a=$1 # meryl hap_a.mcdat -hap_b=$2 # meryl hap_b.mcdat - -a_fa=$3 # hap_a.fasta -b_fa=$4 # hap_b.fasta - -hap_a=${hap_a/.mcdat/} -hap_b=${hap_b/.mcdat/} - -a_name=`echo $a_fa | sed 's/.fasta$//g' | sed 's/.fa$//g'` -b_name=`echo $b_fa | sed 's/.fasta$//g' | sed 's/.fa$//g'` - -echo "Generate qr dbs" -if ! [ -s $hap_a.qr ]; then -echo "\ -simple-dump -m 21 -s $hap_a -e $hap_a.qr" -simple-dump -m 21 -s $hap_a -e $hap_a.qr -echo -fi - -if ! [ -s $hap_b.qr ]; then -echo "\ -simple-dump -m 21 -s $hap_b -e $hap_b.qr" -simple-dump -m 21 -s $hap_b -e $hap_b.qr -echo -fi - -echo "Start counting" -if ! [ -s ${a_name}"_"$hap_a.counts ]; then -echo "\ -simple-dump -m 21 -e $hap_a.qr -f $a_fa | awk -v name=$a_name \'{print name\\t\$0}\' > ${a_name}_$hap_a.counts" -simple-dump -m 21 -e $hap_a.qr -f $a_fa | awk -v name=$a_name '{print name"\t"$0}' > ${a_name}"_"$hap_a.counts -echo -fi - -if ! [ -s ${b_name}"_"$hap_a.counts ]; then -echo "\ -simple-dump -m 21 -e $hap_a.qr -f $b_fa | awk -v name=$b_name \'{print name\\t\$0}\' > ${b_name}_$hap_a.counts" -simple-dump -m 21 -e $hap_a.qr -f $b_fa | awk -v name=$b_name '{print name"\t"$0}' > ${b_name}"_"$hap_a.counts -echo -fi - -if ! [ -s ${a_name}"_"$hap_b.counts ]; then -echo "\ -simple-dump -m 21 -e $hap_b.qr -f $a_fa | awk -v name=$a_name \'{print name\\t\$0}\' > ${a_name}"_"$hap_b.counts" -simple-dump -m 21 -e $hap_b.qr -f $a_fa | awk -v name=$a_name '{print name"\t"$0}' > ${a_name}"_"$hap_b.counts -echo -fi - -if ! [ -s ${b_name}"_"$hap_b.counts ]; then -echo "\ -simple-dump -m 21 -e $hap_b.qr -f $b_fa | awk -v name=$b_name \'{print name\\t\$0}\' > ${b_name}"_"$hap_b.counts" -simple-dump -m 21 -e $hap_b.qr -f $b_fa | awk -v name=$b_name '{print name"\t"$0}' > ${b_name}"_"$hap_b.counts -echo -fi - -echo "Merge" -echo "\ -paste ${a_name}_$hap_a.counts ${a_name}_$hap_b.counts > $a_name.counts" -paste ${a_name}"_"$hap_a.counts ${a_name}"_"$hap_b.counts > $a_name.counts - -echo "\ -paste ${b_name}_$hap_a.counts ${b_name}_$hap_b.counts > $b_name.counts" -paste ${b_name}"_"$hap_a.counts ${b_name}"_"$hap_b.counts > $b_name.counts - -echo "\ -echo -e "Assembly\tContig\t$hap_a\t$hap_b\tTotal" > hapmers.counts" -echo -e "Assembly\tContig\t$hap_a\t$hap_b\tTotal" > hapmers.counts - -echo "\ -cat $a_name.counts $b_name.counts | awk '{print \$1\t\$2\t\$5\t\$NF\t\$3}' >> hapmers.counts" -cat $a_name.counts $b_name.counts | awk '{print $1"\t"$2"\t"$5"\t"$NF"\t"$3}' >> hapmers.counts diff --git a/scripts/build/meryl2_union_sum.sh b/scripts/build/meryl2_union_sum.sh deleted file mode 100755 index 2a97d9f..0000000 --- a/scripts/build/meryl2_union_sum.sh +++ /dev/null @@ -1,119 +0,0 @@ -#!/bin/bash -if [ -z $1 ]; then - echo "Usage: ./_meryl_union_sum.sh " - echo -e "Merge " - exit -1 -fi - -k=$1 -input_fofn=$2 -output_prefix=$3.k$k - -## Collect .meryl files -LEN=`wc -l $input_fofn | awk '{print $1}'` -NUM_DBS_TO_JOIN=100 # Join every $NUM_DBS_TO_JOIN as intermediates, then merge at the end -JOIN_IDX=0 - -echo "Set ulimit: ulimit -Sn 32000" -ulimit -Sn 32000 - -for FROM_IDX in $(seq 1 $NUM_DBS_TO_JOIN $LEN) -do - END_IDX=$((FROM_IDX+$NUM_DBS_TO_JOIN-1)) - if [ $END_IDX -gt $LEN ]; then - END_IDX=$LEN - fi - meryl="" - for i in $(seq $FROM_IDX $END_IDX) - do - input=`sed -n ${i}p $input_fofn` - if [ -d $input ]; then - meryl="$meryl $input" - fi - done - - echo "union-sum of $FROM_IDX - $END_IDX :" - - JOIN_IDX=$((JOIN_IDX+1)) - output=${output_prefix}.$JOIN_IDX - if [ ! -d $output ]; then - echo " - meryl \ - threads=$SLURM_CPUS_PER_TASK \ - memory=$((SLURM_MEM_PER_NODE/1024)) \ - k=$k \ - union-sum \ - output $output \ - $meryl - " - - meryl \ - threads=$SLURM_CPUS_PER_TASK \ - memory=$((SLURM_MEM_PER_NODE/1024)) \ - k=$k \ - union-sum \ - output $output \ - $meryl - fi -done - -if [ $JOIN_IDX -gt $NUM_DBS_TO_JOIN ]; then - echo -e "\tMore than $NUM_DBS_TO_JOIN intermediate files made." - echo -e "\tRe run meryl_union_sum.sh on $output_prefix.*.meryl." - exit 0 -fi - -if [ $JOIN_IDX -eq 1 ]; then - echo "All inputs merged to $output_prefix.1. Renaming to $output_prefix.meryl" - mv $output_prefix.1 $output_prefix.meryl - echo "Done!" - echo " - meryl histogram $output_prefix.meryl > $output_prefix.hist - " - meryl histogram $output_prefix.meryl > $output_prefix.hist - - echo "Use $output_prefix.hist for genomescope etc." - - echo "Cleaning up" - rm -r $meryl - exit 0 -fi - -meryl="" -for i in $(seq 1 $JOIN_IDX) -do - if [ -d $input ]; then - meryl="$meryl $output_prefix.$i" - fi -done - -echo "union-sum of $output_prefix.[ 1 - $JOIN_IDX ] :" - -if [ ! -d $output_prefix ]; then - echo " - meryl \ - threads=$SLURM_CPUS_PER_TASK \ - memory=$((SLURM_MEM_PER_NODE/1024)) \ - k=$k \ - union-sum \ - output $output_prefix.meryl \ - $meryl - " - - meryl \ - threads=$SLURM_CPUS_PER_TASK \ - memory=$((SLURM_MEM_PER_NODE/1024)) \ - k=$k \ - union-sum \ - output $output_prefix.meryl \ - $meryl -fi - -echo " -meryl histogram $output_prefix.meryl > $output_prefix.hist -" -meryl histogram $output_prefix.meryl > $output_prefix.hist - -echo "Use $output_prefix.hist for genomescope etc." - -echo "Done!" diff --git a/scripts/build/split.sh b/scripts/build/split.sh deleted file mode 100755 index d43dbbe..0000000 --- a/scripts/build/split.sh +++ /dev/null @@ -1,56 +0,0 @@ -#!/bin/bash - -echo "Usage: ./split.sh [LINE_NUM]" -echo -e "\t: fastq.gz files we want to split if > 8GB." -echo -e "\t.LINE_NUM will be generated." -echo -e "\t\tUse it for next steps (trim or meryl_count)" -echo - -FOFN=$1 -if [ -z $FOFN ]; then - echo "No provided. Exit." - exit -1 -fi - -tid=$SLURM_ARRAY_TASK_ID -LINE_NUM=$2 - -if [ -z $tid ]; then - tid=$LINE_NUM -fi - -if [ -z $tid ]; then - echo "No SLURM_ARRAY_TASK_ID or LINE_NUM provided. Exit." - exit -1 -fi - -cpus=$SLURM_CPUS_PER_TASK - -fq=`sed -n ${tid}p $FOFN` -fq_prefix=`echo $fq | sed 's/.fastq.gz$//g' | sed 's/.fq.gz$//g' | sed 's/.fastq$//g' | sed 's/.fq$//g'` -fq_prefix=`basename $fq_prefix` - -mkdir -p split - -echo " -Splitting input file $fq" - -if [[ ${fq##*.} == "gz" ]]; then - echo "\ - zcat $fq | split -a 4 -d -l 300000000 --additional-suffix=.fq - split/$fq_prefix." - zcat $fq | split -a 4 -d -l 300000000 --additional-suffix=".fq" - split/$fq_prefix. -else - echo "\ - split -a 4 -d -l 300000000 --additional-suffix=.fq $fq split/$fq_prefix." - split -a 4 -d -l 300000000 --additional-suffix=".fq" $fq split/$fq_prefix. -fi - -echo " -pigz --processes $cpus split/$fq_prefix.*.fq" -pigz --processes $cpus split/$fq_prefix.*.fq - -ls split/$fq_prefix.[0-9][0-9][0-9][0-9].fq.gz > $FOFN.$tid -LEN=`wc -l $FOFN.$tid | awk '{print $1}'` - -echo "Splitting done: $LEN files generated." - diff --git a/scripts/build/split_10x.sh b/scripts/build/split_10x.sh deleted file mode 100755 index 2992720..0000000 --- a/scripts/build/split_10x.sh +++ /dev/null @@ -1,56 +0,0 @@ -#!/bin/bash - -echo "Usage: ./split.sh [LINE_NUM]" -echo -e "\t: fastq.gz files we want to split if > 8GB." -echo -e "\t.LINE_NUM will be generated." -echo -e "\t\tUse it for next steps (trim or meryl_count)" -echo - -FOFN=$1 -if [ -z $FOFN ]; then - echo "No provided. Exit." - exit -1 -fi - -tid=$SLURM_ARRAY_TASK_ID -LINE_NUM=$2 - -if [ -z $tid ]; then - tid=$LINE_NUM -fi - -if [ -z $tid ]; then - echo "No SLURM_ARRAY_TASK_ID or LINE_NUM provided. Exit." - exit -1 -fi - -cpus=$SLURM_CPUS_PER_TASK - -fq=`sed -n ${tid}p $FOFN` -fq_prefix=`echo $fq | sed 's/.fastq.gz$//g' | sed 's/.fq.gz$//g' | sed 's/.fastq$//g' | sed 's/.fq$//g'` -fq_prefix=`basename $fq_prefix` - -mkdir -p split - -echo " -Strip off the first 23 basepairs (6 illumina library + 1 padding + 16 barcodes) and split the input file" - -if [[ ${fq##*.} == "gz" ]]; then - echo "\ - zcat $fq | awk '{if (NR%2==1) {print \$1} else {print substr(\$1,24)}}' | split -a 4 -d -l 300000000 --additional-suffix=.fq - split/$fq_prefix." - zcat $fq | awk '{if (NR%2==1) {print $1} else {print substr($1,24)}}' | split -a 4 -d -l 300000000 --additional-suffix=".fq" - split/$fq_prefix. -else - echo "\ - cat $fq | awk '{if (NR%2==1) {print \$1} else {print substr(\$1,24)}}' | split -a 4 -d -l 300000000 --additional-suffix=.fq - split/$fq_prefix." - cat $fq | awk '{if (NR%2==1) {print $1} else {print substr($1,24)}}' | split -a 4 -d -l 300000000 --additional-suffix=".fq" - split/$fq_prefix. -fi - -echo " -pigz --processes $cpus split/$fq_prefix.*.fq" -pigz --processes $cpus split/$fq_prefix.*.fq - -ls split/$fq_prefix.[0-9][0-9][0-9][0-9].fq.gz > $FOFN.$tid -LEN=`wc -l $FOFN.$tid | awk '{print $1}'` - -echo "Splitting done: $LEN files generated." - diff --git a/scripts/cn_hap.sh b/scripts/cn_hap.sh deleted file mode 100755 index 8c79ecb..0000000 --- a/scripts/cn_hap.sh +++ /dev/null @@ -1,62 +0,0 @@ -#! /bin/bash - -if [[ "$#" -lt 5 ]]; then - echo "Usage: ./cn_hap.sh [asm2.fasta] " - exit -1 -fi - -readdb=$1 -hap1=$2 -hap2=$3 -k=$4 -asm1=$5 -asm2=$6 -out=$7 - -if [ -e $out ]; then - echo "$out already exists. Provide a different name. (Are we missing the ?)" - exit -1 -fi - -if [ -z $out ]; then - out=$6 - asm2="" -fi - -mkdir -p logs - -echo " -Get spectra-cn plots and QV stats" -name=$out.spectra-cn -log=logs/$name.log -nohup sh "$tools/meryl/scripts/eval"/spectra-cn.sh $readdb $k $asm1 $asm2 $out > $log - -echo " -Get blob plots" -name=$out.blob -log=logs/$name.log -nohup sh $tools/meryl/scripts/trio/hap_blob.sh $hap1 $hap2 $asm1 $asm2 $out > $log - -echo " -Get haplotype specfic spectra-cn plots" -name=$out.spectra-hap -log=logs/$name.log -nohup sh $tools/meryl/scripts/trio/spectra-hap.sh $readdb $hap1 $hap2 $k $asm1 $asm2 $out > $log - -echo " -Get phase blocks" -name=$out.phase-block -log=logs/$name.log - -echo " -For $asm1" > $log -nohup sh $tools/meryl/scripts/trio/phase_block.sh $asm1 $hap1 $hap2 $out.${asm1/.fasta/} >> $log - -if [ -z $asm2 ] ; then - exit 0 -fi - -echo " -For $asm2" >> $log -nohup sh $tools/meryl/scripts/trio/phase_block.sh $asm2 $hap1 $hap2 $out.${asm2/.fasta/} >> $log - diff --git a/scripts/eval/bedCalcN50.jar b/scripts/eval/bedCalcN50.jar deleted file mode 100755 index 440510f33bebc7f15028b904816ed79cc8ee225a..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 35204 zcmagFb8v4>yCxjlwrx8*$&PK?wr$%^c5K@_w(aB>+t&9y-_)FX=hW1huC;1)RadR< zf9|@__oW~W3Wf#*1qB84)G;9k^dAEa2nIKcEA3475?iG3;*p9`2Vb08M+yIFd8};+L<#NnVK*P z8`>Dlu`@9k+ZZ}KhpJt;;|`<$#Flutm!LyC0oFo7*yD-F8%=r1mg3vFFo2c05!3)| zA*)(AlNZ%)%iF2MfxUl1)5N6-`I351K?go^9v`U~faL^UK!0F=fc9rB@F&wiJ`CNC zXL_9V`@GG+f1b+pfw%(X(Oj4=bvV?g%9={ed5$hY7){G$x;!rT{sf(tPM@t`ZYk$! z);Pr%IBz*s8w3$2NzJl~lx?BJ6NAW(Vz{u_OluiF^kvY$UiWP^fo^81m|RRf$p`TufM!3 zKFK_}Pa_y_BEn86XU%?;={o9mtpAn^vwT92U z6~zzZ#%i-URd~RG1;;GF)bm%s6kHmPSE+`wgiz6ucQEZQYqYuQmgJ}vV6BJeS}cm= zv(gqE&b=*}e^DyhkdOuDnRTnXm-7m&n=a>K6o^}$oj)J-3IaoRlqtdxTo`w2qJv#3 zR}QE%Tu77o3MXF|jq)(-U~vxzRZ*b$K+h*W!6;2Cn{Aqb?!ic~H+m@M9BYs(xYy*J zwRTi=%@>+!Bk!*-Os_oPq9WJHHoKU9_0ClW097#V&hJ!8t3uo%at)2ieCVLk{{`2^ zVct81-!_1oLuom82Sh+NmFg*YN@~DC=!!3G+~FcOC_F~{+CP$#{{X5t&m3%Q!u0~x z!ixi`VfV}>%Km4~sfnd@&qG6Znpi-o$z~ugiaQJpeKw%Mpli74n~K;FgHpsYw_B+8 zShQIF?3Gw9noHxPBL=VJ7LGvqmTIu}Il{ds3j5RpPdO;dV2~sH0=*9_FsrwMdjn-M25c-f&iBZWa2@K~i?fLnyJBBJ~S zSNH?mU6ROa_uoD^SGnOP6`N$*BPQFhv4QeZ0HsP#2nF*O_TY8*MM?)@h)hM{1(+H) zz7wnlf7yX2G*KHy`78Ng5elq08AkKme#PAW(r@;fMBM#eil@m6#XqPmPZdHO7S3JQ&nozPs!3&{3^6dE;*3E)Yv=_%wI59ADBdeDf?7h{q!I z&>^BW*9lm$$+9Ia6&$#5p*Vr0fC!xH9gn-5s#@J&XKwcgHL^BAGcsr(=aZ~v>?zA+ zP=tqWs>Au%VRG+gnejNkS5Eag_=6c5ZV_;iCLQOqAZWH|9~OYFWXoS^2aUn2lv+BF z+iS0Iqif9*N`|v+A<#euHE)2h6T4#)F@t-61wc5ki?&K(AXK_#*XB=4({RK?Wx)-@ z*r`&B;jjttaLx$%mITn?EMd z7aJj?7|6+_lqU?_D`f{nVT}P%X6?lX{EBDzrs{M8bB>DL)P<$$6wO00cGI*@z^t8g zE5Y4Hh4JmANct&ek4OcI1*&P^G2hzHx`AGq6!;PIYt)g@LzYXpaF6nH`TG<$wou}3 zsIIh$zxtwV%j$@@ak?zY{GgpbE>5$?$-FPtT!C}bnKC7xYo-y+lH8XZR~x=hZ)p}7 z<3v8%*eG)ObNT1*BQ{WP6r)M+b~NHUqipaY<>wtqG3$J{aJ7dGNk-L~*^^#5OY4lD zUGDd3sUAnpGxUc^UHLS{E8^}qN@cTyb}iT6649cU_J0tAUHO*dPOz!JviQEhPnO;2 z1Tg3#-FO93VTx^l}F?893MBYI=GPGD`b$dWvf;YuKC||62NrKvFYW1_i ze8{{1#Hssb2Dquy?41yU@S#ET-RON` z1<${W_W@5lyuvHdEa#}2KgtSbX_KQqtz6dBk z`<2~XTGl*d2^y=Bb_dcp!dBrP8lDcQob1k*V{UE08~!I7j;l5qvbi!0Z;yV5EBN6q5KYqsfO;Av>u9{|__c&aPUvkcQ3l50^?0u%bhWv|e!U zKBC-q=a+OVpFjL3OUL+fs`o z;Ue@|iA$%DN>nbJNZq*ec?dn2LfceI@Oht&ap`E&;ExlAOY5op8J!51 zLSN8|Ep+BvYp}i(tconytV*Y@W;OI7xaZ$jLQaM64c`J(_5Pd;iEB9m9;c8twONKG^b3x`IWAu%dcve z!kT{*6K3qLu7?^8*8Fa}1EBd-x3%LX)EIWr9tIG9W=Oz9v z01|@`}q8Umch*RQcH!sGoQ)83*S&ExJJ4qMHxcvw0o3EZL>7_A}CON z8Ze%KRU^o8bRax~?vkHjKYbP5u|>?a_kvQCVL8$96t;d1OynPKm{Q+BE{ek z)8QDUSfC^TB@Ry{#F((b{sXQ$l_Xl|>6NOm*Am%U6(mx{g%T|-%Ntf#s~f}DHQT+d ztu43_K0n^O*_@2yBGrA*AJ9{}v)nIqKhD>_z4za&A3!-oh;RGIc)mi1fU);V3BHA< zRx8=&qDn^Wi^?p*tJ^`O7_}x73F>W=V#XCbwt4>xbmdEQ}D{E{tL0o;wYv!FAvHUMltFI5u-d*$Ow? z=CU?kE!AD@o2zTt`U*`GMejS4-wJz$jYLVOHgVMlkL!cttYcoqr0maRMb<@*EGl%d z7V?nT?_OTXir)P^g)E2z>lMVKja+F~x*gB-%YY;1#9sgLCk}X=tFqY6UwqOf{|9!oQ~DYJkvvDRj!?p2MQ!( z)=^I)+V7pKA_vkW)}5y*Z?ZjCiM6JQ9x_i?LZs{1MZR&CIMjj_XX7QoX1JJYHIZaS zSE^;BorTdwiE-sg0unpZ6j}}#2jQbqA}3X^cR7TZGJF6(2m`wLwI;4kj+ISPENz~x z^Vmf-DgtR4I`HY7Hjb?7n{!za)N!?=ZBkn+mFyYj20YqoS$KZ@vPG5W(uJ68mUdU9 zs#|7(8QIJceoFvfwOYNQaEhBKNP2of&Eu~a(6W^$r7N?KKd6WIZC}4PHaO(m%*hg! z#@l6rlqlCD_J0d7V1*Rc;G&FMv_&$wB&)I!xEs)}RDE`Y;~mUPyx0hMDfeEJ@bxjS zh*p6)AotrJS;9Z0=Q-L*+S*DY@k>ytmsS_R-BU51_777q^m~n%nr!nF2K-@JTrd>o zSX>9zky=~sW>8P5MM6vEi0js)?_L;EZO_7xc-Be|Dv}Df-`}uys3FW&;$qUkrKnY= ze5g|v*5#+!tL1bAwJSzA;|ij_4%W!DniIM#AhH{v_uG=`X9WsHs4n>&9jbXP{*zEl zPioC*P%AA4Wx;QaNsS1dZYIOIQK!D|hU&hjjP_=9^&YW2+U$Ij?=@qZBm~!|Fbd8EeJ8IF4H+VG?i+>?~oJK{0nBk!?9GV!Rivw?v zZ7?wAJrRGxz$42$OIY;f15HcOy%BmJR!yO|MhbUAM9^R4=acssAHYiCX2Ds|!i)9` z;ug!-Rtbr>VKSZ)Vl=@G8&w`a{f}6((574G4!Wp=G6-mHM6cne=#4O6!d`J5_6V!U z1L)BM`?ny*YQ(|18BFxAP@&JZTjeq{QQ+=Pf zP2lx0Ev^+wKe(bk6WbKhlB(|?r+`14L~y%6uw0PO5@#=$;ion^w2Ue&X~gT9RNdYwC+ubt&02 zkD8zCxxeZQ_|0xH`>YHR-iniJy7)R#O>qOMM7qm-V*495Q+i(wv^2YzLI_MS9glXA zb+eGheB%E+9uQUtN*19k4$R&%x)#dFq)|sAWO>LNN~!T}YpPtHbTXG;%E+f!< zl<^elOrERA;oFi{xM~`lpS(e0YT~Yxz|=G%>I=Y{`3`Xzs*G9CGzU?Uy@{Yu8w5F` zR^7x5X%}HWiF$7o;fT0NjiofO#AQW`ueyqEP+wS*USlOHQ7~nb8Y;&rzNRFDG&L3C zX0#QLr5Iu7)1Z0YBMYjc+bjJhN>o%42rTBt(YM|V)>Na&WYIjDy9;)vLLK_UCUuic z9U_%kiE2Ut2|vbme)ak%MY=?vW{ke4oIB}~&9qZxUa^aGiobV3EH}U%#vaaH;;Ei5 zQwg{t*J46)Wkw@D>A*zU45kibM?}RuOlAR&&?uF%Nj9_SqGPAmM9^aw`-o;*jyrT* z`E;_da<;bo{vfZ{u7>e=7Vl_unOb~C+C{vxv+&M^W4Z*yQ?ni;L5ky};%_Dg%O_$N z*z8BWwfMbHpK3v}tze|d@X1#E?`7BF8MiIt5FW3M306y0{Ary0>JkF>^S0r!oMBeY z?M^d`tBUBbBK;3Q!H2RmYjJjr&?CwO6_~w=ERTMy(kemf9PX2@G1Xot$O*WnJo2p4g*LYow$qC5&$s6&oich-I4H-Sf`UuspBW~`* zZ;q(4nG!vu9xvfSQ?;Xx@!jV!+0SX^Vg%>IW|$JLdkeSK(MW6HQ%(YyLq zC9&Y|D+l;U8{u2l@sL`)PZW9W7WZKk4*=sNipG7hP1{R)1OF){MGG_ClAJ&B+Cp4X zuyxK`y70pDc*x}`n=J552=K~Zx<&g^LD~3p?qo2hT7haWvuma0{Q=Lwb$bfFJZXQY_~@wprh0zR`I^6scYJzSXZ))ed>XFv zMf4Ll%(eY1%Q?u{r7OWdTH*3>^x+09h3hNw)B7^vn@ruG<%fXN!VhIrnJM)rEvSu# zp7hpXeKAArcsusKv8P<;o2}B#7$SRlGPZ{+@hhv$@!oF9$pXRq!R;W$rHNs?TAma4 zd+;!C4w5VQXQs?6kOhm+iv!UGEuK4dwEDRJxUnTb+**`-MhJXElTbhAtQ|+S9h%ic zQGj5x-3gdFV`Y9mg7s4Sa%fQLDqKcXxUqCBx^%g{a!(f+i&yh9u z$ihpY@>#U&*vVp4Xr2`4Pf&_)6= zxiKm~jkB>~SD=@0pON#c_9=Q8@ddb8p`dniJVk?I5u|3klL~ zUvxxZU|b%y)i`)8PVnV4H0oUf)2fHGuG$TH4}#uWVWYeV+K8_(;1B;ke0d}d39RZQ z{)7Va^@3}Dt*IsUc>~F}sa#KPl)kqq2Vtk~oC9BON1>EG_~2%@c;5{;x#f?fL*Axx z*saH{^mtf$nqrLhs#(9%zVOGu_{Y`}q7H;{+}LHMmr_Rlvh7db-kyL0ETNhVU6XMJr`>Iwo9=Ub#Vs8ITyeQkX-P9;4~)uRk6%7tCU={ za$-vcy{GsY0Hd~kviM58Lvpe#5=J;53O5uXEiezEZk88ebK{M5H9GT>a55@ihOd;b zsh2bWI>r-F^bH;gDj;2goAhGQ6-sR0 zl@nDjkdhn}ak5D>2;9}BT(^`Axlyh7ZKm0m9fSfC)^w)8OILBm623u~NP9_fRHGtG zpeR!RbgNOKCtTGQE(YGZF!rL*^3u509X1!QbL~Mz^ z=>iNLRUn-+4UaUFLSY+SonO5=2O!A&_eRf1kz&Z<{e=|3!U(Ki@JDSfP(XB6t5 zOz}h&PZfrsR*x+wp%<6g>|tjdm%THh9j;2~W`i7;bu^$IE=a&Dmq7OgI!hTjRUUf-QNqEwPDKc_QXJDs+A+r}8u@S1Z zLdHBSKbA7DW?^43Dk#J234IdBk+6p686cou;#LCHq#<8^ALyNUlo7zO;l%GfN^v^6 z58qINg}7N_8~oKmbfsa=;v{5qhb(@NrDoXu*Xmd_?}avFe9{fHSBmTdXJ&BHrB;1T zdQ~%CJyL0xh`={F?g?l6(zSCAcf%565o{?USgs3iZdbHTr>N?WF54TrUE(vHw+(Du>qxV9xtCa7S;y6r&HBYthF;xh3M`mm2rNQjoPWwRRROrEy;Dh} zrmv=_=uv=-f}^>nanAA*kGhpgW!~OGF?%T{yegpKeQviP+9&#& zo`KTts|`>2O)0hjFN_=~x;ZY3u1E~{49?2*2IHp5&K6kc1IeQ6awKvjBJb;%# zWU?OLx%x?Mn0mBA^V!FZvD^z>sd9tci()eC?4?!?7QpSurM)idNGNClHF(9v$trxY zs@Wz!u|-n8Kzg`K%Bta8-e_FTl`V6YZX~MtVpZP9Eaaf>(x)eHO+z>0fBai0V-hKW|XXNoozOYQNK zu=Vr(!RYrRjw9VhZu&caWpS2qA=~w%F-2O@K z5E2w*m}y)9uXFt`oHXvnZP=jwKQP3Q$idg+cv6|EJDUmwGuU7m-G7_wzc_#|SEqRT zh}*y7{=9jc0}%0QRRtyYu1x^;*9Irs>pL_mL(uU!WVp z@8n7R(@Z86>F}h;#7m>8olR46rV}&p=+h-MuG5SU9jQ}JiY1<<8@eU_epcYohMOjt z$H+XVOgQFgWdgW`rTERII?N+gJW!O*I)0aW+a=S?^0N#!Nyoi(lP%}19)Zoqw_hjZ z98M=??ntssdM7P;Wt+6V7_*Fe7fe0q(Uf`Iq$&5rOXRxuO}+FIKkej?S?=UZ2OCN> z?WbcHH`B2VG>&DEnB$I48WpK&#j=`2HfW^b(3)n>6WfZ_29caZ$tYFIrl*+7&BEL1( z@;h?5o(ou|?1D|zBgIFBb?MD|PMTERhiJc0mo=tpXTBMPY%%vmbeRr>FoSr5PJh?X zs4k)bD?YUuPs@so<$N6?gj`B!p&omPV4_zDZ#HY?JB>r`qEV0YsUv4+9ybi07&t-^ zNOqzGEhE+7I3oAG!ab+2#l193Qr(Bay}Pf*{kcDS^LSlDegH++7^2^UGNVuC8~-VQ zbAA8+fBb4EN+0N3-uNX+E1t@#Aalx0HI*~HW-6!|Co9$J*i1h@Hx8VVKjiMk$Y|ioe21Ji$Gml1M^kM1{QV0{RIyU32&O42Uy^TugREO~NGGpA@$CC@UMIVky-kJ{#@3F!Nbfdv z8ckKht{~bl5w?FNIDT|sw|=nL`cn+LCC1+++1b{u+${37Jpr-8-S!@C5Z2S5EB@Fg z`U-^E_QWn6`XVs0t=Cn6JLRz~DwpFYM0}jJr(rg*f0Xo8_#l2C)5TSCp7;i@z>jIs zC+X~!Tv49ZOq}amE~X-6I~7#&r#wwse&v_!R#+_!OX!hQ`=*+I<4{aaW)RmQ!Su#a zc%s+b!&fKR1Nc1QF~=@6W}lbu!zOV#cqn_Myo!TIbwUqxc923&kY*mFyUL|C>r^&3 zTtH7_2BsqJbEkkRNfm)ld45-m+FH`grJ7}@LeS71p{tceLOfmEym&f)3LP7?L4RO@ zom0U`1CksP1S&LD*X{O{b;6`RjnffdmGVuBSF$Wc1nvA$tt z_w>|{^^{Sco?kThy>~iBmaXAnxzUJ1{&^{her-w>mD72i?Q{hTQKfkIO~F|Fc_HYz zuS%`inGg~W!6=>z@}u(Dwli>`*p@9Y4(DC4xRCLs@{iPSUYOR zJ}7Kp7_!o{KwLtjeaf^On|6!l_kV9X*>1UC6$YlahN^c*HtMP$10s$|6?)C{|fR;Q1!9&b|avmSw$@uDq zc1aH62maNSyl;@<D^vY<)Xil0_58o?- zq6jlT#`Ae?6VLK3k1lnx^Hp3AwZj@})u0jeWvTfq$L;2m1<>JWL=U1NQ%P z`YoblWd1St|7-l;np^(M>HqhZ?g-b~c~RoNsbaz(t6S9|S;O;kJ6sk}p8#)jE@6PW z_5oBuYZz~f-*}j80@XI|-*@fVM>a`pqr0t$rcrdAmv^1#7vt5M&$sR&n1}Rf?tRct zLf{)!mh#E^ji)c}h-|)iN{TN&Nh68KL14otaU9)UzUIw=#ZMdcdQnwejc~4)R=_#x znYEypdZ>B(ROiUe8;x0@ZXZR)m~?(&FWp(pAfI%yxH>=DE+aM{O;&=*Mz9DVg5oOC z{8cu8S3P#(2`!)*c#gIAD6yN$2i6cMuFKB0;pFRxHry< z8Ka7wOyVmuz8@V4M)*WOFhy6t-6b0bq>yQf#`*Vn_=<&@ocsYky}xRTZUbuG1aoov zh7Gu656tQF^9ErO;ONi!fZ%8 zVve7*7akhtm9fo8DCvRFLAgm#IBG3`f{Lv_?3=2Rtv^P@n1O$q90VlHGyamPvrnq0 zdsNtjBMr3(lA+hJcaH`RnW|ZBJ4V43Wok67cgs?Bk%}%NMN?KxHCn)c>S;0mv;N9E zhdIDcqNz8({7C-#e#!doeZe5ONqjl4mw&Y&H`dHkLtpOf53$gZy~69RXnMi`bG{{e z-j;XR5>1)>yByd=QOO6I6au4N-=PuT;cppvg;|PB@+=!^bhiu>TA;O2(NR!bk|c;0 z-6`5gzEe)V;(+tk&h01bH+ zVE$Mo+lA}ykJJ6<_9Sg%pM0l@@>~F+?(>)^Lw>@7P%IephVW}0m$c0zmc_6nM0L;j zB*AV5)WSjXV;k2S)A-^YyTVgqAsO!C#VV>a@kNE?ESOY~q{EH2a{Vp`qfZzGah?FxUle#p>@;PdJt8zzp}%fIPkYS zx$Gzchzbb61Pnpu^SP4AVRdo%xOOf~+7T7Q*z!=0rlUxXuf(B-YMG4->}pkiA*w@y zIex9JHrT&bhc_f}D-LY;^un4VF94S9ct6=#2I~?93G?Cn0n3Uc6h_EEOF_bL?j$%} z8BVNrcGUOs{mQG1@l|J1b~@@;G>8ci-AwWaE?hfb(9tkoEe-M{SJR}tq#oo5vvVG|jcXs8sPk#v6Jw!O^=>CYRY9KDxtyS`aYeF}Hbd^ih-J<%+cp6l>R z>x6wGMX4rgiPwbB(*?ml;bJ^EvF$U`f$w?CN$OG7tgf+uyER~51;LC;W*}aA*KwQ= zp1-%qRSe5pI#aXh$&Z(vjb0Lclju)I zE1@chEy@CIen>tlA%scOLo=3eL6{gJ1Ityk<<+g#Xox70*o(9JiT;kio5BPa9Y{C5 z-xq*^wwzLS!dM5=$hkuXu!YzYoiw&XLD4G))jx`$-(d7BCLv@S+ww#sbpwafRS`JO^e-&=^w8wF4Q=u z$&4G7I$bgqiybT7ZEwsB{V-M%Q;RL1%y^17;Y_Fmgaf87m{ny9(X;%lBNkm!UM}kh zE7&@V2b7ZFG0mC!kT+gv?dj%M=^Q8N!XS@sEYZ7`InSrJmYJ5#%L--eALp(?dP21y9x|aeyoeg->Enj-5d*6EYEG3l3bs zCV620Nb78|du%H^DMX6Vjy!3;A2_?w*JQskLQotBz?@`jN=miX-2Yh7N@>GpbL%qY zWL?w(nMqUjSP64_NPCje{TjKVg%ekvtJ-CeEB|}{YF5X}$JbaE1pqLoz;uFmYHQ-Ei-Y)JOVXhQU5O|Dm#RuoP zamFg*1^?6{RScc-^7#SUYElKXN2D z2{81wWYG4(ITw~qR;{SxeRj`7e>NpxRdiZH5XUDB#DLtuQAh>-3d6!t;~BwRAc1p1 zOeM8!e&3Gu4#sO~aTELIOLs&Alyx=1wP&>~Gup|b+L;cSJF-U~-?_%yu`?!Jh%MR~ zKk~5Ho(DHSxdYX1g{^I=gom#n6Qb^1u!r2h?y$p4pO!FUQV)<_)YRL3daNT-mj#~d zq7{V_j8Wz+Y&RVy_yvD1pJ6zsOc!8R2BU)85_2V5LIPqm{C-nKmz4{gLB(IIc-^LaEYWG3uI&suwWvLb{7 zzW~AP9ZrAf?iVdiLQzOr#@e;+n13V9>62h3wPC#JF2~wIXM63HJ%5VE-yP;0G;w~_ z0&ZLN`~q1SUqC&KsegX_#ys;WLpob^90k_^G9`{zc00A4HC)jRW;;!#c2DrJ`4iaQ zkvGF+(v_~9dZ!Xgr#|=cti#0EuN`Ep?^=f3u9P|50%QDBX#8R)liVj|Y;i1CBP1s_ zvUacVFR#3R*#GYD)^DKLOQ-ra))k%MUoI=AP{cY3PL`LP5`KDaR3iWKH{1QDbA|Hl zjJW1by!tm(5lPXa<@JvN0`svGsa1d z+1YleM!9pG3KmzL{orFy4n~{p+(3@v^5vbgGOnhuGM$_lOqt1=#h6{nc2#c2!nm01 z>_FM)gut}kH4^~bUoY`w+r|w%zo5wP(RjBZvDA}eOgcAJ)``!JEsLE7GPLAcW17XH z%`{p~p}ua`3Q&{0aFaA9dDD6>A6FNH{9_^4v|DW!3M+LV9&ylqlq&8? zL@hT?I8l*wZ_-=di=z((Pab<))(}7BS;q5MaZws}^5kCk@8`W|&wivJGmEDiVwzhLHT;4;pI3<0c<)c5FMp)HNnBgx2X(J%`YlQ%ce!ZQBnsRuPWg=&d?0z%H zm;CHQ?hVs~};pD4`4Aa=-M^n-ywXs$K5cj*pW>L7E(s zJG>$5=hj8*TWiyhpM$4ni9BAZLZ&U;%V&O0hMQpKe2AY<$bZ%4t58BLDgqFYjOc%= z%m3a-M)r@+{Lkowq&%aMp|hniqnM?QsjQ*3snh@JCW}@7aK>Fm|6WOJZtQM|1(;ug zWv1<^lfXa+a>)|O$-w<1J^%)w^o`A}oQdSbh&G|DSAa?>fPpuqmAD`Oney1@{L}6#}H(thn`gAI74+KLDYxnggdiKwmr;!WjGy9#w!A+gPm_ zyem0KD`NC81?aQ|acc!CpRgJ@)^dZW@ZbjiM_G-?))W}hfGm%eI{8E^(=E73i zS()tM1DsZj(U#=(Z9o1LAdaH`cx^TdL(MkaQ$b20(M3yE?A7`jwC@WM1i>S_?7%B1 z)~}9~nk9uzH5M28UXwkJR8Z%ciQ`7}wLo5%OC8_8Z#FYIEcLMNEaGNiE~_0+4jHvh zE*4^ConLD;{j^{qLoWoN`WSfcA(_cEx5k8DVTI5*qbI^{Ls!6C!36jp<7UwL5C~9_ zURolN#e_XLS_s9zA1cn~s#9~}j+6?P;FvWXo%Jn>tp(*A#smqkjnA~kjmmw;A0J_s zfJm_B=JnB9By0_k!ZiayR8mFtfe--^aRun!sTp~a9Lh1aM2Cc?>CI1QW$UW$lh?)& z1OpLaQCP!+O5_vH@uvhFw@?Jj3YpTPN=;^^lvyP!!tN?Jf7Rjn!l=MyH02I`uLgOH zaL%ROOAxLRFgQc#5sBiUEu$ql6t#=j1Y3r%;-a6y5H;s`V(r=pSPw!O1JDRo6ZKMK zf}NRJbDNYTL9x$P&YP!=*xUe6i-ID{$h4LSGV-yRnWXa#AqFwUM)N}sY{uhVq>qOb z^c==i_iqnT4<$jildV#w$GMuxDotYYSxxt)1y)bk+eR9~h~Be0HO`2_+Yqi{-xeGF ztF5)r95m3e2PtHant3j=o#LhgmPN(&A)GBQ7aoeEoNR3AGQzui^Es(Z4mga7lQHA5 z$1$FN65?wvaB#O#PS8Q?qO!BlIeTW-5EAa}CxG)3m)lHL#}yOVogqz*KPYlU$uWDb z%j+^yqmXM$3U}BXawV@ zKNv8JFVIv5?4~-{px1~C6QxO*UKt@TdqenDg3vJUn@>^4A#iWZ!qXLMnXNFwch?;o zZpmbDSBP~bfe2Te$;g<7nI2*6y9EZTbyA(kYR9oyBGyhEhAGX$JZ%ZWQpE2<{P~LH z(3284(-;Vgo1rZ~JM9eGQ(cO4@d?NxUqZB*EIP}8sx@jtVK^8e-@@C<635jK_BotC z(%A7La87sZ%7C~jL>Jv=gSu|1xo^+j*b{o{=8XTCr8@}R&J%Uv^t{K?8!7ivXJj7- zS6TVA#@5#qiKfrg53omuAJ$F<g6jz+RiA`ByrhO{^4a!k~%vl&I?1y@NVrx?X9{&L2jic?5l%?*u-#h-<8=HMg z5Pm7T6=QntH!6Kk>H;u!*=Ov&a(e7mH~q(lFosa5(md8vuHvqpz&cpX7)2d=VDs`2 z<&&f%e0V_HoxEe$+JX2|hIl@oa|v0}ubBakr58Ku6IBm)Qm33v5oi908N3u2>-y9F zoH?~Iu1(n1#h=Z%m$H_;|7H{>`us!v3^Tj~GOYL1m}O*d8D+a@6qU`Fo$lZidKln~ z5PsYTu(!784PM}te9$fSn>Da8vvlg#R?EON$XynS$*MA7Ri&eFd}R)~mm z0K-##6>e>b?vpB5S&~vnr@o1edkMBG7mOOe>kYB5as5V;6;FBa+FWz(?RTu*kGXJS zv>lo#J9n%Z|;b?C^MR zWBM*Bd$zO8BXK2ZOwCTt<@hd)FZd)^l)&B*>B*fN9mdYujn*aEYMz>aA`+L%Rzxc{ z7pVn2Olv#14bXMIM?X_erJOqWQ$gVN{zF%ZuYyh4CQCji!lTQ^ju%X8alWj|G^v+a zpt_({gK0(8b%fUvJ!!%y+cj*4s#FW!aVsy#sa4#Vx9Q>g&Sjay?Q@&#pbhyje{ zp78Q;4e(*QEZX-@_%BiLVMc;og#-e6NBmEs&}^0``(GaPKNHgbH&OVPM_j3G+b;^j z`6DcRLnepQUrCM7RzgI?N+plHGa>y393^#L)mWnF*2>Z@{owDa_pU@0z#~CbSfOkR ze4IS2yK_fl+JrfbuGpZm^eA7>_L5#H1Jv}MUD?j1Jmp$y;V z`3<3NG4hAB?AWuWK;B>@q6vpTszKw+ z?GK6OU!eBvHk1}!DB9KR+bx596bc0{eNu4ZhU6AM0yWNzq%RHuPLqoXd}pa|7rW1u^)y^Ieb~VfNpl5 zjK88E{Ruw3KUzKY4S|Uo38&l^6&<-jbhK6J7ZJL>dERru&$PG7MrUg_l|xCbfU2+8 z$UcJ^OCF?ZRnT?@L6+B+6b>1oW#Al(LWA2j0q!Qs*TOI7-)Scjezm1U!!Y5`ydlU% zEyc4&^IPndNKugpo9m(R7Q4VAnrTGRwy%AEa`QJsd4GldC!5nvOS^kw(}S!(zLh{-Lkbk3O5C%7}n(o#IM?Z|H*4@%4>-L%75- z)DI7v-5n17HCP8gMr$*;?bd!eh+mMWFg8>5=h>Y!UL#O;73(wxf(#-OI zlk`v(U3FYlv~5m8I5^NTPVnEsqHu(=R`W)(#R_yB?1HulQM<|G;1I$T)-27h*Mw)N zb>~*X>D5PwxB9Qw(Kif^yV)tFL0Ouayv{QnAD2zD9S;}!0zV)I2;LCCK@Pm1Nk5oO zKioViQ+4IhFDoyK)wYs!o}#W>X?;PF16NfvDvN_wV7(Wh!3KJJ=$UdB309@~>nBC+ z{l^_mS}LE*z_r>sYS%xM{!qT45EH)Id4_+d#b~<)g9m;OdA=-_3Yx6bN}(FMXjlkS zCR3DhYzqFufNRD%MaFvP1c7)5tH#tbWuA`)4mD{tk{~)1;Fym!G{r`4Pa0!~EMqJz zb{eG%?sKi9}^(4up^J_A#8A)4{b--UQ+ld z+Z{;_8mV638WL(b1%(xVxe%yR_MJREv8@LNX~DScU0&kW(R~CfsMckH2Lw2wHcIb}FqSaiH66zE>)E174{iF=*n5 zZc39XOY_+q3cdv(8WT=u#gsHDgxuVVy+UwWAWUKPfft2E~#vrpFcsa}^lcy^@E5p}K74^Q)mgGW@cl$@)VY z{wPiLH7#^DrEq?<7x#3Hkj`o2``w=hf`V|! z^cS~Y{KYZ=Mp?X}>@1(~>LYY#o7`L%yK~HUUca%9SJ*shOx<@#1FDS8G~ehUvzd48 zNE*oPlT@eky2$4qvHXp@g47ae!N*A`*X zSw#jqaI~)pmp2O=(R0*+!Vn=snV#A)%v8}cV#uDgl?vAH)MD1{?jaI#$V%k2Ul$=f zu8mx_2V3j~&>B1wR!Hqt?jhIMTf(C}*PJMyBQ?minYRPTJQte3{9b)EhQvsQDa-t+ zoa0$$N37P$_mz*bYW5!}%lZflTw}F-bt&1xam_<0Dqs_g{A&w>FNkWD)49poXYP2J zW@%%yvrdDq4Uc6i;9wlsj37j3+JW67LF`{3~^Y2W{XNmOJpH1F~D|H&;WgCRVWU-L=MOGnzbeNM*yWT?7_^8jEoh7I5EUt;N@6i+unXGTB66D)5osn zomL%M%EpZf#7iP+JtfwG+)m?>>Vnek$DjMw{t|UFzb^TPKoh3stCPN?C3!cGOLf5W zH(CzG_NmD}3^}V}LP3#EW;k3pr-?(rqgZR$$aP9qHj1m?CG{b22`E>4NQ_NkR!oCf zWY}o68NlM~(S%L|*d3)1QYnqLGTu!QhD~~5%OM!%Nl@R?A;-1~bh6X0bq}40=*yf# z_~nt|cbHN;Fyvl~e5o~>R-cQOF24$eE2IO21i>3chng))_V00U!n>HdSY+H-DNYY( zK3)u~ycxYn4$cSqj+Vh~jrPiffWK5K)wig56e%8^H>^9{Vj2j+-&}qy#tqb;>46O!|B4aLHXj5bRQDrLES32 zjk7$qwJWwuMbN!Wwh)&{+zKtZ^qhEiUfBqvy)sj|dLbtTMP6c5?kB#10>=4I-3vBk z^Zd@o?;NylE8Gt#?c_u>y;lz38WEk91*}(hKeJK40K5&BSj^btMS$0V0{nik{MG*8 z|J`5_lTmt2IQ%d7kf;2yG!h>&*G!_724wtnViX`GLF1QNUKF0NV&7;cnzKA{Kl&-t zPhEO#{l|%2gnL+bz|Z|W#4&wfXDQ+5<0ta;)F)9b&L`>a6N`+Om$$cgAM9Lmx(1`U zw+9~uw@_Jc8+sa*H(~N5hs~FYDZuIahSLL=RK6K!LFQtN62v7sIY z30*lJ#_O`=y(B2*>)Z2W9a9C;4R!q#K+nwst%iBq>Dii7+bJz^2P* z2;^7!-XwaW!LS#y1UB6c3$CV3PZjhv_uD8Q%dZ$?)1&D!m^_L`7b#IQCHV<4+m(gA zs#{Jv|61B3M=OQIA%|`XRULHukw&(fhFrN%rI3r&a!qE-u)3koXYpU1)6x>efiX`; zS7{^=~e@t3AO9N`QQx9d^KDjjj zM<5?cl^rz;sqAyNxGwsVM-N78`!Kl8)l<~ZkTHaD9kBDSk)7l7n@^I1PDz;0g2*1_ zp$qeXw*Qo28WXW4$DR4^(jM4{)~D1N7ui5p;@uyDndNi18m{+I}V z)tVpw%|v+9ml^U><~_XdTq0*Y(CiiU7E5(kBz|Crd~Qf1$&%BSTO>!)>5MN>Io6S{u{iOsXuJiiDx2rrtwjvqD|G) zU@q>QWkyS5Ghq%2KdVsqm-CzJHo7a=Oi7mOOR=%Vnw3k`Dg13Bfk+TynQw16d9k zy3>b3!luW$;XcY%R>;`9Sx6)A2@&^lvOu)3G?#RF@|R#(WnS7ITnjcitP@)b^m11( z^S-A3sUz4k7{)bCk{?CrQv;GGnaDhHm~8SR-pm!iAPm7u+pZz^#3hsT6x<|~p$C!{ z4iED^#z@B6n=%qBG5N~!s*EIGt=DhM^&cxEVQV`}U5EcxMslBJFkh8X5{NSDqt9HF z-w$;mU(ETwBBKaeM8wz#P8~_ftRTlO`cbP|szRPJPYToNXP@Wdfi^1gi35#MwOGMw zERgI?^x%vrJM&Q)se!KMxalUJnV| zo66v&C~|w!*#zk081j3Go4Ms24SzJ-52$7?qZ6}3FhdM)+#NxTzrL?7>zlppUIqF_^V<0C3Xq)oHrP%=^}9b&1z@-fI;rzfFGS`@x^ z`t|GN``zjP@x=S@PJf$tl%kU3(vmVt?+Z#@UcJ_oq8gWM{JFEaLqR4rGAyYq`c`*J z*=Ce>M2bqT`1J>1f?`EmF+GTANE-5yF$t9+%5PHQ)D*D8&SFyIV^U*NqQhVm5+I@G ze#~~L1&Q&fJAB^{3w$GHW#zw5F^Vco1C#-&kEuz7;iR|{3#vx@nQ_7cygm#Vxaz{K zMbxX#cz(TpKipqOH1XeQmA3$lyrHha@6Fwz^6TdA85S5Zi2T+MIVYEzIET1V6h&b# zkUt_wq9p)OKAi{*Jyo*Dz|ZIYBd4$TeQ#?(^U*7D-Fo9o^YVj%9l*d&IdfCp?d|cy z^I;axGca9H)?+bJbBi@LrigAYOZke81AaQ)g%Ja%dErhhll=5eFX`lf_y$uM71531 zG(xe&0$>I`F&Uc?M|(OGmBbGO0PS&`8(Y~l;2nsd>RfMFV0Jwg7`UDwsZze@IoAF^zBq$T;8AySr!)<{7FL@E* zC}6;8@tt|e%68-Zq4N^Eh>$i#tP~38MI>J>xHiLcM&D<%F;XL%M?Jz;)x5nyd*4h3 zpxqAPVqMG?%uGGWUVSe~f{_g`DQZoad1^(v@%@k*k4Q3eeO(wD_kd4Ue7sc+MY@?0 z0o8qzh7|G%sXJjVDCSXHLfZgb%73Ua^H03Z6?hj9p&Aue$2|MJ*f6Ni{9yiU%B4cIOrzWRT73G z_ZvtM#01v@qM=_WMI4v993<8%(@TN+N+Mz@BJByK`%2Yj8189>ERzZlx+ zMN}Fa`b3=IWZ<|R9~MWQlK1g8{84i*u{>7AT+e?P8wX?(){NxbBY`$${O5|CVFe-1 zpk_FzqqMcgPhv|DOLmRwxDmNhMaOm#c1`H%kl7i#owT389#c+o)<0O+nTcG$rYXftj)2yJmd_ZFq%vr%?>ksw`so~<`QxpPTJ zkzt--DD*u?sJ?MlG0iWm86-Bq3UuZ)ZEe{va|oefz*M}1wFz;u_1;Y^DE3TD(A$PR zh~MGHIq(!dkgv=fmdo}@BGd`BlWmp}e3cvGQ*mGYu1d98fC#KdxU1?aqV)lzbef_d zcnc#p{zKDH*%I;~ykE0Fin_1)RYEs7tW6RHvpUjR4%QmLjt&l1UvmlK;bE4`&R}U2 z?AQK8EUV;Ii(226Uj?!K!j5c}s@Ji8qZDVYfL`l{Y~*lyE0Sj?iw%34L3e!`jvwcda9+m zX0|CT37PA?z2XmC?da3D-v?-7`b!MwIafGJY$RK-TSskedDUl^EfFTnQ)`YfX*znY z{imPLwcm6WdA;4Nr(3hhdQHTj}~lcT^M@S`&yNMxbYx`d0K@5GXLa>-Yuqudz{bh z*_taT$9m9VzPR4-^}Qss(--Q;WuV!yoo_2%Y+4(Acp{n^&aIg9#iV04%7srB;0tKq z3S(_l#N=)n$tXCKh*Ku;FW(&S$-^`0<_X4Zko@RqklU}TfSX z%F!JRk|god|05bEW7!p!EQct9R(&-paXa|Mz!&rom&Yq)M(`XJW*D1e{|?`SkRSZ_ zh_L`wDP^Q#7Ua;E?(Y;!Aom9dHrEn_!E`u)o0j7WCU4z(QtUh$-HzUeJ-ey&=sb*<@D&8`MJCSmG{s+gR-HYSL z%In(2;SYDqUrjvq>%@7R4*yf0@_#uN{EkLOhIWPqe>L?H%H|44#>ie;3)Qx2AT6f6 z8l|~n$a)Raki`{fOf2y5LVO5GXQ0Y<_HiH8kEJx+JDU<$A3|xpWJQpQ#;(kbABmq} zq^);7>mVTV85!!1PFyG0Uq3AtxG$GWryrcQ&`?mImuN8S2GS?0!x!Qw@UEoL#Td6g zCBp`{^pm=O1Lx_E5VIzAe1r7wpY0%Fvn8^ds9B!eJi{w9m2 z%77!TW@HgIP7zl^yp@}mjwD}vK(B9P#v_pDUeIBj+>RF}V+WH^JS7*QZnu8}T~@4b zLUB493JPkX`h`R{rwR}xTOWn#e70O5P$91fhJyvtI}F?>heT(<#hZhpyCQBtFyKqT zCO#p9g=(QoQ=qun4@m|WKp8Q!CZseB)wUoO3x%47fDO$Bj^#6Wq(}K5Aw)^^F_ipG z#)|zn9wU@t*|{0r2|SmBi2KMbKZTi`+ab9*tMrV?ZxRf`L>!NLtSD+Q2$y*=;8rl5 zbfP}i#+f)IDV;-`pKTK7)FRqqPpp^kXli}1j3QDY zU4l*ICqHKl!xlXbz!V}gL=Uk7A(Qokz!^oJvk!pnY&$%>z-_$GD$6w9a09}={ zts|QvFvmINH zE#lxt0#OjfMVY1*g!34le2yUy32rUB50xURrC$5Wu@0~W(-!%QEv4IH(;$e4${q0K-lEsu~&kn{u& zE{vnB;HK=l5Pt_8T!dh?DCDi>KdUK*}iP1 z`Ri$9Espgc=VqpDW3-@tRs+n7+LM1fbV{`s6vDe7)|b5~vsh%V?5T>@`QYy~$K47g z8XU4`up!v(l4N;=JP<`!zr3noa7wih^_cavXZLVZd!rBXfYlc> zv_R=j_~kN}5qGGJ(*GPfe8yETkPXO7&yBkC(}-3SCsNQ3WRnv5lriF`LPm&j79cPBX+^22b-y^~C8Cg745wIwxLjnN(6w3hR+XAol~$LuYY*Cbbfe< z$iLQvbcKGuwxd+!_+xEH_5We4-i%s-l9crl8$6dpr(QNvg27_0inUL_9xD{kN};Gc zRW++P=Awmxyo9{;bk9piFPvK97V^Tj-Zar#+F;6{9PNR74gVcQ<7NyL;cOUfl`C<1 zIhrdCYX=L>kzL8JQ-P)2?w_DS{zkU#p2VLtOTBEG#*$8o(JuI=oqg}39CO3%qU_o$)jG7#5_h=2DVc1VAgb`(+ zW*DNf0W@aF*m*@I)n%1c-V2jddbLq!0k)tPK<;OD#4}x?)~&Lpg+W}7Tfkk8lfWrR zHQEw8Kj7~BNq#r;D{7?{&{7*D-3lzb5-&pZk*Bd>L=F=!%E;x((Gk_eI`qjL`0=D; zPI$J)SwSGsWj#}%77`jJ?#E6H2Y)$$)>sWi1WKh9n21|3LYO(}w4yH|%6k5|zpPGc z42n@&DZ`#p)`xLJIb6-0bw(jB21R)Wsj52fNcmsF`>3Uh_~*0?E>)fc0lmc=d6>4s zh-y-GR#aDe4V961^b&2dl(8=L+{+CqYka}%c6Fy6uM_$vW{$~LG89Fbo0v1qvYFkG zk^&udPYFKO1MEmpEJ4(xq7DgMr`kvZafPZ{K+Egp$1xMgUdEuS?$2EIQ2RIC$fUvA zPTeB=D+S2_Px)4zAE}~F963OLi!@g^UCEJsF7-C|=Pu;>?ifo_E5m-l!MHBD%$N6! zccE@(;F3?l)YUrUWKk#*CMvaHxPt{ewsu9m$yNV@6MW_lLGjEd!s%cit0`x*3H7t4 zt$7|E6?PD(DTqY7N=#x*WS=d7SEV=`-R>_V`PGJ+%03|yE?p`V zJjHlXquho#g%C!%$|p9K{l^rQBHHIPi?q`VytkcK*CH{7;H&R5;E$cxUzd`98;p@R zG!}BP`EUPcw1i9#FLL0-Svd&sd1%PUxXMva!9YF}x@Z`>xJ|E-Iw@N|>I`ac56{j> zUm#Bxo?$yRiPMhW_?Yw0<*|;->!)*wb--_&^3d1h-WyO656fvLD_@N+c{zUYO8)q= z%gF$pK$hJy<+zp!NopW_`e{XXsdXlTMd14*lxhxb$1t(LaI!1kP=`(5N{eiDPc$xA zQNBP&LH1DB&%<;@_t5~^X%%T zXa*Bnn(0vtY7zibwalkpzI2Yqe+0;k-4#@XUfq)$uh(x6Ub%zFAJNp0f2T43yHEcA z-IL(rQqC1*vTyFmfH1TqqQnm+^<(XQ*n#V@xnu1Pfi`dx@bl{*M2uBHG(ZZHll3f# z%^x}Vu`{~FCBVrYZf>s~fyuw)=QAt}(k}qz%V$N4lMa3i`DD&ZVeZ=U59g%+_?~OT zD{$mL9__Egjo|MdO~p>v#^!%Klb(u}GDv*zUWqr&H*W0$a{TL%U+b0FR)KR+U@&v@ z!E5pZ?)KG_?6yzTZ!3WX?&L6q^m5l8H)+pR(2b+w(1^8a>*_u+PH=u=Y{BdN>H$Re zWsu?rdeChJ(BWi#z|Qw?Xjn$yFOE}NsLPGwpipdstoC&|roYOq>^At|?}H_RH(Dey zRx=0Pv_(C0f8ixPED<_Pe85>Oewxf@N?*wGs;IH z0I%=HTS{(0+am|>Rg)FtrxzP7;)p&J&4cEr`;o}#2aLS}Xd+NIbyIW~Tkk?y;x0ls zf-I{3Fk31I%T0LRzJ4!rq){YH);!blz-JU^bPn3&L`<#EBB<=3*Ijn%9{Nc+v~7u$ zpRdVHnu8%wZZSXwB)Z^$>@yy|D(#lRY^eK2iY3n`teo{#CwBz!X@-8@aNIv_M74Wy zGff^S55v&2`{wUBwgE5&^JUjXG0_}xhXT=n@1CT1aZrN?rnp!UIkgmim)$uP zTD|3U^YPVe)O0d0NJjk3C1`+L*iyINq)>D~tcaa3CH|P8tlsq3M15>lc4v_S7KfNq zPzn8OQB?OKTgtn2km)KkG$q)-9A&(+^ z?RWwwaZ_HmO={>6qYGY`xrM0^Oah%2)zBs zFF3j%K?>-_$?+EW0`5bG;N>ui2DeRJSm}Krr}4kD54EgD7CY&{z?Fo|>TvcFn9_%N z1Q1Unhh^~kG9`Ec^}}byNEvB9!@TKNL7nZ&oY$?o$iIy7;s3h<|2CN&<)q93>EXF3 z@9;o#sx_D@xm=1}sE_#|EpvA13)@WNERqFls20p>Tkm2FyYF z88Kl86qqUaU5Y+SnE_CeCF1Cm-1xcI1bR|8c_W*SP|0laXBQo76BB9iIKAIJ*t9I5P@iMVXf4qEvob`B!E0i!?p#`KF3Yo(J*>iP6XBe<&I40~`zDR|UiS`rd!*Pg@^@KS~+@Rmm6| zT3K5f(*C2y!Ec5^hANbf!o2Hqq(cGa@caNcRJ3Im5vRrCz<|#eJUTr|Ot2*)`kfS( zp}BtLXkw9A*Q*#37&VDBmV!g5Aw@UAwPQ5cDy8_=SkMj31kQ!J=|?$F}vKTt;3yyyw1m2viATX;UHI2jjW<9 ziVBen*zdW2*5#Xd7BVd)BAcx@DY<3)Pj^dTr`t%ij0C_(`@YK^c= zoW?7IVH&p2G~w(JXbm_h%I_eO!KE@GuEmKFruw;Hn9@XPmt2{St1Kd1$WM7S9x?PG z{e;afbE9p9syore;HJ@PfDrYbAaH^}FZE3Iv~@t@{%c*KLiTgm^ND1*JShQjhdWX^ z9!^?OL7~vL09Ye>GwLZxVu*UJr6jAP+f-l$M`e~%^{FU@B`nsSFkA^@9y^{O5)3QR ze!*9VX=pgc!A5$3OW48|_62&(2%)ijFjmOs)_LYe}mKE!;1zyOj;c&^8N$RvZ{BW@H`=8}cbb&<9yh|=Ue$bg+gF8>|B+r>6?_oM3!Q2)r2R3Y#Xhd=?B}8cv z05-*Tfi5(u&y1<#h^WzxA_NbfK(|}S6+N1w$Ki!)c`H-MfpkKe`vDAb!t|d{q?D8a z&Nkra3b)`e6xlPmHm;(5d3x~5{9Hqf%RVV3AUsjn&N*R=rIELcx!3o0n0EA!c9cw| zik6%^R6&~``?vg)on+ zR9NhP`sQ^L3w&n}^b+Kq^iAN0WKau5TZ_jJ8uyM0aI|v?99l+0sA-1Va`Q&2E^;QPykf!YuC*JQmTw)%YXBu%i2)tFM^il9WMQXadaV_~tRc(b3JWMv&ywDhUibcXrW zRY_~sWkesS3?XbvDGS3Hy@>2)f6aA~e!_0%-(CO3Knkm6A-lwj@B$Dd)MZScr4_Fe z>3IjPCSHiLBa>af)S?UB-;gj4arY3(A!zt3Wmv_0u;o~hYBvtVWubqo7A|Gf7~YFG zcSKJh6-0d2gK9eyc~gXa9(y^CwQ`I4&Os59R+p+ruaZ_Sx(qozm_5!zIglbm7f^y# zaTQFvK;NWSZg44~NLCSJL}$me8QWMO0b`aca2gtO;v%3Bo<>;{HO?7D>9cr%%9Q%( zuSUPR1J3-iYEbh8W))qW^QqZ>YJZM+Q0QJ@n=@_EoBBftD&COQmuxmxg>oCv_7P2k6MV-=7Onm;K1LlR^Cqe zwSlp&5>0OdGx!a7ZP;gl|Lk*d-!<~QzZ**^Of$Vq=y#L41l0_s+*C7Z6@K||vOOB& zBQOieN}JC4h`Bw&49&IBX^8Vl_08hy4xyjFFu4y_JLcd*x}ZEqLUVwrUj@2-LwhB8 zKAU|xq-;LnbaQ^NW(!W+-tb|K?RIBa$mn(-ly#W>gmJdYnxeMnMpSlBQ824TviJ@7 zMhmbrQ2Z#5LAP{SrWj};0^HepgRQujH8xeY#XLnPnh9Z%<|2GJB`$q6_$hr@aAouf zJb|hghQhnu6*QJJ$O4EufMwXSDbkG8yZH%}1E&h92u*I;tmv1(_Q|@vz>s$rPdcAo zSsm9l-naEh2$v_k1&*M8{FmUPdALhQ=^^hTH^+z5ADF3KD}zgv<3>dd}&w6_`sP=y#sVgGYtL5L!I3v#J5SlVkHs zsduz%{c%gHfh`-*UW-q?WZJ$I4LoJxs+@e{y-ko=XTBAH^1QN}Vh-$;&{x0oImG;KHN=k`#Z-Wc+QcO2)x<0D(N8iDCFLUbke68$s zx-mxW0(E;euCIlVx=3G^^3vH|y%5-`WLrFF${vpuwK$r`6O^927e|b3&1Y@BT;7e7Qm-0EP0vP3nE z3|SfLIlr5{nkSUm*=~BvkU4i|tf8`y0SKoxKAB*V3sTT)CR`rFg!V#ZYuZtMTN0hL zMnm>bsx{o%dPvwxao7_fY?8BxI>z=(0=TD;QywGlYb4b9tlH`Jcr-4Fy*8v$MCu$c znT3)kCSzFU2Nq*_o0e%6;R!|sTzzcQFp&rA(;6WMUfwh#Z>UpRvReYnUkrvXu{RGW z4-qOCXN|_ktEmm{0&^!{YMpKuosw8GfhQ zVW{m6`L0Uq+uQbzTgeB) zsY)6z9M>NHQVVHE5UZORS1ftXi00_nI3A}cuvZYIK>c?Rg>;gy5N@+*G`btQVWI+U%5g-%Dao%bVWz~?q zi|Wv@?oJwRZtwd9!!%vR?XsW}59yyX5XB+qx1$GZCP;qZ3g8IN(-`L3EElgE1)i-< z)fJ+X9QyGDb%$-BDj+Tw@fHddvj%+)kt#YDa&2c3t?#>w7*+469=jG2<-%*aA-5d$ zW3Dh2JiZ*cS--p-X{u1iLf3AakbkopQy2huuuvQu(Ok$0uai zPeEt6;R?hwvww1TFe{u^NiYvpU@Npy->0;-GulFqBu`IQ+qsvDm#&5VrP%ZHB~rCh z`B^#_RVax-t`jZg$PJnJL;D9ODu@bucwOetY766 zTBIbB{USkIYH4*f`zdGwqmFhpOg7EY+;_wf^*R;7gJ`<-!NkSC643yNUX+ks>f9R&z& zDi{Tr&5IN>M-|Y=#g?L;uo|Z^wsAZe6?ED_%qCryi7y54 z^)=>2!$<<;wB~4E-I}poq#nu4Sc&?DXs`gFO@8l3iF^0k!2&uCH71|QrC`iZHI=B3 z2WCJ(NZl^U15UH^@L&muy6R~B9Y|2g2OvG%6F+ZRq&vc+2McK?>~=p}!r)55jLl$d zq<`ZUV8%AON5VMYzTrPJXuZNl7Mv9N#74I09jhC{XgdPTZM?~dDnz(#=m&h@5R%18 zXA7EMlJA7i!w*=3OHJmha^*DUmMGQlyMk1_DEF=F83M|#mG5pn`XCYTP>|k_r1Pk0 zZKUMWFw!yMQb!~VO`crfehbKK{a%57BJW{TrSHoDWHX$tZtTm}@^NiK_(K1$C6d_Z zhrT9I6%3~9H>{OwR()EsKFjoErHK1dfduNs1ruLo0u`k3Dmo|x!OTF7!MkIn^4k*C`Ipon8p4f1P8G}@YSbEN z^!o*_`li2$;{kr$%9lGb(!iXh3KH>U;^sMG5Ae*{nJ_E-(ZFJbO)_jm#Bn6G^&M07 z=%5;?Zrr}&=%axkKDqjScXz*KY2Qc`>xedKwP??fd1quE$jyovx{P(Q2o1-5`Sq2* zSsm96g|xEfuPGyD1Ti;^PIv`~xP4Vl8D{bd$&n;!2xt(p7h$J zr}e_(-;7iG%Gd&(A@-bXI1*-o7t==sKy5a)374u`EtU;AQ)>{MN^~X-yQ6`-b(4q* zK5yOsBJ0uUY{s61%GNFhwneM~{fvO$u^|dSg4O*a+kTbC>f~w7zxSuq13RF=)-yn| zepdsn0?dUwKLg%!kOu{)!}pU~r?iU@*^fD_fJ;?FejK&)U(JlENBA54`5wTRb`!=q zk?}9I=J%oAWh5MK>`Li@3=b1EXy9a0SJz~=^5|Eq zMKuQ*<_jSl>*hC2_o@}xR1+%_)zJ%xb}Qp$SlF0-*mKIhN?DdY)GOTzUgbQsu+?Ha z35u0Q))LN#HA`hqx6J;^FAv#|))LQ+WAphP>?{34q_)_QY-~eU4tZqT9!fv-k?B{0 z>v-?;6eT)Z&s0^BBM$?{Q|YF@=~je@k5Md!xC29^bhH}#WH2}oe`rN6b|OdFx}Ql% zLx&qcRzXT9yAtkuhT+MKHSyE%!Nsy=ODo@!!L>OM;~HJtgai|YyKVRa&!WZAL}Eu= zgGT7B`;~%@WdWCZUTQdcAeL=f>L?aavg01sveNKRUedvDQgr!qK{u!F^R9TWXCHq^2WanK zzH(|(H)nkUWqtDIu)42I^IDYZJQ#TWiqAidG%~5{8Zo{DSc$gsyo;WMGO4Bbf;C(5 zj7}V69*&m(bOJf&#W*bcHkTbIGre8!d31PYkcC4U%|1fCL^{~3X;#FqvI0txj4}ek z6#hsif4ceePsQeNZErGvW~&Dr18sN5)vN;UQbSreGIcwYzDv^Nd}^+!Fk&nBl*pj| z+J!mOEK9J433t@h{{1AId*XC2h@wsqQnC$T;`HnQB#Um0cpFgSbUp}US$_ewYlpo!%m>x+pb?+R5T*Fx()rtD4*~R zL>8KQx)4B>*=z!UD6w0S>~p2lZYB#GP#=ZNsqenbtDn78&xtQL9v$-}cXT`5taeYY z=TLoXMw_*(LjUYl5a*4PXO80`J+X(fE(NGprRAKtFfwD@tk3bi=L1>BWx~gr6nVM9 zydG;f{ym9L`Wj_MAIIkTfR-Dqb6k;)>n7-|q8m-!nUmOmYT!uDXa*pEUD%270X4}> zbzhQ~^qVkdw5~kW*GhRke_%SF_vu<#iK%U@>i{iM+>2f}&<}+(0mq`1yxsL4rP%{I`tC?S6zwi$Fd@_Deb{ez;qdUZ8)x>(~rBg>kV zex)jAu(LXW?z3;-c%|N@2I)MclG)gRb6SM?nDFag%p&7Ur`^Y~ARXk3GRh|~um2*e z(lLoCP`t7WtY7KdzoRjK4f+1H-1IA1Ro?l1v^P-U& z{$o;~VHWQN1Cg z17CgCz{6N}Z}CT06sY8wM^+}PEs{`cC>e^ST^#WdhZ-i_iuStFoDf7}?Vkf?^yM<* zqsgYl{lUd1s&V1;6R_zNEH)CpSRWvD@88gYGfP;VF{{pXwnfFzW&rg_|jn_Wb zC`dvShSr;e4x&}*gwQkcNX*b@xvVOgYp_Sm$keAt9dG_Yb>CmQJjv8K_}RVQ^)AVL zYA_5**$F-8cdXUx}=oT2Tw*@GB`E_&MgR9K#XrhCoG>oe@{rw*f}#N zaJv!@!zwV9V1svlfx^Sc7%jb8u9+7Hyy2!B`5>QnTg~Mi;zQ#*wysOxG56ifHoY8E zr}i67H3Je^9I6RExia_1;EvR-Y}2+`F(sCOm<>$mjGU{~!04SeG&~d3J}S*VFw4|y zAPuqAqH20bT^N-ZLK4Vf=O{*IO9=Pc8~7Z;BxPg?%O_(kjO~!q!>|1mNnTG<+^zsG zY$R#V3y75!ELp9$LQ6_X1f2PBfSUoCaUk<0k>~G`nz?!(&+zA)Z6#(CB5mb_XYId4#fT4 zT%i81r2j3NuWzAiZ~tD_TgmuV~)2efEDV@t;9_UnkL-bzgTZVi0@ ziO_!r^nGgfn@`|(+Z6a4neXco{1xK+EaQ*AWxd`V3jb?y{|wyw^wmGXS&ROSwD)<+ ze+BnGVf0UM2Csy(Hvs=L$9|8S|4MuQXW-taa=j&){%(A)gx`NI?Jp^(?*YE$fW1$5 z`Ey`j%lvzh{|w}NPV!rt%~KgfLV(Y;OT_Z*%- zA9h0TZ_xdZyuZSGBmTVi{{0D0z~E1D|HtF^f7JGy&+z@C_D^`K4*y2pU*Wx3l=th9 zKdr!&=bz%!;_{B$oTW{rCR>f?J2c diff --git a/scripts/eval/spectra-cn.sh b/scripts/eval/spectra-cn.sh deleted file mode 100755 index d80b6a2..0000000 --- a/scripts/eval/spectra-cn.sh +++ /dev/null @@ -1,225 +0,0 @@ -#!/bin/bash - -echo "Usage: ./spectra-cn.sh [asm2.fasta] out-prefix" -echo -e "\t\t: Generated with meryl count from i.e. illumina wgs reads" -echo -e "\t\t: k-size used to build " -echo -e "\t\t: haplotype 1 assembly. gzipped or not" -echo -e "\t[asm2.fasta]\t: haplotype 2 assembly. gzipped or not" -echo -e "\t: output prefix. Required." -echo -echo -e " will be generated as output." -echo -e "\tWhen asm2 and out-prefix are given," -echo -e "\t\t and" -echo -e "\t\t.spectra-hap.hist" -echo -e "\twill be generated in addition." -echo - - -if [[ $# -lt 4 ]]; then - echo "No args provided. Exit." - exit -1 -fi - -if [ -z $scripts ]; then - echo "Setting \$scripts path" - scripts=`which meryl` - scripts=`dirname $scripts` - scripts="$scripts/../../scripts" - echo "Done!" -fi - -read=$1 -k=$2 -asm1_fa=$3 -asm2_fa=$4 -name=$5 - -if [ -z $name ]; then - name=$4 - asm2_fa="" -fi - -if [ -s $name ]; then - echo "$name already exists. Provide a different name." - exit -1 -fi - - -## Remove this line if R is properly installed ## -module load R # -################################################# - -asm1=`echo $asm1_fa | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` - -echo "Generate spectra-cn.hist per assemblies" -for asm_fa in $asm1_fa $asm2_fa # will generate only for asm1_fa if asm2_fa is empty -do - asm=`echo $asm_fa | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` - - # Generate meryl db per assembly - if [ ! -e $asm.meryl ]; then - echo - echo "Generate meryl db for $asm" - meryl count memory=31 k=$k output ${asm}.meryl $asm_fa - fi - - echo - echo "Collect read counts per asm copies" - hist=$name.asm.$asm.spectra-cn.hist - echo -e "Copies\tkmer_multiplicity\tCount" > $hist - echo " - meryl difference output read.k$k.$asm.0.meryl $read $asm.meryl" - meryl difference output read.k$k.$asm.0.meryl $read $asm.meryl - meryl histogram read.k$k.$asm.0.meryl | awk '{print "read-only\t"$0}' >> $hist - - for i in $(seq 1 4) - do - echo " - meryl intersect output read.k$k.$asm.$i.meryl $read [ equal-to $i ${asm}.meryl ]" - meryl intersect output read.k$k.$asm.$i.meryl $read [ equal-to $i ${asm}.meryl ] - meryl histogram read.k$k.$asm.$i.meryl | awk -v cn=$i '{print cn"\t"$0}' >> $hist - rm -r read.k$k.$asm.$i.meryl - done - echo " - meryl intersect output read.k$k.$asm.gt$i.meryl $read [ greater-than $i ${asm}.meryl ]" - meryl intersect output read.k$k.$asm.gt$i.meryl $read [ greater-than $i ${asm}.meryl ] - meryl histogram read.k$k.$asm.gt$i.meryl | awk -v cn=">$i" '{print cn"\t"$0}' >> $hist - rm -r read.k$k.$asm.gt$i.meryl - - echo " - meryl difference output $asm.0.meryl ${asm}.meryl $read" - meryl difference output $asm.0.meryl ${asm}.meryl $read - meryl statistics ${asm}.0.meryl | head -n4 | tail -n1 | awk -v asm=$asm '{print asm"\t"$1"\t"$2}' >> $name.asmonly.txt - - echo "Plot $hist" - echo " - $scripts/plot/plot_spectra_cn.R -f $hist -o $name.asm.$asm" - $scripts/plot/plot_spectra_cn.R -f $hist -o $name.asm.$asm - - ASM_ONLY=`meryl statistics ${asm}.0.meryl | head -n4 | tail -n1 | awk '{print $2}'` - TOTAL=`meryl statistics ${asm}.meryl | head -n4 | tail -n1 | awk '{print $2}'` - ERROR=`echo "$ASM_ONLY $TOTAL" | awk -v k=$k '{print (1-(1-$1/$2)^(1/k))}'` - QV=`echo "$ASM_ONLY $TOTAL" | awk -v k=$k '{print (-10*log(1-(1-$1/$2)^(1/k))/log(10))}'` - echo -e "$asm\t$ASM_ONLY\t$TOTAL\t$QV\t$ERROR" >> $name.qv -done - -if [[ "$asm2_fa" = "" ]]; then - echo "No asm2_fa given. Done." - - hist=$name.spectra-asm.hist - - # $asm1 only - meryl intersect output read.k$k.$asm1.meryl $read ${asm1}.meryl - - # Write output - echo -e "Assembly\tkmer_multiplicity\tCount" > $hist - meryl histogram read.k$k.$asm1.0.meryl | awk '{print "read-only\t"$0}' >> $hist - meryl histogram read.k$k.$asm1.meryl | awk -v hap="$asm1" '{print hap"\t"$0}' >> $hist - - echo " - Plot $hist" - echo " - $scripts/plot/plot_spectra_asm.R -f $hist -o $name.asm" - $scripts/plot/plot_spectra_asm.R -f $hist -o $name.asm - - echo "Clean up" - rm -r ${asm1}.0.meryl read.k$k.$asm1.0.meryl read.k$k.$asm1.meryl - - - exit 0 -fi -asm2=`echo $asm2_fa | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` -rm -r read.k$k.$asm2.0.meryl - -echo -echo "Generate spectra-cn.hist for combined $asm1_fa and $asm2_fa" -echo "\"Is my diploid assembly having k-mers in expected copy numbers?\"" -hist=$name.spectra-cn.hist - -echo " -Union-sum: ${asm1} + ${asm2} + shared kmer counts (asm)" -meryl union-sum output ${asm1}_${asm2}_union.meryl ${asm1}.meryl ${asm2}.meryl - -echo " -0-counts in the asm; only seen in the reads" -meryl difference output read.k$k.0.meryl $read ${asm1}_${asm2}_union.meryl -echo - -echo -e "Copies\tkmer_multiplicity\tCount" > $hist -meryl histogram read.k$k.0.meryl | awk '{print "read-only\t"$0}' >> $hist -for i in $(seq 1 4) -do - echo " - meryl intersect output read.k$k.$i.meryl $read [ equal-to $i ${asm1}_${asm2}_union.meryl ]" - meryl intersect output read.k$k.$i.meryl $read [ equal-to $i ${asm1}_${asm2}_union.meryl ] - meryl histogram read.k$k.$i.meryl | awk -v cn=$i '{print cn"\t"$0}' >> $hist - rm -r read.k$k.$i.meryl -done - -echo " -meryl intersect output read.k$k.gt$i.meryl $read [ greater-than $i ${asm1}_${asm2}_union.meryl ]" -meryl intersect output read.k$k.gt$i.meryl $read [ greater-than $i ${asm1}_${asm2}_union.meryl ] -meryl histogram read.k$k.gt$i.meryl | awk -v cn=">$i" '{print cn"\t"$0}' >> $hist -rm -r read.k$k.gt$i.meryl - -echo " -Plot $name.spectra-cn.hist" -echo " -$scripts/plot/plot_spectra_cn.R -f $hist -o $name" -$scripts/plot/plot_spectra_cn.R -f $hist -o $name - -echo " -## Count k-mers only seen in the assemblies, not in the reads ##" -echo " -0-counts in the read; only seen in the assembly" -meryl difference output ${asm1}_or_${asm2}.0.meryl ${asm1}_${asm2}_union.meryl $read -meryl intersect output ${asm1}_and_${asm2}.0.meryl $asm1.0.meryl $asm2.0.meryl - -meryl statistics ${asm1}_or_${asm2}.0.meryl | head -n4 | tail -n1 | awk '{print "both\t"$1"\t"$2}' >> $name.asmonly.txt -meryl statistics ${asm1}_and_${asm2}.0.meryl | head -n4 | tail -n1 | awk '{print "shared\t"$1"\t"$2}' >> $name.asmonly.txt - -echo " -Get QV" -ASM_ONLY=`meryl statistics ${asm1}_or_${asm2}.0.meryl | head -n4 | tail -n1 | awk '{print $2}'` -TOTAL=`meryl statistics ${asm1}_${asm2}_union.meryl | head -n4 | tail -n1 | awk '{print $2}'` -ERROR=`echo "$ASM_ONLY $TOTAL" | awk -v k=$k '{print (1-(1-$1/$2)^(1/k))}'` -QV=`echo "$ASM_ONLY $TOTAL" | awk -v k=$k '{print (-10*log(1-(1-$1/$2)^(1/k))/log(10))}'` -echo -e "Both\t$ASM_ONLY\t$TOTAL\t$QV\t$ERROR" >> $name.qv -rm -r ${asm1}_and_${asm2}.0.meryl ${asm1}_or_${asm2}.0.meryl ${asm1}_${asm2}_union.meryl ${asm1}.0.meryl ${asm2}.0.meryl - -echo " -## Generate spectra-asm.hist for combined $asm1_fa and $asm2_fa ##" -echo "\"Is the assembled distinct portion bigger in one of the two assemblies?\"" -hist=$name.spectra-asm.hist - -# Get ${asm1} / ${asm2} / shared kmers -meryl difference output ${asm2}_only.meryl ${asm2}.meryl ${asm1}.meryl -meryl difference output ${asm1}_only.meryl ${asm1}.meryl ${asm2}.meryl -meryl intersect output ${asm1}_shrd.meryl ${asm1}.meryl ${asm2}.meryl - -# $asm1 only -meryl intersect output read.k$k.$asm1.meryl $read ${asm1}_only.meryl - -# $asm2 only -meryl intersect output read.k$k.$asm2.meryl $read ${asm2}_only.meryl - -# shared ($asm1 and $asm2) -meryl intersect output read.k$k.shrd.meryl $read ${asm1}_shrd.meryl - -# Write output -echo -e "Assembly\tkmer_multiplicity\tCount" > $hist -meryl histogram read.k$k.0.meryl | awk '{print "read-only\t"$0}' >> $hist -meryl histogram read.k$k.$asm1.meryl | awk -v hap="$asm1-only" '{print hap"\t"$0}' >> $hist -meryl histogram read.k$k.$asm2.meryl | awk -v hap="$asm2-only" '{print hap"\t"$0}' >> $hist -meryl histogram read.k$k.shrd.meryl | awk -v hap="shared" '{print hap"\t"$0}' >> $hist - -echo " -Plot $hist" -echo " -$scripts/plot/plot_spectra_asm.R -f $hist -o $name.asm" -$scripts/plot/plot_spectra_asm.R -f $hist -o $name.asm - -echo " -Clean up" -rm -r read.k$k.0.meryl read.k$k.$asm1.meryl read.k$k.$asm2.meryl read.k$k.shrd.meryl ${asm1}_only.meryl ${asm2}_only.meryl ${asm1}_shrd.meryl - diff --git a/scripts/plot/plot_blob.R b/scripts/plot/plot_blob.R deleted file mode 100755 index 932f6d5..0000000 --- a/scripts/plot/plot_blob.R +++ /dev/null @@ -1,47 +0,0 @@ -#!/usr/bin/env Rscript - -require("argparse") -require("ggplot2") -require("scales") - -parser <- ArgumentParser(description = "Make blob plots.") -parser$add_argument("-f", "--file", type="character", help=".count tdf; with headers; ie. (required)", default=NULL) -parser$add_argument("-o", "--output", type="character", default="hapmers.blob", help="output prefix [default %(default)s]") -parser$add_argument("-x", "--xdim", type="double", default=6.5, help="width of output plot [default %(default)s]") -parser$add_argument("-y", "--ydim", type="double", default=6, help="height of output plot [default %(default)s]") -args <- parser$parse_args() - -blob_plot <- function(dat, out, w=6.5, h=6) { - - dat=read.table(dat, header=TRUE) - - max_total=max(max(dat[,3]), max(dat[,4])) * 1.01 - col_lab=names(dat)[1] - x_lab=names(dat)[3] - y_lab=names(dat)[4] - print(x_lab) - print(y_lab) - - ggplot(dat, aes(x=dat[,3], y=dat[,4], color=dat[,1], size=dat[,5])) + - geom_point(shape=16) + theme_bw() + - scale_color_brewer(palette = "Set1") + - scale_fill_brewer(palette = "Set1") + # Set1 -> Red / Blue. Set2 -> Green / Orange. - scale_x_continuous(labels=comma, limits = c(0, max_total)) + - scale_y_continuous(labels=comma, limits = c(0, max_total)) + - scale_size_continuous(labels=comma, range = c(1, 10), name = "Total k-mers") + - theme(legend.text = element_text(size=11), - legend.position = c(0.83,0.70), # Modify this if the legend is covering your favorite circle - legend.background = element_rect(size=0.5, linetype="solid", colour ="black"), - axis.title=element_text(size=14,face="bold"), - axis.text=element_text(size=12, face="bold")) + - guides( size = guide_legend(order = 1), - colour = guide_legend(override.aes = list(size=5), order = 2, title = col_lab)) + - xlab(x_lab) + ylab(y_lab) - - ggsave(file = paste(out, '.png', sep=""), height = h, width = w) - ggsave(file = paste(out, '.pdf', sep=""), height = 2*h, width = 2*w) - -} - -blob_plot(args$file, args$output, args$xdim, args$ydim) - diff --git a/scripts/plot/plot_block.sh b/scripts/plot/plot_block.sh deleted file mode 100755 index 348f762..0000000 --- a/scripts/plot/plot_block.sh +++ /dev/null @@ -1,58 +0,0 @@ -#!/bin/bash - -if [[ $# -lt 4 ]]; then - echo "Usage: ./plot_block.sh [include_gaps]" - echo " : generated with trio/phase_block.sh" - echo " : output prefix; automatically appends given and " - echo " : number of switches allowed in " - echo " : interval to be determined as short-range switch (bp)" - echo " : T for including gaps. Only set for phased assemblies" - echo "Arang Rhie, 2020-01-07. arrhie@gmail.com" - exit -1 -fi - -bed=$1 -out=$2 -num_switch=$3 -short_range=$4 -include_gaps=$5 - -out=$out.${num_switch}_$short_range - -if [ -z $scripts ]; then - echo "Setting \$scripts path" - scripts=`which meryl` - scripts=`dirname $scripts` - scripts="$scripts/../../scripts" - scripts=`readlink -e $scripts` - echo "Done!" -fi - -if [ -s $out.phased_block.bed 2> /dev/null ]; then - rm $out.phased_block.bed - rm $out.switch.bed 2> /dev/null -fi - -echo " -java -jar -Xmx1g $scripts/trio/bedMerToPhaseBlock.jar $bed $out $num_switch $short_range $include_gaps" -java -jar -Xmx1g $scripts/trio/bedMerToPhaseBlock.jar $bed $out $num_switch $short_range $include_gaps - -SWITCH_ERR=`awk -v swi=0 -v tot=0 '{swi+=$(NF-1); tot+=$NF} END { print swi"\t"tot"\t"((100.0*swi)/tot)"%" }' $out.phased_block.bed` -echo "Num. switches / Total markers / Switch error rate (%): $SWITCH_ERR" - -echo " -java -jar -Xmx1g $scripts/eval/bedCalcN50.jar $out.phased_block.bed | tail -n1 | awk -v out=$out -v swi=\"$SWITCH_ERR\" '{print out\"\t\"\$0\"\tswi}' - >> $out.phased_block.stats" -java -jar -Xmx1g $scripts/eval/bedCalcN50.jar $out.phased_block.bed | tail -n1 | awk -v out=$out -v swi="$SWITCH_ERR" '{print out"\t"$0"\t"swi}' - >> $out.phased_block.stats - -count=$out.phased_block.counts - -echo -e "Block\tRange\tMat\tPat\tSize" > $count -awk '{ swi=$(NF-1); tot=$NF; {if ($4=="mat") { mat=(tot-swi); pat=swi; } else if ($4=="pat") { mat=swi; pat=(tot-swi); }} {print $4"\t"$1"_"$2"_"$3"\t"mat"\t"pat"\t"($3-$2)}}' $out.phased_block.bed >> $count - -module load R - -echo " -Rscript $scripts/plot/plot_blob.R -f $count -o $out.phased_block.blob" -Rscript $scripts/plot/plot_blob.R -f $count -o $out.phased_block.blob - - diff --git a/scripts/plot/plot_spectra_asm.R b/scripts/plot/plot_spectra_asm.R deleted file mode 100755 index c0ab59c..0000000 --- a/scripts/plot/plot_spectra_asm.R +++ /dev/null @@ -1,53 +0,0 @@ -#!/usr/bin/env Rscript - -require("argparse") -require("ggplot2") -require("scales") - -parser <- ArgumentParser(description = "Make spectra-asm plots. Line and filled pectra-asm plots will be generated.") -parser$add_argument("-f", "--file", type="character", help=".spectra-hap.hist file (required)", default=NULL) -parser$add_argument("-o", "--output", type="character", help="output prefix (required)") -parser$add_argument("-x", "--xdim", type="double", default=6, help="width of plot [default %(default)s]") -parser$add_argument("-y", "--ydim", type="double", default=3, help="height of plot [default %(default)s]") -parser$add_argument("-m", "--xmax", type="integer", default=150, help="maximum limit for k-mer multiplicity. [default %(default)s]") -parser$add_argument("-c", "--ymax", type="integer", default=0, help="maximum limit for Count. If 0 is given, ymax will be auto detected by max Counts ignoring read-total and read-only. [default %(default)s]") -args <- parser$parse_args() - -spectra_asm_plot <- function(hist, name, w=6, h=3, x_max=150, y_max=0) { - - dat=read.table(hist, header = TRUE) - - if (y_max==0) { - y_max=max(dat[dat[,1]!="read-total" & dat[,1]!="read-only",]$Count) - y_max=y_max*1.5 - } - print(paste("y_max:", y_max, sep=" ")) - - - dat[,1]=factor(dat[,1], levels=unique(dat[,1]), ordered=TRUE) # Lock in the order - - ## Line graph - ggplot(data=dat, aes(x=kmer_multiplicity, y=Count, group=dat[,1], colour=dat[,1])) + - geom_line() + - theme_bw() + - scale_color_brewer(palette = "Set1") + - scale_y_continuous(labels=comma) + - coord_cartesian(xlim=c(0,x_max), ylim=c(0,y_max)) + - labs(fill="k-mer", colour="k-mer") - ggsave(file = paste(name, 'ln.png', sep = "."), height = h, width = w) - ggsave(file = paste(name, 'ln.pdf', sep = "."), height = h, width = w) - - ## Area under the curve filled - ggplot(data=dat, aes(x=kmer_multiplicity, y=Count)) + - geom_ribbon(aes(ymin=0, ymax=pmax(Count,0), fill=dat[,1], colour=dat[,1]), alpha=0.5, linetype=1) + - theme_bw() + - scale_color_brewer(palette = "Set1", direction=1) + # , breaks=rev(levels(dat$Assembly)) - scale_fill_brewer(palette = "Set1", direction=1) + # , breaks=rev(levels(dat$Assembly)) - scale_y_continuous(labels=comma) + - coord_cartesian(xlim=c(0,x_max), ylim=c(0,y_max)) + - labs(fill="k-mer", colour="k-mer") - ggsave(file = paste(name, 'fl.png', sep = "."), height = h, width = w) - ggsave(file = paste(name, 'fl.pdf', sep = "."), height = h, width = w) -} - -spectra_asm_plot(args$file, args$output, args$xdim, args$ydim, args$xmax, args$ymax) diff --git a/scripts/plot/plot_spectra_cn.R b/scripts/plot/plot_spectra_cn.R deleted file mode 100755 index 5de5f18..0000000 --- a/scripts/plot/plot_spectra_cn.R +++ /dev/null @@ -1,66 +0,0 @@ -#!/usr/bin/env Rscript - -require("argparse") -require("ggplot2") -require("scales") - -parser <- ArgumentParser(description = "Make spectra-cn plots. Line, filled, and stacked spectra-cn plots will be generated.") -parser$add_argument("-f", "--file", type="character", help=".spectra-cn.hist file (required)", default=NULL) -parser$add_argument("-o", "--output", type="character", help="output prefix (required)") -parser$add_argument("-x", "--xdim", type="double", default=6, help="width of plot [default %(default)s]") -parser$add_argument("-y", "--ydim", type="double", default=3, help="height of plot [default %(default)s]") -parser$add_argument("-m", "--max", type="integer", default=150, help="maximum limit for k-mer multiplicity [default %(default)s]") -args <- parser$parse_args() - -spectra_cn_plot <- function(hist, name, w=5, h=3, x_max=150) { - - dat=read.table(hist, header=TRUE) - - y_max=max(dat[dat[,1]!="read-only",]$Count) - y_max=y_max*1.5 - print(paste("y_max:", y_max, sep=" ")) - - if ( length(levels(dat[,1])) == 6 ) { - dat[,1]=factor(dat[,1], levels=levels(dat[,1])[c(6,2,3,4,5,1)]) - } else if ( length(levels(dat[,1])) == 4 ) { - dat[,1]=factor(dat[,1], levels=levels(dat[,1])[c(3,1,2,4)]) - } - - ## Line graph - ggplot(data=dat, aes(x=kmer_multiplicity, y=Count, group=dat[,1], colour=dat[,1])) + - geom_line() + - theme_bw() + - scale_color_brewer(palette = "Set1", name="k-mer") + - scale_y_continuous(labels=comma) + - coord_cartesian(xlim=c(0,x_max), ylim=c(0,y_max)) - ggsave(file = paste(name, 'spectra-cn.ln.png', sep = "."), height = h, width = w) - ggsave(file = paste(name, 'spectra-cn.ln.pdf', sep = "."), height = h, width = w) - - ## Area under the curve filled - ggplot(data=dat, aes(x=kmer_multiplicity, y=Count)) + - geom_ribbon(aes(ymin=0, ymax=pmax(Count,0), fill=dat[,1], colour=dat[,1]), alpha=0.5, linetype=1) + - theme_bw() + - scale_color_brewer(palette = "Set1", direction=1, name="k-mer") + # , breaks=rev(levels(dat$kmer)) - scale_fill_brewer(palette = "Set1", direction=1, name="k-mer") + # , breaks=rev(levels(dat$kmer)) - scale_y_continuous(labels=comma) + - coord_cartesian(xlim=c(0,x_max), ylim=c(0,y_max)) - - ggsave(file = paste(name, 'spectra-cn.fl.png', sep = "."), height = h, width = w) - ggsave(file = paste(name, 'spectra-cn.fl.pdf', sep = "."), height = h, width = w) - - ## Stacked - dat[,1]=factor(dat[,1], levels=rev(levels(dat[,1]))) #[c(4,3,2,1)]) - ggplot(data=dat, aes(x=kmer_multiplicity, y=Count, group=dat[,1], fill=dat[,1])) + - geom_area(size=.5, alpha=1) + - theme_bw() + - scale_color_brewer(palette = "Set1", direction=-1, name="k-mer", breaks=rev(levels(dat[,1]))) + - scale_fill_brewer(palette="Set1", direction=-1, name="k-mer", breaks=rev(levels(dat[,1]))) + - scale_y_continuous(labels=comma) + - coord_cartesian(xlim=c(0,x_max), ylim=c(0,y_max)) # xlim(0, x_max) + ylim(0,y_max) - - ggsave(file = paste(name, 'spectra-cn.stacked.png', sep = "."), height = h, width = w) - ggsave(file = paste(name, 'spectra-cn.stacked.pdf', sep = "."), height = h, width = w) -} - -spectra_cn_plot(args$file, args$output, args$xdim, args$ydim, args$max) - diff --git a/scripts/trio/bedMerToPhaseBlock.jar b/scripts/trio/bedMerToPhaseBlock.jar deleted file mode 100755 index a86fe09fa3b69342ab9a2b533d3c66d463ea073f..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 31354 zcmb5VQ;;r9(5>0FZQHhO+t%B*ZQHwhw{6?DZQI74@1L1CapGJ|#8gy8W#mOwUFNfL zm69wd7#a{16ckWTETt0A|2EKoz<}h%)r9G#6(kttg%zYF#8uT8Ep( z=@}N_W$9^Vr{|hfn3h@hPafr`r13v5*J#J3=w`td$;Goi zyn8Fei9djalGMs_zmIJV7z`MIlw={H)SH;F6#lzi@xP_?zuvS}l=eS~|1}8yi)@VD zjXfESosI1+8BNU17~Nd0Z5c((%@|B=ja^)NbfA6JR#yJyXK`jt9ES-V2AZNk1tT>g z#o`dt;h3aap(Fw&jZ7xSnz6yo%Ax61lIWmkRI0tI?#cZ)1Ou4O@649aQH(a2PfNVvTmZ`u0keoRoLm zD5Yt*EZevryh`pClTjWI)eTQyRzbmyO%I#FR^F<)Z!MtPGEQ9OHgU)Uh*LCTu;(gr zGM@6v2J>o$34TR_h6ZU6132(u-|TP02rKwlcC%yG#f*{Sf)!BE*9Zg>*hG^jH-T@| zG_>$1nHloYWTsosK&U)5lf^so0utAQDtX~FjEzDC1ziG@E3@>saTo8Rhuqo z(J4kvatW6zg5cj@pf%tO!U!+5fjs`7F|Q8qK2PSq~R@U%WLAwNS;yL`jGRE~Z9J z6q(75df8ZKVN6j{e0j2<)b0$0jw8lV#F&iOS=IYPE+M8I-=RN*5#7Ri6IUn4>J}-M zF3!cDY~`%8kea5kUs5(83yAlnJZ0CH@e&D#r7H5ekL@?@@EJ9iGC#zbs3O z#-bcc8^C%p>nq(1nyIx&XlWer-3IjCi^J;e*%(rVCvfttt^{);p+k-`$2mDZMgwfpfH5$lJBwMn)l+r ziN*A!07j!)SqUgBegGylB6Nm@9Oq`8=7Bq^$G$4shslXM?3|`BQ}CO{Ap@$%(UKnC z)kS2jW{hao>{9uvz4}?NWR2CdV?WR>DjMyWRWshu&1f9{mF#Ic6$xUdr_xARQlvf( zyitzP;JD9Z!WjdPJntM~(T^`QEk*Zc*kgD#g~2*0+!+zUK#{+1-g81AD}}oiXF&@u zngGN-mY&B#FmzhN>(_(& z=j44N&%bY~rd3&+6QCjp<|D|Do=N}@84HY7H4G0-RJG0G3 zpj{-Pjo$z{U&uZ3Y!rqZpCe`SjnGk}Kg_+BhhLd8WU}?8Rf_1k5}29`d&xbR5^2e> zEbAdvYDoi+ZpBNBhs=VzDgp|J(OZmi-zrX1>Y$tY7X{Gf)ot?MIDN|2B-kt3hT>IS zDz?pImS_7O?}mc@bKA_mtHXr%lH}U1e$G_W+(2qk9&+E<0mjXgJ~x9c&93GUf|E?A zV_jt3ETpmD`2U`d2rGo7i%^yZ=k6KZ3gzU|siP3GJr#{*H25}me4NwIHj=GlF#^uHi#>_cj0U8)sg)c>5GpDS=%b0AZ$ zAkci3@f7JzU8*VI+mTkdX&YUheL!Mr;7Gb`K`)n5Bh`Gy*r#7&}XU9mcxruL5Us;phVkN0iFy)XMtHvw8r>1~3 zH5KA!wiS=3nqcSCpm{wa3#p?!sE80HDXR$v6?5Yl0=7c5H7K%Jv`^+ALR_d&hySw4 z+@(;5%4Aidno&T)kF#Cgy#Gy=Eit4Sr|&7}PQGR{?^Ih*?joJ$?_HGm6X*fs0Ouj~ zQqPyA0$lOOYEpW2Rx2U-$V}A&rVeCROwBS}ZV`^qB#p93KC9@eW4G5#$a4?-gl0y8 zJ8VMre5$Z=uD1N~D6iMPhVgU`?__L+T5?v_RkE|Q@WGX1rUb-GyB;G^hU2TEKZ}Fq z8?g&)?yDXk`RLoHUXWrZ6lFGYww)ld;x;ntzHJ)HC!5fffy_;pzDr7Fu>oD(bjj50|D=3pkzV;HBhMvyj-`=W13wciPH2Ci*A z$;~sVa(Hi>j>o{$$s`39>ZX%4m47W6by4X``==DP(}wTwuh$skwqfHV)Expa_g*dX z5?afj(k1RS3SRo;!T#t9*9*X1w8SBe5Kr2ZdL zSFow-16N}pm*EgntuB+HWJvy~;%+c;5H z-T9};M&EuCodkDNw+zb)S4)B_(!Ve=El;#F}$HQZCb7D zQ+=p}&tdZawOn*#9|B{-jgPCDX&YyH^ zB`Gb`I`1P}c;$6E?D~>J7IYze_|9LtP5V|s+4yqlY&5T4f$AW)XQSiub$FxfW?E1* zm*6g^$zoJknFMr93`f|w{#3Q50Jc8o{t|M1*8WNH)lvIP_4=guvv8f@^zyX9*smON z9-;R`^e=vdYeyj4CD_!pD={EO>H2x>=?*NF>nG~l=Q{D1Of!Jx9|5P8KgyOWQ`*1u z;5Hh1(tF2^rA&>}ow&!wo^rikwn}$Xh@6$FxE`*gpX@THNBd=GD+He>_oGzTCWf7A zMNZt`q2s)HNUo57vt`~vELeQr9Eh%H3EXL8)u#iejV*za0CDbFVen0DLc`dLb{zF~ zXjV^UL4vJzXJG2g)rEyf)|aU_v|RHFa7870f?tMz4&^1jK=KzusSB_DQIcHK(tL7p z7~a+7ne@#z^ur|OI{62Hm;m?+>~DNS4Fk5oRj-k^!LuD@H4Ov4z*gYuZHWAap)S#$ z69D$;;#-jFMU30?<;QcHx3_-CsiuXp^Lv~m7@+xZRJB7!-v84MQrAF?`!eE;|9Fy$ zftwIOA4JtqfTw?YBXVZ&OgCrv#&)%8oYajD;E2 zNFbpwPUWw4F+So3^cLYedU?}6O%Efv2p1qb~W6< zNSn&ss+FJ>Zf4b;mSFuemgAZ3!8!%^_!~iBDm|gQFU-HSu^Bfu_zfvNx?9p765Feg zAj9rWPYede^=U_ggU9L&UqMT&-Zdz_dRXVC-Kh5{_@fmz+MA$__yz<1xc}+LGkI8W zO)u$R7%*Qixc1MwMpB;-kYbzK?bK%JXPas;cG~VG@bykKO4*YyZcdBO!=ST!{%8i| zeHw@TM*M1zr?ru^oMN0mV5V$ayHe=n27TN+?jV!ZW9m zQnY$zM+UvG{2d6Rv2nKaPP|KUwjvfzxDW<694RZf0HJA-7ioLvgLN}D`<8e%rdWor zlCN!$JP10@6Ik>M9tJ8XTY{VXX4Mr};6K$icq|-XV1ZqxZL3O&_1LgA0 zJ?I5XceC6;(eSlgk{H*V;3QERalEaebK=*oH9$8t2=cy4N9X#c{x0Ud#z@e8)y7d- zLOEV5r=)QG;$x%A!yEX{MN2&ge-_3PabNN4L9?~#ZQp^E`FR!WtD&7Q7azu^O$FEO zX*N4Ix?V6fB{=eIi)IM8t4pO{wLralJi3`@x@Q2y{3RvLX1kF&n!vW+4E1#{!h} z*GA$UKK4fU0jhQX3C{~n(HPQztk)Mm;b3bxlhFCs?sP`b6!YqYawOr)49*>}R8vH3 zjlSgy3>{q{o4asG93Sd7;Ud(y{-~19>E?=?dTy=)y|j{Ducy1PoagmMvs#?6}nyi-iq774I;*`O_P6I>cAT(+# zTx)}jd0c)fV_D6@zG_lXhSwAJB8ejffae(`pkC%y0oA4$L% zK6#ATRDp%KTV@*)Xd$}MvSe`oX6d?##jPdjto)g!kga{Z__KP`m4|Ofpb|KgsFc* z15t5iqxq2ZO6Ow>8{azGtXu9aQCHS+^J2Sk^^$2&HzgWy(iV3fN*zh^OR}kZyap{`^^G$H@`chpq-yLpcMeE;7~Uw*z~iOBx(*sQI9>Z zB`Gpc+W)lSseY)$72t)FxAU`Ym2_?r@xoUqAOy{? zla8%s5-&Xw#lw6I0@tY;*^N0PIE+aXFKN}(fBf2IVz5B9mNuOd$ExNm!>sj2c}oWJ z5{OOJ(9t_2No^YJ#H*(WI1gJ`~Ody{N>%IRQLQS|!)@i8|I<2IY`Q$isu_yg$Ux+_A z`!kFL1sQH0AIR%cFMyNI-M9lAeDD{B7!o<;b^=c(D{Xg6iC`8REVH}6x&DU(_fxsC$jn)ICsjN5X zhKL76l7MuxDP=l588XSz7-|>u)ZCe*Og#DwDXrUd(_<&<)U#r#SJ{Scss8T@JlY8J zWXo8&*VIX;Je{mVZc!P2OPLPKC^b(Mm5UCMQXl&inmK-!p(fe*w{EhPytNasxrFxH z#N6YV2s9(~hqy~HoO`Qz5R`LZF# zQcVXL7^cm1Y=e#CnIx9DV^bzY8ai>TW>F1V={U6JSqsE=5_P$z)fp$F?2Uu%;1k35Od1o^pNpJRS`=lMy!{oUA{9#+z1~IAr%!7MX{pO~ zG0h5hJtjTHFJ2BS6_k;2yKyRYhD#hioa5Hp#yG~$+PtwkbPeuUv`)b-Ls|(9s|bHu z)%5>zGl_2Dt;;x3$vN$7#fxj!G@VzvaqJAD>9;6Q_L+~vV&X}lbaZNMS^pUSTr)h- zRAi;ZSmya6O(7n3hP3~zhnx#UhSH#bioR7m-tm+!L0~(q!M!Y($;^0#?4gHmdy`*! zpHH|r`{vob}GWgu4Kk>Ye!(5c0*`S21OJ5TXo{4SEAz z!(;l0My&YM7CbGha@GrVh!6^?VTA_lp+ZUCp}aY)m7g?@f0j&oT+W?1JM*|<@Fc(y zia>Ibq-Ytbh9(fXAC(?C{j45kS(59%j2=Dwv>q=5(3>af8u9}vy2cUxpHvxrvp)FG z4>>mu{_p8uJ6ZZf-}1pPO0mXCZ%YC2fUhy#-Q5s)A@g{xP z(rYwV55IwE!$jB-NObz@z;69wu?wIWa!*QlNVd1DTfJN2YkL7=g}d)P-XyH2zf}I` zkmM&AZr2mHcL{dWJOnB->Rkt6itI zwdo3a9y>T4`S@oVsFG9}_?$Wh8_U6BoZMUL05#mw>JFc^D$_2NiW7oJ8YuIlyw@7mq4Rh&Cf71-U?8Dnp zI}Si$gTj$jUIpV58y!+-+}RvTH&OLnZ z2+Cs2{E!RYiM_PWRjLhm^?-_*k3AbF1At?#$lN3Uqx4|er}rXf@jg}}ugZ{HRz!rF zpoWUw#h4jn_FbbOB*u7|Dz!(EWrA&uoaN@jq$QvDy>~S8BrAdP9GI)OO!-(;`vB{N zdZnHFbK6D*G{@ABPJn@4U7wxS9!=9q+oVGJ75VZI=;i_8x8VOM{Z9@RW_g%EKt}BU zyV7qJEhqQCWd8pG{Qobd|9@HlC%6EYC8@`zipc=1Zgry+EwAgH2>HYMM0ndvDI?VN zFQ5uq;{-eW#^V$-sJ4lKzFV(8vMFL){T*dAt)kn!yxTngSnt++z70>IJfv?+pQ8a1 z0>9|8)Nj^bJVQw*WXq*_z3YW9|EwIwx-47|a4qhiEd!EIf;cBvJy--f~7+uD6SH1 zfwF~%>hUu#XhGdz{Ie(U?#f4t&e$Ul6;E8m*0A<~_>2`g#&K3K3e2X^!+N#BN7L-s zajLkfWWF-f$Fb26gm3gCb9Bx7J+g_zR5ESx_<$ZyKZ$U&v%kRSk2g&*Z9vUiV6M); zuz}a?LAibYJ|Ikj90R#O5J&;mr>BS5IahU%F92Fgs+E|OjIz49f5abk`o-Wgn2oei zm=kC1g~z6OWo)xjDh42QQ0`I`PCARApb{HT2j=Qz8_$uk7U17zM}dh8i~@3X4$1X& z&q`Zxq+wRUa`bu*9x=e7(=}^trzp7MOpWIC?%C?DGBHJ@XsXI-CW{zQJuQ}hH{N;Y zF$eicv<((kp2^=ouUS8Rt{8;2h_4q6@^2Ou#+!L+=*wLKAQn4vR(U;?&CeKMF1O_` z+wzWEVklEY%7M+4ReYgIAu!qv9UJi-`^(5HEK=oC=GaJMy5*SAf&eB(C&BT_(jYo? z=V+t(&bj%@gD%^{tE@(-_sdhw!YBy>un}`SWDb(n}LW|KwuW4CoTy!K|2mm0G|&OUa8`YrqTjx2@d=GVh#{ zNVasR-$42nwGOOPPs!=4WF?%XN}G8&B)s8ZnKLJop)a(dBFxJ6kgJwoDR!Ev&Ic0ezmAJDl)$AK|#ioVxz$=W_+S&c|TRQF6w z5$t6`Egq#jw{d+iO)NdIE4?HYlHo30t)a>iUsXuYfyo3*JKpK4*6(pJ`i4^w=Q$I9 z0*p-DB*Mh9pk48@bY&eO6jLFLYV^1xL|xbu2sYj$L`DA@1DeJ+f~ishR2?5Afqy;7 z<;M^})Ibi+z!2oVUn`j$*Oo?3Yv;pdolr4Mt&bIGI*JtdN*rsbR@k_}Zq^JJqdTOS z6V}^mLjr2`cteA>m`|5aST-bKFv3PUN>ax2XCWD? za1yn1V}944H{NB8?|M`6GtmMuAZA2#vngM=aP54-CnJ2dG{{q2O~>lVQT1VZe$tiL zANTJOhzZhg$V^!lUcllMkmgEXW)ZNOaKTyEsNVcC8RU@zp=7-8$=A2niNwX7Fwf_r&k>MWMeD5A8#)Xe!Vq@g&;pBWDmXwQnMGjk4enF-p^!!f_9#KbHrmQi1iMf`MhC zxS-)2oa{KVdZbbl{kNGQK=O`BR2aiJCjms*NrhS<8_AK?ohu1Sy@Q}RWs@Jy%r$q8 z7sxQXij6LmlRNf%FZJZ9^kaGngG_C0;%XBa>P+&n2@3)*{09O z0x-~)bIMK_KoFfG4cBdla_q`V*7UPmEYa5~bp)DrB(EMu8z<#R7P?3&Fx@&8j3gxo zFkA;xvv1hIr5L1mA_S(bBbe)O|V@(1XW4l#(8aK z+~~BKlIb|?IN5FoQ)cL=@sikDY{eADbG%6xLKPq!FioNCDm#dtl~+B9n3D2xc_&z* z);T<&)I`r2&a|hzi9!c&VcIbt4hwCUcxhJ#MRZfC-u0{nK7;kF^c-GRFaX$n@I+?3 zwFM5Xe@w6e-^C}f=VK`i@}G0pwoaNJ>KPq)vOcLJWN|2b5(jeZOlrB%$zVHh-~u-3 zBZp^N7puKfJNYSLQjB)wDa(VPxy`;Nht*Mn;`l?%DYmBMG=TQO*Q!ox8#bGJmpLcv zk`Bmhx~k`Dxbsu`i=6(?=nXBLr0RUt9*aVG|Kqzw9V;JiFOBAU`Ksezdl;2))@Dmb z{UnHIur_-BBv#k<6zvWJI`-5R<+LT(eJ&y zb_fdLdVXB%x#7ddQxHo@u`s7n9InwI%%BuJgnDxEZDsOyNtZ}Vl>~yIQ-m8nIIqnM zRxxLnBW4$}`)wzUP=NBJbJEfU5ePR93a_<`)idA&T25n**qNo-8p>kR$+&nn^j`-y z90aG$v_FZLaUb*{*qGk5w-u)vM0&N`?QkPe_48!y5CVvYM;FEWE{*5AI6S-NS zv5z%_t}o7|sCh2YL=pF1XJG|_92_q);AlX$-z5Tc61|oG?(4{_F zQ5eBEW$xlm({ZAI$oI+xhD+*90d{2wD!3gnSCTd4VXT(Fh!wV=wwT`9ha+hR^P#c< z62Q!mOj1s;#8j6j@*qFx(y{?l)2r4N z7k1)wqT#*v!oMWB7rok)(;7}q2N7VezfRO)P+jsWOwA;=D%k(k7@|=2O{o^f0Z%&h z{7><;$nd-LYYS#TEjzjAtm?ESMymB=-q5BOn-U6p<}RF+5$Zp=QtB`JbBD=vO~;0^ zI-)9d5x()NK&te4PwhA}nI(*D{>Z-GO*88ty^r$#ou_oqdUNEFN9`wRl3qmA z3KN8r70Higz2&_)hEVY2arb2n3Bz7xJpIZ`vanNUkNW-Jk6t|oQ9{ftnlAtGH#lNQ za~;HyMvk$eTS@q4o&{^@Yo&b*n?0~~!6!uTBT*&jdJ$4>l~!&2njYj$M_A(7QvMkE zOC4=(=TL;Ar+9H{#$e0WKn{(ly1OVtKLcswVaU%a0Tlg4I}J7E$QH|F?mF4yR<0lU z#h1aqko`&&YCt+^4zrir>R3Oug_tflFv zqOt$mKa4Ea3o=_((9#Sm&Mm&$3N6f9c!%CFo%W9)8~+Ky_Ed#>Z_EJ=ta*X#jr z81UM<1h@y7hyFWyX_m_4l__M}#=U;!=VZ7GVa|v6_YL_!bonlv7>9}g1SBW^ztd&& z|KV+8_&?~<(aPAx`aeIJtFyHOqrADZnuF4R8FLX^2UEcR<135P@^V#M!Tz_y$(rRg zK{l8eggx{)9ELOwR;a>4haM}%79AcgaZwFjP2g66A zSJkX^-sxfj0Dy#mm*?l_gjcn%t~YrRj%B>lox9zy9S6hPlM_ZndrH+wn z=C+kgY}OU+lv|c2r}z%J?ON+b7Ql@~t=1H(d-WDsroUU0}jGM5rN~V>IOf?+s%N z-uXesZM~ zTaM{WI{s|rro9*=p(~B3hASdIPb~bCJwYQGo}Q(r7F)tK6LSY6hhw>wV7~ag+pH&3 z?8~K3)J>N&LrGmt9Oux2Yrz2^-$L<2mcD}XGOIM&;8EFNk^wjPHNoK2Q<8*wd3lUho!C`~!dWQ)fIm(g{$sP3AOLmXAe zu=ALD7c6LWQnPPn)-hqHJ%k3tnYZd<^Nhr}fRz%W*G-U7GwQ&2apb-{`tk|mEdIIr z)=IcrP~x{Baj3aP;;AXDw8? z2VFjZ0H#%hS;Tjqi{P0@WoO+vvhvf8WENvJ2&q1YBac#NCaG1>2s^=)p05S}{zfMW zO~S55Vec6k3fhk3G%Pcq?nTs01LjXE2JiPt89#w?b7^ z*8Jwd)?R#60l`KaXAVL1p>QF>FZE$~yc>Zr4-Yv|A$)2bFR^4lFt+RSK>IZGhy?pX z@xuv9afzTh)U{1$hl{xLK{pMNFi#oe&$e=`)b4=%;ktZnDwhhi!yZfW8ot2m@MNO2f~#|QfDux z^z`OUQfppX*+~^sW>7O19eqWv0$vVF*xG%y4@Vaf1A%lLTu9xWWb_+ zD`nW^u7v@n+^Nx&D;8({bSIdx%s=Cls@o+B0^Os|MOW^VEY}EOi^h~F?u>7R5d0>~ zSk^D#>pz{EGsFtnX;U@v6dH%M!2Kn7U-n0s?8~Iy5kbx_ zNkzzQ*s`_`Cw!fyaI%&xXe=rKV_Lx)?#--Py zQ<>aUgL2kISpe>06hE2#B+UADa>ajXcM;|);->PNxzU@ zf4b2rNxDNCL|WHB1?%2yj!z8Y$ZhBjD+x=6 z8IP{n;vehy9H;%|TyMt6-u%HR4@-VBAHmkxLqFgX89aAEKh&V^U9rY0_3wR0kXmr! zlk>_hfl3Lu4T>?`+=lOnjCd@{*+rvfhfHYD6McLt5q$0jb7?&2$txe135MmO)zpuH zRP(~E_yZ?tG;e?`-b?{2#(Vlq3E}!DHSlVPe=RsC$vkNW^tt>Adszd?zai2Q{~E;8 z<|b7w2P{$%l!Sc%X~>l5$yBIk4&SNKaK_2VOD4j|m+`oi&r@ecQSH%}9f7eOitXMg zh+tX}$H~|Z+3Lm8>>8XN@~9n2h~2AlM(7z7)1kwrDg|j!O?#<8@WNx<2#3Lg}i~X$|=V|ql#WQN3z`+^q@B`d`UGIY{b z!lWZ|(X}T#n}kZq?oG+v7UCyNYSFl^88RPbCh-{X(Y|fSh14@CVJKxmd6y^X`6?n2 z7(B_hj7Auyai#A&lfbA=&p}GCPOipq4X1hteIm5>6zxw2qIalD1%mw73m%pJHZkDA z-H6^3q=pZEF`Yba-S&U(9~+hX4VZ0~ScU3~0FNOIBz^@8_wlw?Vo{I47nPLnBBLy? z@#8#{CPgkr#OCR`|GRSx>JS+ArRozy$$xuZ79OiCUsy+fx#>a?(V8+nC$9Ldtjkv87-=2IO0$-ycK8Qb%)uKT< z)p$azMJ0ULzN2wQ&b@2B<8ekYIyRk-)aWp~Tc51ggw;#tg!48ZIhMJD+^(MP3)n#u z%wI0l)eMR1&=Rb293az&OST~Jm5wxM${sq=$5A!|>J_^PAe(a>mAcRst-BnhaQj&| ziK|!hhFmtcJPlR~JhW%D3}D1{M7pe44H_c-c&pOFspf8pj|WROMYacM&PUaC;BSjZ z__Zo`g1g!^LYF4)?6-$_T2;EDIr<{V`;@!z-kj>dW;V^+^AiaHEwcAO1uD%pn1%VC z=3ItTN0y-!L=qQ-WN1Pet8o!fX+^qMgF|$oDM^cTd?Hpbtu}U^sY+1%%>$;xutkbn zfj!q==#_c9@%@@`5q~iUV?r@!zA3A&`7&Y%&+=meM4&JM51QTS3H|C=354hY-Xq1P zXI89~X><9=2tSwrh3{yL&(>)x4FPM|nrA~a?bfRs7c-g}dNUNjjaLxm_MA~b72s9) zyZeL)VY1}lig>O(uHFOc1QS7KzKyl~P7Z(g;J{OO56J&Sh=-U`JHhqAWaZYdpWIU! zd;-U!cv62f`{zwx|($UeVgL!Ul| z17;urs(1qb{_h4DE9pNSNjM-N9gP1kM4|uxgDC&g5S!7qQ#_JF`Mv&*9$dUYp}=C{ z^;irG(%0I@u*Q-IUGM`6<&mbm7;*o)=3wNV3fQZ(_vreQX{PaYpb$KKo8etoRz-`Q z!&lE8VYor+6R$P9=#@&ZMSQP&B691&Lo#7ji$gohSX5P(YmPnlwwFmr8Yk^aVd6hP z&}+GYZn-3XRr~(^%ni_IJzXkSBZ-gyKnt6lbfTkRiqcl@2cN4u!Ca0e9CLmI`qqAM z)d^fV&@DqX>+m#-7A?pej3ZSO_ygIfX6P({IyV3u`i^p2SB)o&?UD{r!8gn81Lb%U zI|4qI_Q%M>$EORmpTn0S1?m_snzDcb)&MBR^kcZ}FIg;iXc0IKn#K~5iV+iqHfawp zRzEP@HJ!aOECPz4vsOJ_`X%{MU<#}ebtKt=Y5_Uh)!4ZyNmCa|my|jZ+J8#8=m4aA z!JP804GMzs{F?0_2&*|6Z|Vt3B?DhQ{L4WCW7=Sn7cs>8 zDNADG=l2hPEvm$VU5as%*NXBZsp5y38!#@&7|N{1&{N{*4u;bGP_y<8!u0YNZ7^+F zk|>a@@SC$Lk*ZetWthW4qnZICXQv2S0*v1GHK<@F5c=LRCV8_9vXGGX1Ji6IvFY1z zx3AXQPm^Ai@408nqmRKyF%)i3zqOi-&5bRLKLtA3k7Xx%=dXnwY*nw*m74olJ@X`w zg1?X?y3FAW2wXt@?wuNsPvjH_ULFeLYydE`jMuP5^cPBrlZ}KE3;svne^2RocUir zEqC5-Zs4=G82M)51j&A)-`|!_NTW$gx4~QE&%;t9cXx398mHqv9p;OA^d~ZXgd>pG z+LvAhgYODsW!X~Jb68j@&4btz8QFR{sH`Rp2&T9}aaXiv%8GX&t5H{UO_RP#I1l2yvw zV)CqQiT9<^i_M(}oY~BURYUdN*F}o45rcR1Jrt-bwDdxr10$;Eo(ENVyc!SCGih9r7Li`j%g#3DpPE>i=>N-FPH28nIM8u z?EX7j7QdQSeP}jNPzDZx?`r84D(p_|@(UeOajwHvFkd*tpTNXh#->SL4b7Zeh_0M8 zBH|`rf}&4|ubHR`9ZA};g%#Hg8p`5(R@@Go5AQnAhHI}DwFo?mtjj?gfzIdJ{> znk40$v=00{oxD^SEBv#9=eq{%bBElf&S_y+84#Wdj6PU~I3k8{DoNX$^YniR)KM5DWHUD$9u# zL%}5^tcHqbCS9W)(0~8xp706L6))QViRxD35&F%~Z8Hnt=~KFMvoMg5cVwk)Z)+Qk zFmPon+7#r|3#~4;8BKi+LqwP%( z4B;QR(HQnE?__FHX(}p&Xx-y)3wLk>>5lM%nXXqM=W*VkYc4c5>!!S+uU2y=cR<#Z zor2!Xtzn@jP}4LeqmiwQ-;<^Uo7GLM)Fnpmi|b_&p?#9p|6r^4O&1dW07%C5ttg96 zQtFb#Wk$F|(2xP?y0%-Zn~9X*|H34xV^MyBs|d|<9-)MxecptuH@^maA2EpO-WVUn zUaT20JlVA6S_84pel)mGuj1xeFMAvKGj4%V$B2*?Q`gNmGvk+$tb{;KSpBtLF=YY_ zM*qMZa-kcIRk7FP+l1$nvD{a2Us(ianL7}8V|q{5Q_PM2bZg{gT44^Q-k#S^H=UfR z5J2$K7g%0`>DOh(LmN0@XmzQt*qOBInADVsKV&GvC=w>YkZaNTAAPRdmM?N$Y0~d` zUvkpd+W<|2s&Xa#nu4WeC)DeyO)ixMuX5s_dOxN)vv=dJul+wE8T?_EfnMrrVcU z2$fZ%ixdshJ-}6RW5py>!}nNGxbT^E9qbkMXs>g|`s8R72<_saSkJ^0T1mF?#*h=7 z_AYn92l)m6dRar_#>~LG$*YcD?OVw<>sR(-*v>|DiX-s$HR9dbFO?W5i476L#tL@N znm7LSpCtAjX6l%89RoL1+Vsgmf3HkAF?e%BJQ#J~Mwq>oBP^9q)NPqIfgA~9Ab zEShq=%Qkrf^;b00qoAHGtXNyqs< zQh6qb6vUb!hr=OY;)GNLeY;cP+fJ&r55De#IL@$HIeW)#U~Pe)WNTnWgh;G=5BD2j z=}_PFHV0MjIna=Q3XST-mLzQzUos0gL;N3+7lW(;@eB$GXddr>ub}?V&%$K?M+GIV z$Y}DP%xlUhVQp(JZwxSZ{(q%i%`X?+74+ZL^ybFyhPXq^Yp|^JeN7S==pZh6A_Y0P z{~RWVMxgYK&8?hCeo7e}cG5CtZCMx1Ly9-goEDRa<~&Om zuB9Ka|81%MXSd$uyuse1i@W@1+b+F6y)Lt9|K9iZzA^f5-i0Ky!s*9~SQk1R*-+&b zLcbt}Oc$!rVc&G1n*_ZFr_b-L%8_#%Md%>L;&U2u!y6m4;p31?>QxSYmX3}rMq>%5 zw76N9Y6|%5&Dp?6!~TeN#HQu`ofK3!16`n24RJ_FYj*Zw-(Ga@IeAHbHHOd2Q`?67 zi>D3&Qf^V)`g#Ck)jkl2&{xfYQyyq2`3J%j{4X9=pbp!3tpvOqIY=vF%m@YOj1_Tf z1uCDY1~}GAgShC>CjM7hjo9`y7}B6TkB%n!WGmA>xUX`E0IIxPPRJ2XE5=w$%EpdA ze<~11(LjPOo0YM48}7LfrLg#tH7oX7eGS_0l^BB1iG5Dc4HWB7M{3Qo(v}8`D}Aro zK1Uj;%k1Q7Bl>y}uj{p*U;l^gY%WVZtOtvvMY!u)$BScTt+T6@L|Nz0dd&bWSm^K@ z0jMDc-e+i53eCMK;ZJxWG|t$WsQd5@@OB6R{?~*BbUp+ERFt=lSX41#4~`B(vB*=! z#e8+zAGi~hf@L^nZ6_B)t73qVg5$Ul;jQU~&V)(1-^BAX%rX!O*8GAYT8osO5mJP9 z5QtiuxFHbYVPt#(x=&hWo-~JQtR2xYp?OC03tHKRy2sS5DFne_WOy{z$dC&8q)Wm% z0mnTQ!HQCrthh>(MJZ)=$*QP_+FidUJYP5!xSY1avER)Qj|tADtVapLEdmB-C_N%k zJhXL;G>5Wo@w!mUFjjobD;T2oJWrf`8v*N47-JwB!CI0*T5O05GwYuw6=_iHi`C2K z851`5L#QPou@z)mYXmvPxU4MFg@#b0*kY4~VMjL8i7wLTV+wi>Q>w?0r|74WV7sYS z8S~RW+9_&H5{lVPkEI1RFW5UKTB3+Pb9yx{h@v|XZsEUHn*(dDwa^?i&~ZnpWY5}p zuJWCd=7ZKn#r2__EpJzz%43{tY#DN*d;1HyX-tkdj7d|m6LF`pUVjr4YOZi_cTmpI zLF=M(ve7wvX4erC9~>ru^O9ED%+)89lh|D#%}&24a>dCpdv43?GSi}wYfB1u*&G#k z)~kyB_4XS;3q7~`c~vNPABJ7sRjfGAoPp`2%>DY>=^^!5G}xbvm?al!Dg*b@9Bt8S zBt?nRrOa>K0_ulTWkM(KLoTpG`EpeWG z)vlnZQHqlNy=K2PrKuLlWeXQZ3do>a_GG%y<8pG)DBcd%dy6zAvC`2i{Xky}$U6f(DyRz9RhGMr( zMdZd#BNJxmD6UKRftYjQS@J>!ybw=zEiH0EqAwxcF*I!A(o|g6+Xf%oBGb=sf{%oj zA`B0_`Xw$195(bF)@fRg?XFr?3}1@P^g(6I)poQL%Gs)h&^BjN`Vo6B8QoljIK?P& zFE5c+htHYRS3sW>KyMDl9RjAbD@Oq$sf4%Ng;av=R7l1WMCcwOx{r7UT5=k1Qb!gB zHE$)dnQY z`t8&cz|$3%Px~eJGYO6?Bd%Ups!sBbIrX<$tTe-E-`$X6bX@W-t_TCflR}zL_)GcV zs~-xog2srHG=^P(I6(^YmWoNeSBkKL_r=oM%ku5fY0fU~h7q({3evPM;AKOiY=E(p z?gZ+;gtv?5&P<8RCX<~;M&0_GmGFiQp0@>9S2{k!Nr@y}x-HH*wsqQ8tw$W%(OC8j z6(87FeiBK^5()SCp{dQy^u1^H5S*i}>j3*q9KHRt`GRa+e<)HHlhumRqj+m)ZQ4mZu}kbs+<>x`w8PGM5U2lcrVx&`E!?#;8#1Joxf6v$yx9a94sj?ZsimNLWG37f zkRbJC|Avj0gEgwr5>kc4?w_SNPWL}4%di#ENgAYx#s%3wGBRQLQUMa$ZNWsLZleE<%cn^enkA549msZ9I|12-scYAa;4PLt()}6ftAS3ro>C4q0 z=bcLmA8Cb;ktt}Q$m=C5+T+iTVu8m)2~xloja!L0U<&FIdXbF2g$(F=%;5WYcj-gf z907GsCdO#-jz=wkL=S8h-bSoM9$-}VhTVGqv(Oh~R=&)4w{+X`Aw2p_a{$ta4Zw5L zvHiW*l8&W7qxL+w?TxBla!*#IglkyQ4mih+J{spJP8Tflxo}QwstcjF2*W+p!1Ncm zM_Y3n4~El1?f1B^hJptj-1}Luckk|D-x>;aMp;rXw!wdd;N|LVSCPLk6kfkRR$jLL z&I{!YGx-b_A53*D-bYae8WJfU-*-<7$G@>1*L+-Qitw>Yiel;~_Nq$TOjs^79C)c2 z!XnSZ?&Z8QTR5skkWK%L1tNXtq?!v841w~Qi_6u~9=i5UYK%JW1*JzrnzY>7*5@^Z z-m@&P9>gUYu7F02_SK#giL63)owhIp9D$7Lymugv9Vn!0>65NaYv#o83ygT=!O+{4 zNSv9S0pVPeWUj5c5<-*tt7`45#b8%_fdECY$Dr0`l+E*1pxp_(hL05nl)1y zW7Fv{*l=yeloj-v)S}A3mZo5K2iZT0sgu5JDS&6dE}ikmMdP4VtsK|KG(yRR#FWJh z5kh-N4F)TSP-NsPKKtPedgz4mGSB;+Bi}QRXsINlTYN@FF3-_3jHwhGCb&+mCFPG|47$OPTZvmUZmJvW6Ej&)N@1@2Yo;|A-@BTh8BuD!ButLzb6OvFNcv>Qt za%`QDtF#mO8alZ%TqW@b3@sTPMu{aNGJK6tf4w|69(oZ{Szt=#e#Lp67=*Z{zT@<=HBQSon2QKclic?FNxPz6RYiEckij<` z0r2@>3dZy%mYB-m!a`vd=L4fnRsjUnlJF!f@7i-Rvljz7xPyLXKEXQQbXX4dJKyED zvCe@I>$2?%*)Tm5y7H4D0MEMYyL~8;?f`JG5O(^`ELE|&-Yv$(7RuE*i=NHv8;8Yg zBt&4RD&_N_^F9T=5k&4TQPdl_WZ;Y3x8X5q?UIA}3KlTvoq*yJc5&$L`};&8=Q{k# zALOGDvU=O>^Jj5{bl;HLWx z!2_+~gOe*s;zt?eqq4&S<)t{y>#&o03Qr(--&sY~vI4)E_wJJrz+G)$sc194IpUHa0EwH09UD?>Q=D}T49U6R8 z`q}(J0(v`1=`Ja7EFo{AzYTM5j)u}B^b?v!+;*mp3T|yhsFH_%r!r>?`87^W1+?~e zctVm|ac);pn|0V9B>mZP?$7H&BYBb*08Q1!k`xnj56YTbS%4p0hgjHQjY_*HC3WN(jwZrD2S!cofG}zejk170WSKIC~Ak+N*8JPs_(Gr z9N@$xC_F7+M-z(TF6QWTav^-2%vS(3fG^cmBZ7`Jd{h|TwW>_k{GLqMywN#8R2p8M zh~numfZef{)pB!*DHlSOeaH;1vD`V}1bs=MpZ$ad;jyn0zAE*s3!eQ@joa(aQ?*AJ zzn7%go76sjGlK@uKyukl|nKb zVdLmI2kjU|WO~}3--+&yWGNJ+4Wk~YkQvIxG_06L3I;4o*3O_C$fAdQVPwv>j9}pW zimSL84L}Zu;ll$8HWs9rJ_{wg)IW? z_dtDgp`Svr*(7q=uUD-tvG4$t!m#?YGBY!m4dTTRx`CJCh-|xiF=+@Nw@)8C8n>Hu zXejE}D-kXUrS+6p`tsV1hpG$9wjOu(t^FlFfB$jGHvpPAF^L)Xsy=D7~!#j3j%8==F;rT`R;*n$@668VMD!7HS zBCfe3u0uu8ys#o5C&X*iCUo0U9r|lu+Eg#xEp3Eg(gz7}@W;{{?^(JoVhEpw95d`)M z$lv3O(dA{4teLJz47d(UF7x|PH9Oqe+)X8L@>Q>cn-&8K@!p$o3mue7S0jiOEn^BS^|+5zL$0lt0GKDy3vjcw=_9p#?E8ytb!^eVre)sPM1BDAhZqtbOO+5-AO zewFD4(Gv~&y^uw)$yOL}HElY|ptihk!??^pVvS9YCd;Ao$?ILDM9q}s$HZ)x=k`9` za?<+O(i}NjDI^Uzbdsy;pxF=Av(|jcljA6ZoUfK^Fk6Dr4ShaKXmd_aPZS5nI2mR+ zU&AnVLgF;obB0NCSWLr0?xgRky%s}qtr)jov=3q!abZ{(w{n&S)NG|1$gzEL`vM$^ zbSPDR)F`B~&)wuY?@JcbAEWI19MUiY1pK8)*tjenK&`~$!FI2q`Kg!vSR>`^|N zFb`+{e?qqb0+qw2|=PFQqVY z8<@F{K`AI9LbWgZrOTF`+!#mX`)vC6FP`|-`RMsfi>MA6It>)0Ez2?FKg#*uQW=QN z>=jm|9F*C?vO7Re1qqyl`A1XFO1Bm30=2pRb7hzz+G`5QkDECoMOLkEAQ)?^j!JA-?3P zrguxbSg#C^3^S*E)fYPp(pVB>ERkT-7+_LQj12Z)=e0~*uvsJe9<@FB0a+@>R6QN~ z;?7xSxGXLU`k?5e3b}s;zqxL`ixZ)T8~#4p@IV#iJn2Wj4=Sb-h=!z_3r9jCj02wt zK1s&G3?}-Dh#IorC8P8N?D%^_4@50& z9;SQrq0H4+G7>8_X=8pNBgvQ5>(|QlKa-KLwVkD|!@tN#?xPIG3mGMYD5E_3%tZSw zs0;aG%=Q!;MbaQ3#6@!INJ?f0Id;$uTh&q)@sxX#n@&FaJeTyfP?BB5j&^e~*pFqn zpJeiQdq0EP_{y*&1Wd-eplD0MgMTt+)-V!*O3-J1%%S1O7V8ogq+KRmKHwB1AuOOt zEbgN4Q!ZVBkDGn=XUzE!F92!U)m+Kx%5datU_-%PMTlgv3V)IbuIP}ZN!%^Mgm*Om zvqq-|z`BC5@%nI5&T}1snc}*P>v|dnX#2VXc?+WSVt?Lswgc1Rz=cD<{RFnj3Pxj& zcyFu=dq~-tkJ3mDbS1YSpEVy2bx0T0oip%c?!&sh;lwbiyO+fbLO&;;RTS^_fUv!( z3~s6-wyp4omtHB&jQm>s+sLPR|e6xO`##Mb9?n=yQPm=0t( zsz~QYUo1!N)QFyxvZszA=M_R78_8%ICIhw|J-|(JeNc1NY|i?WWStRI=D7BpW%GUl zqy>T_k90)wc7!2)+*1;M+-Hcw4b62t6!IvWM#+H`#4=ihGJWM^kk^-diT+^t3Q z3uipP{QSDP-$yjj-?7ST0!H3Y*Wj0Pcc}d0+&#kpBLtD%`XS}!Q4!@5)r%r4>;>{i z21zsp0Lo_&f}y2Jb{Y8j+~aZjdf#_92Q(hNbSJT14{uz0Ft7vY+pb`0_Vq{xt6 z0dO|>(}usTUiN5FL6BSBfr&7t91Pf|mjk$idpzhk(L3nlJ)?jXU^`wBM8@al-0o-} zVQt)lJ!5WM4!m&;1}Zk~0m6c>dW4g4Ktsfl4olTi&(5%(qgkn_Xc12`I0OJ5_;@F# zjgcH)P#;Cj20vnVgUf?%@pX#s(ir9>+=|+CN5lKxP`FxmaRoe~eI0Pma5fO`FK(1U zW}w=jluMhpnHoaIqNj|M6chm0&&U`fW$Ml%?_@PT%#8wN1U&;O^mMol@c$t%0vrtt zI4Qn8D_PZQygzVWY8M&OqKKJF?mUm^s|DL)c+TMaY&Jq-ME$5o(5#xj*KhBeMGv(5 zRk%bKV;Lh$PqJIzOA>!*-Ajr}6MB|Pk#=-Hq{btP)LdT|n%X^p!-|jh)0bl1tjK`s zo^eBR`NXvCFc)NVFIiM)WijKZIQrcbHHm(TllB=bq6jROIBR7K0x>RPgF+irct)Dp zgpG_?`-d7=752}84sLSabY^YF%^@A--K-X-WYr#2yPZYvbiuD&Z)_ZN6YZ)9LXmn6 zB=BQ{YXMQwF31~_)wVGW^;QN%>&3C#`wJ`$Cu?JxLx++!AhrV9CSv`$`l=2Mc6cWM z1Wm5*IYjx5^AqFj+^Rv+7Z`!|+=k6f`z3ZERCMUdpJ6RRoUGkDNl;gL_4e+VBuS}~_ZWJN_>k;g#x{7GMLob^oFAUy9 z&r5jMFi^gT)DP#^=#Q-KYkrm32@YeEOwOc^xSET(3b3t%joH&!ig0+C?Xulp77f$Z zn}lhVvSLx|yZoatu2Ob<2N#lQU)*s%Z(xTNOEtIE|&m+(}3CUgp zEGei%_sj?0SrqeN00C`-^r&gJIFiAA!-4kbxT85+YFAndc}`A3*u-8B?#vFXC4IWg zDb9y(74R{SZf*2g&8oSJsKTX?BRlp~>7Hig1vehVFi)#6K&Bm@nB5W@*vHwtuFaXk z3d{!`ri<%!U*AhoJAI*E9D3?)+u4?q`G(cuhbO}CgL#!Rz8JJjMtN{40(=3jn_(>V ziWuBYLz#t#67kCT{uLX2KKZyNojk!9UnKDy4f1;R+9@hiO~9M{(iSH@oBi!YNZ31r zL6RkQ{P9pR7|O3OW!XjHwd$%-h+4rf`ohsdTpq8G7{GH?n4oQr{lEGige>siBg6rG zN-ZY|vmk@IB;z^R6Xu6ZA6YtNuKx68&L+HC4nKI!vO6@_J5&%L6iDZ;VLef71NOXv z3rVPy*d>%^+=<53g#l$&QZ$<=tllS)Z;#Kw%1f?PY()8V*JtVB8N79$S8;_sF0>)? zYXAG#l6CZ4(-llEw7LG0y+r(HHHg94StxUg1mL3~USa6OgjnXYmr@&n&z@^e70LYw zxh7@g-~z~N!pu{K>q`0HaE=rO6p{54ll@gaK@OY#fg@$a&X9Bahxh>!ZGAOX>NOSP zS_=2>fO|zNdCqS8g7h#ClSJsfRHt6ey@fPXHvUg^t3CE1RRV)6X$jH#n(BSu_K2q6S0nr%AOKYy$Rt==dlvkrHPYg-#%Oqq; zB`PB`9GnmzeDWEnvYmZAp8BzrhI@NM(#k_9jhCzlV)4k8x$z^>6STDTuIFb6hyn)s z&qpV&V{9*C5+)-Nuh}`Y;mN(1UL1PxPJrZ z>5LS!DwRbu3#_HVSM*0jbtUbUrPHOWfE@s0(4rkdcVr-iAYl3?i>gYGEv{x{5jIL5 zUrMx@m!E+sUvfaFZ)ChiY{|2hOMBjw`bSe}S)I>F$ zSU2|*AV`irGUNGFg+QP}elavVGo*JIxKA#zPM?c6J9}qke4k)Ic)-R70(uM8BA13h zakB+UdKW+$F|!85bad60AZ81NnlAzCnse++XK;v*@(Uq^$#k(4{0+v6y*D016k$1e znVpF|m;DI)NKHG!jLxl)+?=2EjLB{i4Z=hmkGibLYtRXnc+p{(F`RUw@oM8u9FmpJ zq0G-Vh;r+abG^g{i#ix6Q%+~rc@q&IN&>JIR6ss}IkZNopk;6}r*c|`QIju%COm!! zioC$fq$g-AH|lbBnBtK|3JfEFI=mRjdnXi(0h2G*&381h)?ZE@rH~=PDzd}R8B4!O zhYc`+zy#4nq(H!EKS|H8=+yWAF+pAS{imcbPvI`B*aNbx11S-FAUjHor7S!D-1te| z?|BP-l@u7#es1?j0uuA}7hp1SGwjlvP~Ga`@rc~R?{}ayhHI2zl>3C2t}x}OWycO! zNit@3GIwB0Nd=b*&>-o=new_&9U^pPTc!8Fb>rM-L+^mMU}fm*7HwQ8}04fglpZdH%HB$ zgKl?Fm9rN4YHAOb+HO!`sE3gF8TdTI__l10ToNq zta|=hDYDJPhkAi`n`{d#iB}Lvs8rY?Jn9b@(zt*jr$C6fVE655D}A{+ekW8blrt`< zHBJRp+~Iw}s=Le2aMW@2OZl6SoB}`72Z#nFJEhIGn69+^jixr#T#1)+68U{R5UT5R zdK>ycAL2U#XBR@9#C+o4mMy`T{8{h+?=MY+X1Q&rce7@c?iPBYxi??r<{_N;tI$jL0? zwXQS#juY1fJ6)13kC6JJY3r6&6bw!&=b|68pZ4q?ZfbAzK^`!BVh83Z+zG-j^B8aj z$|?NMp(4J!>IJd_d8yo}FeGw~Syx%!#MPy5aS<*COR>uK6weBK6$QxLL$ertW`q6? zW1@b8B%2++w!z7`9n}%7yb`Nsk{~rATt(oWLz5FlV6Y@1>#k zuUqMO(~@xIw_J6n(Fsk4((xp|o7jqM1enPHC=|Y#g_!B#57CP;rXz71h6rH|DesIi z2fw6joj2cm`L1p>#&l{L7T5Ynx4?XaN#2C5Xo@y>v^F}s0c57xm=S#-XwB~JrxpRt zHaK$y)xKU;RL!qlpY_Jwg=*LGnJH{U5-)tJ(!>@(qu!b%v&YSCu)ksy4#+H)Lo0@& zmfb|ZdLRuDAcTeE5dhzW@|^E0k;h-vPe@dy`Qb?$WlcvjkJ0aZrTF2imyDE(SjE*E zSRC1wC#Dx+s=o&tp+J~H3A&PZzLpT?!~4Vi8{jGuo4ZiPiI)u+j$sRqX4Qv;H!7x; zWv14>;ux0*LT0g?Vqe7NqCN{TDC||6+SKIC1r-%I0V)goWg-TM^i>z>567709h=~UCc`DXEy?QL~ftHIz zO>5o%^#k|fPP$7RK+g48^4$U$zn}*SxInG6Gt)D+s^-l|w$8L$2 z$Zk3ZqrBs}ejcuF@2W+#;+9Fz$U1kaqanMI#e%i%&TFzL3TDA;l_6hsRqktN5vA<{ zAFUbe&(k*~B|Ws;ux&7u7U&awbMhS-sI&pspu!B%fH%$RN*P2+3ATH*hRZbUr4d4p zHc&GRQCSBXF=S}JqLk{e%C6{!PA>0M{ z%`n26G@WJDm2N|2q;0(~y~ccTI2 zqo%ESJ}xCz5T_}KM5{_{QfyR@Er3^*I4j?#A&HI6IuTOu3!4F9>Yk_-_}kyXab&3A zY>0M!)gC{iAjKFR>xaoILFR6bj`iAU=s*zX&TNE|(I*_3ECz}UZ8}Dxu$p8Lv+l{U zo+%U7A{3U7*x(*@UT6784Yljem9b?xH&yjHC+g&`NIlxJgMT}@@501<~a z?L8d%XmP#Vx;VKIdWH%IEA#$is!B1_hq=hZP!gf93J`waM#=k>cO`M1d! zc|&6%7n^^@KVu|hx_FTS$IdE1fX_oihDKG6x(fRW7|}$-(8O)Jjnqk43sAnJba(M= z5A_7{bl@7cQV~0C>y3^$?<|dcy}W)phgbvr#wibVP3FD+KJsBH-DJ7V=#rOxfmd=N ze3z3RDv>m&X~JK^%h`Z6Z5*23DP9CU0uq4X1w#% z%7$nrBWk+oQ7lR_0AsZbM>k&v`{N%8GGliI6`_~#>O(slwS6VmH+JhCi7MC~D>JVrHI|etq_D;lD1w;d+FeOFL zlF0m#ogXW+LtFx!%;Dzt+7Xy+`a=QzTtD3$P=S1Qj5tX@UI>RdFS)sE(;vY}|It0y z$QR&9f4a5bha3Li-I|J>u8qyVBa@zrmNJNZa9&9_jW=$s0&@IokZpBJtSi8|$j}&h z1>iLW0eAcA$#z>O>bF(E0(Wv4LV9_tj~g`SDrm;h@u);vwVyw8FpP0>Ff`$|w|M~3 zhWC>%patDl0v(Rm1#C}$L&Y?jzBo>6qN*^8e~)Y%WVNr$KG`O>y!*xH!#-FNc)dk3 zLp4*-P2TZ*C%ph=dO2-RHuH)NG+lv1i6~*`c&-2+hk(l9&E0U4#^Qy^?_0x-s5phHth~Yu?)5Rk+T7b4!08Ikwq-u!iVC|kOPuhi# zfR{zlA7o8qXTFKZ-`DSEiZY6V&Yopl>idZ7jK)rrl7yl4Q3Qnz^t!`N-9tY)m!>6& z;^Q@$NnGM^t!rllvv7C;__Kfb;{QOKF!dN>yG=U^(b~f z-Aq&ZD#Fn9?7sOsj;sSrK!@+TC?=UBY?C7x@ZA%ajBix1PaHMY-ISClBBhnV?Xo$i zL8&*rI3I0h!=~f;K{DcJE#!lmz1S(+yh^{%rBzKMv2K9b#6t(nwGl;S%=5ZFf9USe&!r$EdRhdhex zwc-k##7}rZwNRjwp(aSGf9TzzxVSm?y9Xfns(-8vzZhDWsIWwS+QYD3RG`$9jqmM8 zcER5H2vSHVPKGDeL8n!fKN{6!>-;^%Q zBY&8ejdo2{-O^TtyzAk@ zK|@H33`ze~|)DKh{H(j;u1(i=bbnm|vg25%(OAxfDIKJDxKre^&yrJ z=a8IdIsd3?{3%TbbH1S>lkb5vL~Mlj{0GTsA7EP$y$}rV%X9z5pTCpL?_F$PEKYgn zxBL>QaA<|VhsFhz5NHS=FsXhFQdh6JP_^JzsBDRt4>40@N&Cdskx)h5V&Hs!0ofM$ z?9)Tc?E`E0XoUI&AR_%pv(0)ZCW!P2_2hW8J3p2xK_ zD-~8qe5gx;QBdc|qGw{(PPk=e%I-y7a8wD$pniLBq5v55WoK0wtCVTVc--%L&f}V{{cxUldh3ZPWyW5r{RT zIaX>WJ|z2@1{H<^aI#3+lDbH+lGmohnkm?7b}E6`L+!sO%+{1~VASm91R=u;V?#Be zwRrLYb|DlvlUC~7piM-bKE8p97u6d;Jn+?L2|SEr^A>+}MSh`lXQ+le@Gbw0|MhW~;i-$~$!0io_Yhi`sWb*I_44k_{v$Xn$Vo0cHMn%t-F@%Je zGiV9?fj;!`*yfFdl=nZ0j6xoAXaeR9WuuR>z9RuxTsTKtipGfg1_XR)2d^v%HdQn! z{yt)oMGg!1z#k5$h>_nW$LA+Z#mZ|RXA~qM3QglpObgMhbVA@6btGo!vsC^mg{!|y z%*fQIOC5LiL3Q6>x+2-sIr!PV&h;+Ye4;-LQP~fTu>wv2#567`w;|)2X$vST7qn@i zPhUQ^hurB5_=DS|BJ)7;Q^J%)=9ov&@*k0IXEE6VxtDQ?@+UsH{jE~!W$)DQ-q3%> z<|`F=#*B&g(8BeNY3y-8xtqaX0QQ#e^(HQZjD2@d-z*YKG58sBcS2l)ft)~_oO=4; z$zV$g%L>oA2|)vh{*Ca2`Qyjw#I(%qGjn{mEAcSQLQ@G=IOm`5anUn}%dVDcX2k(- zxM_#p$>-lzb9smOQ2UOo>C%0jnKrY{sKC&v{YG6)k4PGiV)B7ZnHw+oYuaXxY0H$D z5_3T8ItEl`?p0b~%ytVZt_eyHrDhM9W!g26hS*ARH65fbv`Q=iG320gGy{_*gnR7` zTrNSfGLnSlld%^1R>FgL^8|uwgvXoVY4+QOXRmbgOT5tjLD%a)d8_KRL)zQM1g{;^e#z2HS>0P$|Gi(@@4DY^*?SEd z{}OM8{}0{2!+w35Z;Sk2(>K4Qlk?B2|Ksrgwgd1_U|%)|yscg)_@(aU>OkSL1aE8HUdtYTNxkI%qV%tTzAYDeH5h)$=F91Bt zp5@zGx!)nat$e}zrR(MDNbaB2{S~;kB^7@LXQS{p+TPaA{0{DI5yPLs87lp=y1xSV zHof>-`0z{kmH$rLZ)Ff~0lwxX-=<*xJg{%n|GUb+0{J#h_UBnhr1du%|1Od97VK-F z?X8pl8fg0^lLr60>c0a0)>ZxOwVV8n=C>h(w?JS0+qdo^-Y=9dSHNcftnRPCy|ukx zUCLjw`ghv?TIiLHBQce+BQYDe&jjLD2gDQTLza|BsWD-y(ap m9o{aBueQT4S+o6nRsR " - echo "When a read has a k-mer in , both pairs will be excluded." - echo -e "\t : meryl db to lookup" - echo -e "\t: read 1 fastq file." - echo -e "\t: read 2 fastq file." - echo -e "\t : .R1.fastq.gz and .R2.fastq.gz will be generated." - exit -1 -fi - -meryl=$1 -read1=$2 -read2=$3 -out=$4 - -echo " -meryl-lookup -memory 2 -exclude -mers $meryl -sequence $read1 -sequence2 $read2 -r2 $out.R2.fastq.gz | pigz -c > $out.R1.fastq.gz" -meryl-lookup -memory 2 -exclude -mers $meryl -sequence $read1 -sequence2 $read2 -r2 $out.R2.fastq.gz | pigz -c > $out.R1.fastq.gz - diff --git a/scripts/trio/fastaGetGaps.jar b/scripts/trio/fastaGetGaps.jar deleted file mode 100755 index 7687cff411a277aef2ea4cffb82204af529c69b7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 18792 zcmb7s19W9u)^4nd?WAJcwry36iftPe+qNop#kNtgZB~*>UebTRcXRLGx4Yjw{DghLvFPk*Jj1+VYj?I=XFm+8`18WTN=OJ2L;>{^UZA(F5ai4ilxg-`#aS6l5 zl*76n+0z?V<6$@>j+sM;!e5`>{=|FjpSclp`nmN_uW|3M*tfg4+m9!_KCG`lqjsJ? z@QJ1bQ4i!XPBq&wAxX;wJc0EXOqHX2y4QfL<#p>GJ-GNWc#2%n-#RLx&d!~A zW!mY-?rrRy9!zGY(h77owlWw{v2lLG-6qu6rY;}Yj&c^vVjmsR3@`)GY-}2FZlJpW^^o`H zL6VkA^WVm7Kpm)$TUzzvjt5}L>x$B3`l_eih;_h6$}c=)%80R0Re}6+ECkE9YnA4E z3c>ihIlf|6ZdrxFfx6vrgCzmjetdYZ24&Hg)8RzRv-6k5csfJXCuasx;~~X3ae;MEu^$6c;wUZ*@O}s7LShc=3p(|l0XIh-s{$bjBe$eD$vSGTp-`I>LeCu1|nhj(Uhs4=8p>^XFE z12qdkfIc~;gIq6WsS6*h?hB-pAglubw-pkZgW{f$oGH$t7;cHbgJ+afcZ*!Kr0hI$ zW&noU9TF6V-q)i*JZv9zfXi|LfjcLcA}Op;Yg|a0S}-r@taRR~3d0pd1}dd4v*UT% z!=aCPBvtO`JtuGk=k5;Z${iZg{@_RId$_!B8vV z%?=4Qiviiyuj{bu0zb=<1_`6R4D~prS`pdQ+N;7Gvs;W+eKkQi_X({EdpN;WFvp;m zuRpsM8Y&@KC?F%Y%JWzm6(xdM-w&p z7X}sPYh8NLX^K*x^#mIWlnPNn{6sOMQ+nknqS>UxE)-Cyjc zO3OS1sx<#>j-jJ2^o=@2$Hy80rcX0rD?^s|q$!M^AUw)Kv-<(DoKlRcn>IQ&$wdW@ zr#p+CayVw!R8)TAJR*Lcp7OR>0GK5Qfi*W&&t;R zn#s*Ym`j`r@A?{fedLl^V-4)LBG}#0q(i{0PQ^H21hvR+hp=j}ohs=>yeR!sc+Uy% zU~6{MUCQXsAx*r-R_-*~jrhg5&0qRK!gsIacTjz6Kz-V`HL3d6reT&d`eA8YX-PJA zfjd5)u)%vBTh``Z+x@3G#jn5QyRwF!l7zb3(0ao|#}q@~#RRpoyPs`>EJv}54FG-H z_=4SQbh8BoyLmtR8}RoHw@=44&kHH8921^?Sc-1ao;l6f0(Pp=obP@}5e6PbH+Q%p z(pdrRcY>8Xi1p9;nL!i8%9qp&5HTLhA!I!tScN};H5?5uicgx9oxHX zD>uWB?Px7~hl`HvE53^+XNrFF_@$-8!#vZwa16oO-hG5~A%W3x-f~60X)qk2hs9>a zen%G`<4A zIc2&L?UX!^Oy@}?p&KOxP-)m$4bt`oMDM{+@A0wWI2jepZ?y_y{r4vqa}YJFVLe$p zfU@#0RKDB;vEDgk2oY8Um>EJAN_<{&V!i$xs22FVRG|6nQFs+dgQj3Ep;t*5J19V| zr;NVO57$0aE#c6Y68A-o36 zi-0mnjUbD2F#~~2MQka@ShgI0RSH1J_ZVP#+lL|!WIC+{~Rd1#}TJ0 ztJX7o(B81qFW_;()Ta{tl%-%H5fX9zk96rxm7%Un-?3XJ2vFRomyoa$w*g zO3V$Zp;-jk49r^~F?3C=UE#}3Xex(jTO5=xyP0GKEf?)%74 z)yPi$vxdb&Ji*e)t#B}Tr?k#5wagbTl?$a#G}Jv4S@ZYJApE%bbq;Dvy4bd^2>rBH z^=$gYTpBD6e7gx%IfE9Bm3xU`uT=dGd(}w6KX!(%n z(&%9#7!T>8U`0{N^c2Aeg(vZVe@Od0rL6q@=(**lI8vK!kiyDcZnZaRsvfATvA=@q}BE{f2 z$sDSQO(O1+E|goCq^@tO3BO?J$l)=IEs;-J;hFkQGRRdst>R)WBmSP|u@^|(7xStt#P5R& zzT*r)$g9tv(4SgiE`<*ZgjQY9`D5Rca%jN$ zA*Y>f5tKxaeP7sy`Gv?;fE*EI!R65HbFoY(po4{o(+rDr`SxayI5&GJch>?&7N2h{ zHjA+^p`Dtv&y((_WQ--yD@tR&x$OX!GUrXSwpyel4U7AGtjFaSU@%tyuhYlKevZxzql z(A3!U{|4!SN?%m5l#y3i@t~oB2U$S{{DqV*oP9f|v0YD#FZ4@QPvNv5Z+<$dooK#3*5Q2x(uH*g69C$F zzaxC2GkS7zB~ADui*izWoUgJHt92W8)J>qA)2P@a$m&9;QCx1}imtNQV*+M?UB}xe- zJqNX~f~0XIg)GZ_JgCrhmEJtU_P7C7UxivK*vR74=}a*{Ip6 zN+9LNSEZ)WvHhz%G+U(~7R~+(BYhsnmG7xXtC=b$IY)Y| zD!0Ep@O?~RI_6&S9lDx%EMe=KQo!bj#E=sNL-e3FQzX=>Z(!SeCL*fxoS6qtbyaBa zNlC5wlO)V;qVhvx;6!}x)IdAd!6Ge{N9cES&Dh-xT~)lA@=#?DgDw@WR*GBP>T(#J z83ZEI8VMd(F`G@630FnP_&iThFiMjwK*n>Nhe#zb-!o&+I5W z+1I!^!Qe^vPVG1|MO&yvQF_u7T!ITXkS&celda75;V+q823l@GQ-tB)UjnL;q%6j` z26pI;-7Dv#OZ6&)eofjNdd(IajJArkEel~A@MK5*dOXC~d4i9+CXU>rv(iHzx#>S{ zz63n`8uTqSPge_?@*(=v1I{p31P3_ISP`H5D-bGEF1isoj@5a*9z4~0 zMF{#xi@eF~qZlUo%@&I{6?aoBHm^yGI`DEFBQ!j}kTL~hnFNrOe2UihuFUa2#;;UN zVkc@EzvQHwpo~aMJ@7lz+mkARhO%MO2NO0!-I{|F*G$HQW6j(fasyrVP$-DV-jx*! zTv~ILQ16bv4$LewBXbg~(6$OI^NXy>uf)4SCU2*X@SqL+#_|Ahazu7-{ccjQPL5(; z-b@P0n?L^}>3>L}l=$e&N&YL7^6t((6{oDW=-@qM@}GOcl0M0$A2V{7h*2vPvupuuu5 zwpU&5Y0y~>^V28N$qa`p=QIgOcob_*8~F~Y$_5FIhooKvZb6kQ56RIf%!+AHiwqm> zRzp~v1DeojK)aI^LMr8v7RHAu!mvp%Y{n~TrExyXFgpD ztGpYzN)FBko^X)6p)X?@>l=!`R{~9 zTvpk^)Z%~2L!Jt!GD!T$+%t(*no#l6iBUjM1P$M6_)vJmihZM*XfE<3{OG4lIXm@R z`%V)(2@kLyKwkQIiDP=fFH*uU$Ile#sn4RCozK$UCzcqmukY{i0PI|HItQY8b_Sk> zwozH{>$@9NwqWw4hRm0XDIn;2iTKKfw@;^Y$d;kIPT@M9<9@j_$4Uo8^7oA;?4TA`>juVzsnN&ZqE9&^z_f z6uzt`AFUZQSs&y-Bs=4oA}K`gjxa9S#HPz<2ozAA zYZN=vWH<;}2A^(&g;3X_rwVG%9UH-8`4wYgdNN%GlSk3uA}wa7tS}*NySjMr{ho`? zzlQe2(MmCK(4mV$O&8sMxPh&@K3ATz6l$qTzR_$2RxkABBEH=@EiFL;1oLc!^>PE# z#0iNe(yXICL*|mJqe%U^VanyzJ=cJXh43JhE^fSPMT35&52x!UL?h?5o zc?&esy6-aEy!bruJE|tiU(PrtCrRy z<6^xwGCIPN{KG)}B1m&ZoT*rnU2~9GBOxN#f0NHL^{34S@l53I^k)?5Xj6?en5ze8 z*^$!NOqip>FRB#&+%?a=O)FmOOS5%VoLJk`Dg95y3X^Qfj|W6ImV^y2FP; z(x%(F{xQl{PT0h|Nmw)QHzMx!WPw;=X)fu?@ zF@`hN-`7CAN)(=7|DN;#e6eHNuVmI&pz`}em^yY zeK8k$i;N>^5fNh}xOAnYvVt5t=|`+;s0w+@JSj}4UwmGQ`&+5VuVTh~xEbvy(%sK8 zc)h(}z-)YFIS>P;<6Ka6q!A##o3Lmai$cd6usmne3gC!$iwM!JlC2(biIWoL(C%XKK1Vo|gCAL}a17E#T+&g68gJF^@RF#cvhKcRMIzZz_wI zqQv7#XA@w6VMAvsBR6pbaeBYV|+fDrj7 z;fY5Yl0*mMumRpVsQ}&!WWko!CO#@fq)n65U@}rE9b&10$|>mkq9>tBSroo4`t6(Y z{juo(wc`DA(ccx1QcP-GMoL!sV}IoOW4;SYuROG7Klhl~_nE+Ruy&9`Nm9z#-HY?=7OcgvX@Ywzz*pI(l$cwa#1 zg0h~Ak(!#Vu`xyUdRWR=Z5%$Q(_I-eaG4kG#xf~P-}R7A_DgIsl~EDh8BHS;OD+Os z&=Zrf8FRFyGf_$YL;%zox48r1&{v+VF-V)I89KU#-Q;N^a8p@alHceTEi~d#1X>9G zzU6PLpEX{XALLed^i>2)9u|Dt%K^f{Jq~P=7y#zXB(~6((dwr2!W<)?$@1zlh zLPF|pm`Xy5s=mS;ss!aX2<>thI^-p*S~*VSx=A0u${*{8oC5 z{ZqB8Do1UggPZ)A?t;yvIh3P^<69@{k1`{1__*wfqYA&>6+-4(BZ@_$en=pshECL9F42N%^pPM(Lqk) z1Sdnso%pag>XbbE`|xM2`NZ;ARdap+A#5DbNmw(IOOFKFl<~b)d7}zKoB^$HFh?0{ z%{^jEP)l~rs<>hKQYFVWQFblp>5$nOyWO-s5RWM*dF!978_Yy5;8WL_zb$)$Io$38 z!@ffP5PfnJH9=UIm_)*47o^R&`$1@n19q?^9iNS|+k*t<;`3s?1;&$0I)V)I8-~Kb zbC_z3>pRoJ;<{mCJ*;3yPUH5r{R)RL8U{?oYgnr=7hBK6#G+F7#00%<*ptKqZkz*e z;S>4l%yGF~uM|S9a2wfH89}@JAit{n+MF8IRskZ2KH1guRG1+xaydJfh);I1wXR&P@Y;_-2o%kDsF6l{B6B9>M1 znnjK8>aT*>J`qQ@O10ZqzY&Uy79g(;KVaw=oaLlTB<2S`P}73q{e6w`C;PHf8}6E~D)o0hW(Ji5lKN_;dSO7_jCO;F?}UQ3tX!lB{ovc*exTrwtO11E0zco=Bd@Em^44SZ~doVE_L2* z7JLb@kltGq1Pb}lX5}|409}Y5Y%{^^fdPdD&{24%*Yn=NWjDl=L7TjttD9 zXZ|za0;+j1prDRnTI3viEa}jp(O^eZ?8%}njVqmn0v8t%Ttc4*Pev#9iUEDbEZ0+y zD#V0Gj}FFyR^{SVWWh?vi5*9ZOmB1|^xaWo3?(O-4a;zs^ z=BwLHU*BsoI|Jc9Tn3t5+lAKRrN;G~L%9|6zL<2(#<}pxg8Tt(+hME?N|-#& z!x;s~l5r{o{^eW!K6!XwyLf{!>!t7=4RibSJ18qszCtwnr7llX1VV3()42@gOI|d?GLj@i8nH^4HLUwVoSJYWjuLqa@c5uq z0(IC=u(h5b0Sgo?Q$k_obx_F|Q-I%u;gpJr*UII@2IG*oc~ z8WRgVyf8mP(gm1`oqZg>#;LTXdq-p9+EXa4mz*e4(ddo2$use97#Zt*&ss>xd`5=a zlQY)|_P0;71@7zh@;QLhHW~^F^fC=*ZGZYiRrq531m2A_x;WzwXEJPXb04Yu7zA%u zgt%3S9J*Ox4K0C^KN6ZNS)Uxe9(_67ATXmg-6)15BN-$i^Ozi(8Uv1my0JyrI7M6u z@pf)rI+8;15xs%28LwcTd%+Kri_f{uOav8dt#WUNViqY5qaj3v?u|v9TnzErk#II+_c_I z?$7=TN=z9)w?||_$)$!Xa9R0z4w-G}9*u8tNIWB-_MkOKs#RcB`bAc5u;gjvCXU!h z)93dx_TWm$gjVv=q39);bGy+T!u90ZWDX$oV%_G&h>D${IWs8$+>Zi-t0Kvv5hW-x zISC~SBMfy`wtyJP_dS%d#}7g}4Q)bPLQ0d1b+~kRr^31fm?cu#)U@)C8$?@HHJZgiz-Y&48%6*^?vQ^di81ga>8i!x0s z1m`g{`4V4{R4wMM$G^(TVKoAz#x-R6s?yh^_p9^SL7q4P;uerw;V^@qCkbj$)sif; zL7+yOd^_Q(K@eb@{U>X}4J0xeHBJbx#?zGy9#F_RFcKd4V@K**e@?dF8TA^~ybD^j zQ+_4SxBg(Y{nZzEn%IVwylp5h!Pm4w;z6k{8M7Vc8*P8%*)4TflGW@40Ur;t!L2|lZI6&j(DVcjZj6(x;KuCQ5Pt_8 zT!dit8)MtL!QL40{A=4^R#>t3l2H3@*Slhw}{ie zP4_!Z*c9q=NwPda?vJ9YTUk>yJf~WWdd~WNVE1%ab7uhhgw-1}xJcd6yRM{`SF^d`jkgcoq3ttY(1a{e@VinAM-ZK6XPVp|FQ?J|hDjtKqevdT2%1K2 z8{_7QEI^P54xU#KVjJ3XslQl(VBH`-L5=p8Ctaj9J?#?afb)&g=O13O(yHR+Hyhxw z|IwrbPHLWN$ zweFL^yhao@i{TQlmynOF6}~mA+WxMt{oSion|H|d@z+19T)IC!L>0gQ09>Iz#&(oS z9Dj=KsQy2k)w@$GP?ok{W`pOJ?9k6fN-$i?Rkil%(`Wqzyjm!xKvl&mfw^R1s355z zGu{3AqX$ksaT|GYM}L}VJ#8RmK%Vx@-3-*pW|ddG3U99O`5er|yHCY#1IcF*GJ7K6Tiz zZ!wnS8xW5<`Hl`&*g$AfV})oUm}Yh-51}Rp+dW&uXBhR-3S&eWsvCuR124Bi(XGO=EAt^lpLm)GMdUEyq6}Z3ot#ijY(SsQL!3@J=7eW!UK9iZU)M1O zY9pax;^KE;IQYv0w!~^GAy6u>!bIGQ6T-~Pq!oP&QPKCu{bhAtZCH%bLK*fuWn&08 zl*84`S$72TYCw#4fa-h41F67kcrUd~k-)sR;kD{-As}yw241G^Frw;I-Bq=<9wQax zUHwFxEETM4efM%B%4%P5yM4WB$J>P7iJ4Qf)eI#u=0@fWvutKJq@+Mcy>kNmI>22? zie<<;RMbJi+f*AFVD3;g3up!X{5WO;x$78owZoa~ZfgI=JJ~c?+o^kGf8`)qkSX8q zmnUkd6DJPPW0B?>rmHz}FQwk*{yc@;bB?hjHL~ni9E=-MEBtvccvl)`hA#OOOr0$= zP8NltVPevYM!Q%bqw6=+TikU&xgcg95R@)_BAkv6u^Mx>8d1M!*_!9!QDFyhnSx5T zsm3J6ME2SOdR0oW@oyWE+SqIoBlo>oGhi&eGxdCb`v-W=bXD9fvF;x_lNXd^nBx-z zu$jduJS|Z%UVDw5h!Q**P0+FiM1#{sz!9O%r^u96(~RQQz1i0Dr6SryA`0PKJmb!r ztY2uL_uaYEcPtlYD{HeO&mIajqKvP*PzMS+E^ju~RSyY~aOqM%!BdPEHOOyDPzYnB zt8%ij96qP07SX<>S)`p`;k}2FolC?RLT|OtfInqke+?!7wi%;fWFqWh^UwNcw4`h| zA9CQtML8(QWoXFoxavuFL4Q6Ix>y*xgiVjJ1}R%U>I`a6H}CFnZy;|co>3b$iPNtB z_^9*V%IJ^l+uxUv8$e@R3edOY-kYByo>tPnuC^Or^Kty-llu8>pNj!Hfh@au%5gmt ziquf8})V=lX%1f3idpkSRCop8{4lVuridAQ=MKHgGkGLy75+@ z{kj`+B}7ivYv%WpvALsd4uLjs6YvWg0HP+Upqiis$;tYb#OBW&0@xXy5|R+) z4tMvrjv(Z7pYs_O2j~}p^A)n9B}fPGLpaU(D9l}(|52RuA3t!7cw-#-PmlIjaU=M< zM^m-av$6TVl}S$}OIaj-c(25}raQMbL3x1Zb`M~>ZvzxR(Szu1-^%smqPyKB3qKSsm(eOt;Ig?$`T#J_JvMXs}3PtYQwj%ROD{Vi2TC zE2HbpVp+3+q0e_H79;AK%n{`06jU9$dl*U7+|?f`c;3!ww}eZ#C?Qt1UXcLA+}Twmb>u0 zLxUdXNaIMDtOcf({x2xb=p3}kiJ017L{ZtnZaeKXJPeX@Xj>C0zuc04Z3>1&xyJw# zl@64XyFVyRd6wNR{IyojAJCH|D4tj_e;L|trFc1MvS7KgY~Pz;oC@lILxH<`y2 z)|M$G{o72!F;TJ5XY$9&jjwf*B?OL;gtn2k*Vx7Rf@@R8Q+^hL&IY1O&R9q38j=#6YTOeMR;uXFo&+BO+}9ugqlSFa?|LKWq7sw^ovt)*7T zRL4j%8#%u8I#0(_TK-npKyQK5%bq+6cZyWYfB?qV?BEZJz7yC z2T?s6`~Pul5h|KCSgIH=k@Z(cH4@nDcDY)tV|`c34VuXMq=nVXL{JUbCro|f+3F($ z(fi2haoNulrFq=k9V&Egw9LaUtCIjIz}1Ag*n}iJfEBt-Dvh&r4_9@SD?1B9 z&rd`G2c(=d^N&j5{J^HrkRaxBv8Rv@jiu)T=4_(S({~tm9YmC%ps?gg%AYWRgD>@H zc2b6I-R2aQGR1N7fR=IZn5r$37KQ6ii+v67?QhFeOx4G!Ah9*F#zwo1QtmaBk-HGM z${svD{J~eR_$N*($STbLOvQQ1aFr6qICjWB3JyZ+kONo`xTu9x4q5w-Qb1aOy|36 zjl~8K&9p=dS7xY4WkVcKj~16PQ3b={g}85JZKO%UXc#Y=l*tg)#n-6R(jduFHA$2- zPpxWKF&v@VT488l$x`$mrAh^1e9MS4OZUBcInX2Np|Tu8L~294)|aB3LtEFvsu^d; z${G+?__$7t{7~y~R3=Hyk}O9CYkfYIZ&8;&1^DyI);LjW48`cNa&h)<+a; zl$1~0dpeyxGTBG8A=Z)Do(0Az2IXH_7{O^YD-U%r$2bipRj+d~Z|b-AFi|fQoGl3~)kylg25U`|sUYi88>gnWluL1Y>2DOde0!#3t1#ITX+1}ONNh1? z9-)rL&!FzU5J|C5Ty0YVAt{*m@5Txy6Zj2Ee%ZgLxUZUOom2P28Q~EokVQM@BFPs% z*8>v~oaE2H@63jvLGMDe)BVBNo z?SU^>9Z!CQUJC7yG~C}gABNfxB+nvFwJxqN3dXsJ;;h=9VL*h3o82VRAx=qR_@Aj< z{^Y~yNx9=AX@L4J8Ko9{>nl_q3Fo-8)?HYr)*MkSw1XPN#BrmXqWj9p#*j68rguWa zgP2vk&3=~Xm|nV_=eVPit4~a&QtiGytD;YUAX&Gr?w-*+wyJTgT6h^OmDb!Ifu`K{ zK#ESKH+vb2XT;xyb+P%I+3i<2KcUCxZVX}Y0_tb90TtL?OOG=?pw?elJNO-7)fuY# zh)AL-{@yn0ycp*?=?t%Xj+ou@)ffZ4pTnP!)#NMCR;8iHN_=v2L5Z=F z#8dqrkSeo#ADOE0j3h7x^l2$Z$wP`411YpgDp*2^v}VP!sB=HzI_kZk9~;iAhK{mK z`I^rP|8l1YQYFKvYxg~l*MC43)Dduq=)>X&h?L2Q6BJc|kiLLZGpM*vYO3*?dWO4R zWXG0jsGyI5;X1i6XP@Q5xqpDZ%s9S&V(;=w<+1dDB4%o6_zZg7v?)2{#Iu#995rFu z(Sx=J4EY;ZW$$dWZ1xw|L;r6B)J-h5@ez;ifnV`>%=bLw@yIRRNyuiF?fG&u9wok& zvB+De?fDl5mS$Q^Wyza7yM(`3+HqcKLW|edt=qjPiHr$zn$90UzVV1l7>l>g92@i1 zq^Kq5x6Cl2`N zV=6PuLNJ``XflHRqlVHm_tRa+GYy;ObI~>^t?W5wRJr;X1%0qDd z8W@MqMdP)!wx^nXMW5fc`(l%dKA>a)JkH@nu)zLgS5WS_Rsd0RoyH&{#php}8}9A1 zHqfR=me<5Hx*e?8((ZN>%K>$AwxNlkua(GW0yNKg4C88!+No&!R-3Dx?uANyDFK?s z<5Ox?!xvWmT63B94o)m(RAu9hs8&gX;9@z4u_^Zw*Up38(NuDh*yLzm6eBR3gk$A* zRO3d^WrYBH+2zZjEXUCFl88K0jMb;XE?{8|!b42V^nSDF7hj99;l+;rP->2&i_+?L z6|a$mMsw5*CTh0Evn1dW8%qltYZgp4Fk>)Fl~9{bg#CTBcUm+v5=H8rZmb%Y25%T0 zt-Azo=YS0uG-r8u=|Nl(UFK5?I=8oK(7k{t{OK90z_-Z(5)Tkhv{LNCi0Oz3I$0#o zrn}bU3#L|1h8h6_TfYS7`tcpt0U)neXC&y7^{u=2l9c}J@{!mm`3&lVah+9tafvLq zEhGlNo28;hG}#4uq53(t;RzW>Y|?tfF0(SFBvzY#LoCsc$LpK>2@_Tu8vvSA&X!Mh zujgThTdTkW0!zA0(ii)~H#r=u%EbtWo;sgjtg+@M_^tr)C|9#-?z**CWyjZD$Er3< zfFvBO0&zvIz?&W_{mrV46C0`P1~@I!*r7SUz-r=*`@)Se1h*$Oc*Q7Cij|Xj_o{G5 z6shYN_GmeX=v1t#k_~r6=RR4$6}sS<8;d_4@BAV-MTrQ;=l#;to4e*gWXNf+(~Y2?$JpL^=k_4iYwq5Kz@cMbC3%Jj7QxUBMGQj$ z0Zj)`S0|9Xigx_-^ZX7-#M6SyuvWDC>srk>c?6FgPR)3=5r(X{dk4goDb7Py zx$Xfg3t2Z=pufNFW%6&eeA}kqTK`8?R`V++e#l$h2mMcVAK<@i)5iVZd5O5J^T#AS zP)Sp9Mh^Z(CM_1xPYY_7KtMjALCz1=8reD4wg8UB!lVps`e%)BL>6U z=llg{fLmO+h+43yw$^E7ee4jKo1o zrxv3^CcxIC$GuK^{;W3SW*aK{UYcW{oEiV>@iB(atMz3vzE262>AcCVF+!`P$Oh3PfZHAJB+DA@GXLPr3&UzF6S+ zH(R6g^gMWZmq#9#U|P-6DfU>vT$_lUPt4hxys4XAl9ot5=2_Mgk>b{~&lWD!CZv(} zZ*a=PDJzO)!8){=gC3wCr%C9Y=0cZXf-a)8SKK1p2j6(wz#KkTYz8ij1I~*>MSFo` zA>-w9VQPQ=6x!}sXAZyKCXE^Jm3{y7;v?{Tf<6(M{rclAQKSR?m@NKv77*v(DZ;yN z3zQ#{>7j$?E>d{an~fu+c;)LwnF2cX2iup z*06|ZVG|`9iK;w;He&WRT4HEjck9j3CPE789gV}inP}mdb}!M?F7DCzRZ2f2a%nCD z*{A**H4uC%$G%DpmjUVEK)eD6TH+xrS6}Y5UQ4iSX;`yFBbox1A$YUYnP+XJm>g<< z>XCIR7?!duq!*JkxWAj5jH5H8`Xu?VOm@Fy-;kQMYu!0)OuHO)M@z~a2*KC}&BH~K z{1=A1zL83nDx|b`OdP&b>W$gBfNRFs)9a{#=U1ExOis9b1lXAnnbsM9JHjR+x}u~* z&M<<4FIESfvrQ-0XayB;wjFn@!6T|9-^{P#Z^bh7fT<0C|0AsM*?p2Bd&|!Y{*j*p z0wDwZac0HG<0IbBtoSqi-3b^M`~&XostMvhBEH}7Bj7uWAGgfB@3s4b=-;BPzi07JJMTV%zKiU~trhR}jXx*= z{=WhJSB!p~Qh%>g{Xx};e}nmPi^gAB{5Wup{|DCFRSU|$hx=E|eVpX|Gjlv>e}ndM zWc#noeVm>BGjog>{~qpNG52w*?LA2TgWxg$4(%_~w;vgNABX!m`SWMN?tSyse=9W~ zhT}dWzRTgq_U!w_&>wXD_Co(9;@`~=eT07R?|$qpz4v$jplq`L9{OLs*^eEPzbyOP zjs7LmA3Iea8GUc7eC$f#|H0+k)%4pd`FC*tin)&^UxGiORnYtm+MoIwAGv#1>W}5j zKR;|H{eQswD|_#$+sCBt&o%;)@y~GoGLieou{D4P5xQZH97tZ%)hev9)^DO z%kN?6AH>4*XUMm+gZ>hfeq`_c4SsYu@1Fb*;^qH0SpUNo=vxB(wzKi=ll^v8Bk=y; F{{aPp*Q@{l diff --git a/scripts/trio/hap_blob.sh b/scripts/trio/hap_blob.sh deleted file mode 100755 index 023873b..0000000 --- a/scripts/trio/hap_blob.sh +++ /dev/null @@ -1,74 +0,0 @@ -#!/bin/bash - - -echo "Usage: ./hap_blob.sh [asm2.fasta] " - -if [[ $# -lt 4 ]]; then - exit -1 -fi - -hap1=$1 -hap2=$2 -asm1=$3 -asm2=$4 -out=$5 - -if [[ -z $hap1 || -z $hap2 ]]; then - echo "Check and ." - exit -1 -fi - -if [[ -z $asm1 ]]; then - echo "No .fasta file given. Exit." - exit -1 -fi - -if [[ -z $out ]]; then - echo "No asm2 given." - echo - asm2="" - out=$4 -fi - -if [ -z $scripts ]; then - echo "Setting \$scripts path" - scripts=`which meryl` - scripts=`dirname $scripts` - scripts="$scripts/../../scripts" - echo "Done!" -fi - -count=$out.hapmers.count -echo -e "Assembly\tContig\t${hap1/.meryl}\t${hap2/.meryl}\tSize" > $count -for asm in $asm1 $asm2 -do - name=`echo $asm | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` - echo " - Start processing $name" - - for hap in $hap1 $hap2 - do - hap_name=${hap/.meryl} - if [ -s $name.$hap_name.count ]; then - echo "Found $name.$hap_name.count" - else - echo -e "\n\n--- Count k-mers in $hap_name --- - meryl-lookup -existence -sequence $asm -mers $hap > $name.$hap_name.count" - meryl-lookup -existence -sequence $asm -mers $hap > $name.$hap_name.count - fi - done - - awk -v asm=$name '{print asm"\t"$1"\t"$NF}' $name.${hap1/.meryl/}.count > $asm.tmp - awk '{print $NF"\t"$(NF-2)}' $name.${hap2/.meryl/}.count | paste $asm.tmp - >> $count - rm $asm.tmp -done - -echo " -Plot hap-mer blob plots" - -module load R - -echo " -Rscript $scripts/plot/plot_blob.R -f $count -o $out.hapmers.blob" -Rscript $scripts/plot/plot_blob.R -f $count -o $out.hapmers.blob - diff --git a/scripts/trio/phase_block.sh b/scripts/trio/phase_block.sh deleted file mode 100755 index 4c37d4f..0000000 --- a/scripts/trio/phase_block.sh +++ /dev/null @@ -1,68 +0,0 @@ -#!/bin/bash - -echo "Usage: ./phase_block.sh " - -if [[ $# -lt 4 ]]; then - exit -1 -fi - -scaff=$1 -scaff=${scaff/.fasta/} - -hap1=$2 -hap2=$3 -out=$4 - -if [ -z $scripts ]; then - echo "Setting \$scripts path" - scripts=`which meryl` - scripts=`dirname $scripts` - scripts="$scripts/../../scripts" - scripts=`readlink -e $scripts` - echo "Done!" -fi - -module load bedtools -module load samtools -igvtools=$tools/IGVTools/igvtools - - -if [ ! -e $scaff.gaps ]; then - echo " - Get gaps" - java -jar -Xmx4g $scripts/trio/fastaGetGaps.jar $scaff.fasta $scaff.gaps -fi -awk '{print $1"\t"$2"\t"$3"\tgap"}' $scaff.gaps > $scaff.gaps.bed -cat $scaff.gaps.bed > $scaff.bed - -samtools faidx $scaff.fasta - -if [ ! -s $out.sort.bed ]; then - echo " - Generate pat and mat marker sites bed" - for hap in $hap1 $hap2 - do - hap=${hap/.meryl/} - echo " - -- $hap" - hap_short=${hap/.inherited/} - if [ ! -s $out.$hap.bed ]; then - meryl-lookup -dump -memory 4 -sequence $scaff.fasta -mers $hap.meryl | awk -v hap=$hap_short '$(NF-4)=="T" {print $1"\t"$(NF-5)"\t"($(NF-5)+21)"\t"hap}' > $out.$hap.bed - fi - cat $out.$hap.bed >> $out.bed - - if [ ! -s $out.$hap.$tdf ]; then - $igvtools count $out.$hap.bed $out.$hap.tdf $scaff.fasta.fai - fi - done - - echo " - Sort $out.bed" - bedtools sort -i $out.bed > $out.sort.bed -fi - -#$scripts/plot/plot_block.sh [include_gaps] -echo " -$scripts/plot/plot_block.sh $out.sort.bed $out 10 20000 T" -$scripts/plot/plot_block.sh $out.sort.bed $out 10 20000 T - diff --git a/scripts/trio/spectra-hap.sh b/scripts/trio/spectra-hap.sh deleted file mode 100755 index d24491b..0000000 --- a/scripts/trio/spectra-hap.sh +++ /dev/null @@ -1,107 +0,0 @@ -#!/bin/bash - -echo "Usage: ./spectra-hap-cn.sh [asm2.fasta] " -if [[ $# -lt 6 ]]; then - echo "Not enough arguements given." - exit -1 -fi - -read=$1 -read_hap1=$2 # .meryl Haplotype1 specific kmers with counts from reads -read_hap2=$3 # .meryl Haplotype2 specific kmers with counts from reads -k=$4 # kmer -asm1_fa=$5 # .fasta Haplotype1 assembly -asm2_fa=$6 # .fasta Haplotype2 assembly -name=$7 # output prefix -if [ -z $name ]; then - name=$6 - asm2_fa="" -fi - -if [ -z $read_hap1 ]; then - echo "No input provided. Exit." - exit -1 -fi - -if [ -z $scripts ]; then - echo "Setting \$scripts path" - scripts=`which meryl` - scripts=`dirname $scripts` - scripts="$scripts/../../scripts" - echo "Done!" -fi - -read=${read/.meryl/} # all read counts -read_hap1=${read_hap1/.meryl/} # pat specific mers with read counts -read_hap2=${read_hap2/.meryl/} # mat specific mers with read counts - -hap_hist=$name.hapmers.hist -cn_hist=$name.spectra-hap-cn.hist - -# Get histogram from all and inherited hap-mers -echo -e "Assembly\tkmer_multiplicity\tCount" > $hap_hist -meryl histogram $read.meryl | awk -v asm="read-total" '{print asm"\t"$0}' >> $hap_hist - -for hap in $read_hap1 $read_hap2 -do - meryl histogram $hap.meryl | awk -v asm="${hap}" '{print asm"\t"$0}' >> $hap_hist -done - -## Remove this line if R is properly installed ## -module load R # -################################################# - -echo " -$scripts/plot/plot_spectra_asm.R -f $hap_hist -o ${hap_hist/.hist/}" -$scripts/plot/plot_spectra_asm.R -f $hap_hist -o ${hap_hist/.hist/} - -for asm_fa in $asm1_fa $asm2_fa -do - asm=`echo $asm_fa | sed 's/.fasta$//g' | sed 's/.fa$//g' | sed 's/.fasta.gz$//g' | sed 's/.fa.gz$//g'` - if ! [[ "$(ls -A $asm.meryl 2> /dev/null )" ]]; then - echo "Generate meryl db for $asm_fa" - meryl count k=$k output $asm.meryl $asm_fa - fi - - # For each haplotype - for read_hap in $read_hap1 $read_hap2 - do - read_hap=${read_hap/.meryl/} - - echo -e "Start processing $read_hap from $asm\n" - - echo -e "Copies\tkmer_multiplicity\tCount" > $asm.$read_hap.$cn_hist - - # Convert counts to read_hap counts - meryl intersect output $read_hap.$asm.meryl $read_hap.meryl $asm.meryl - #meryl statistics $read_hap.$asm.meryl | head -n3 | tail -n1 | awk -v asm=$asm -v read=$read '{print asm"\t"read"\t"$1"\t"$2}' >> hap.counts - #meryl statistics $read_hap.$asm.meryl | head -n4 | tail -n1 | awk -v asm=$asm -v read=$read '{print asm"\t"read"\t"$1"\t"$2}' >> hap.counts - #meryl histogram $read_hap.$asm.meryl | awk -v asm=$asm '{print asm"\t"$0}' >> $read_hap.$hap_hist - #rm -r $read_hap.$asm.meryl - - # Copy-number histogram per haplotype specific k-mers - meryl difference output read.$read_hap.0.meryl $read_hap.meryl $asm.meryl - meryl histogram read.$read_hap.0.meryl | awk -v read=$read_hap '{print read"-only\t"$0}' >> $asm.$read_hap.$cn_hist - rm -r read.$read_hap.0.meryl - - for i in $(seq 1 4) - do - meryl intersect output read.$read_hap.$i.meryl $read_hap.$asm.meryl [ equal-to $i ${asm}.meryl ] - meryl histogram read.$read_hap.$i.meryl | awk -v cn=$i '{print cn"\t"$0}' >> $asm.$read_hap.$cn_hist - rm -r read.$read_hap.$i.meryl - done - - meryl intersect output read.$read_hap.gt$i.meryl $read_hap.$asm.meryl [ greater-than $i ${asm}.meryl ] - meryl histogram read.$read_hap.gt$i.meryl | awk -v cn=">$i" '{print cn"\t"$0}' >> $asm.$read_hap.$cn_hist - rm -r read.$read_hap.gt$i.meryl - rm -r $read_hap.$asm.meryl - - echo " - $scripts/plot/plot_spectra_cn.R -f $asm.$read_hap.$cn_hist -o $name.$asm.$read_hap" - $scripts/plot/plot_spectra_cn.R -f $asm.$read_hap.$cn_hist -o $name.$asm.$read_hap - - done - -done - -