Repository: Uauy-Lab/IBSpy (language: python)

Issue:
Solved compilation
There were two main problems:

1. Naming the C source as the Cython source is a big no-no. When Cython converts the code to C, it uses the base name of the .pyx file to name the .c file. In this case, this meant that the C file had been overwritten.
2. The Cython Extension, in the setup.py, needs to have the C sources included. You were creating two extensions, but the first was not linked properly to the second.


Discussion / hints:

Fix patch:
diff --git a/IBSplot/HapGMMbuild.py b/IBSplot/HapGMMbuild.py
deleted file mode 100755
index 1ae543a..0000000
--- a/IBSplot/HapGMMbuild.py
+++ /dev/null
@@ -1,104 +0,0 @@
-import argparse
-import pandas as pd
-import numpy as np
-from sklearn.mixture import GaussianMixture
-
-
-def read_db():
-	inDB = pd.read_csv(args.db_file, delimiter='\t')
-	return inDB
-
-#### prepare input data ####
-def prepare_data(args):
-	inDB = read_db()
-	filtered_inDB = inDB[inDB['variations'] <= int(args.varfltr)]
-	varDF = pd.DataFrame(filtered_inDB[['seqname','start','variations']])
-	varDF.reset_index(drop=True, inplace=True)
-	varArray = np.array(varDF['variations'])
-	log_varArray = np.log(varArray, where=(varArray != 0))
-	log_varArray = np.nan_to_num(log_varArray, nan=0.0)
-	log_varArray = log_varArray.reshape((len(log_varArray),1))
-	return log_varArray, varDF
-	
-#### build GMM model ####
-def GMM_build(args):
-	comNumber = int(args.comNumber)
-	model = GaussianMixture(n_components=comNumber, covariance_type='full')
-	log_varArray, varDF = prepare_data(args)
-	model.fit(log_varArray)
-	varGMM_predict = model.predict(log_varArray)
-	varGMM_predict = varGMM_predict.reshape((len(varGMM_predict),1)) 
-	varGMM_predict = pd.DataFrame({'gauss_mx_IBS': varGMM_predict[:, 0]})     
-	varDF['gauss_mx_IBS'] = varGMM_predict['gauss_mx_IBS']                                   
-	gmmIBS_DF = varDF.groupby(by=['gauss_mx_IBS'])
-	return gmmIBS_DF, varDF
-
-def transform_gmm(args):
-	gmmIBS_DF, varDF = GMM_build(args)
-	key_group = []
-	value_group = []
-	for key, group in gmmIBS_DF:
-	    key_group.append(key)
-	    value_group.append(group['variations'].median())
-	dic_grupby = dict(zip(key_group, value_group))
-	IBS_group = min(dic_grupby, key=dic_grupby.get)
-	varDF['gauss_mx_IBS'] = np.where(varDF['gauss_mx_IBS'] == IBS_group, 1, 0)
-	return varDF
-
-def gmm_by_chromosome():
-	inDB = transform_gmm(args)
-	byChrDf = inDB[inDB['seqname'].str.contains(args.chrNme)].copy()
-	byChrDf.reset_index(drop=True, inplace=True)
-	return byChrDf
-
-#### stitching haplotypes: number of non-IBS "outliers" that must a appear consecutively in a windows to be called non-IBS ####
-def stitch_haploBlocks(args):
-    byChrDf = gmm_by_chromosome()
-    gmmDta = byChrDf['gauss_mx_IBS'].copy()
-    StitchVarNum = int(args.StitchVarNum)
-    for i in range(len(gmmDta)):
-        if gmmDta[i] == 1:
-            if i < (len(gmmDta) - StitchVarNum):
-                count = 0
-                for n in range(1, (StitchVarNum+1)):
-                    if gmmDta[i+n] == 0:
-                        count += 1
-                if count == StitchVarNum:
-                    continue
-                else:
-                    gmmDta[i+1] = 1
-    hapBlock = gmmDta
-    hapBlock = pd.Series(hapBlock)
-    byChrDf['hap_block'] = hapBlock.values
-    # put back full dataset of data for final haplo blocks
-    hapCntFile = read_db()
-    hapCntFile = hapCntFile[['seqname', 'start', 'variations']].copy()
-    hapCntFile = hapCntFile[hapCntFile['seqname'].str.contains(args.chrNme)].copy()
-    byChrDf = pd.merge(hapCntFile, byChrDf, left_on='start', right_on='start', how='left')
-    byChrDf.loc[:,'gauss_mx_IBS':'hap_block'] = np.where(byChrDf.loc[:, 'gauss_mx_IBS':'hap_block'] == 1, 1, 0)
-    byChrDf.rename(columns={'seqname_x':'seqname', 'variations_x':'variations'}, inplace=True)
-    byChrDf = byChrDf[['seqname', 'start', 'variations', 'gauss_mx_IBS', 'hap_block']].copy()
-    return byChrDf
-
-def parse_arguments():
-	parser = argparse.ArgumentParser()
-
-	parser.add_argument('-db', '--db_file', help='Tab separated file with variations genetared by IBSpy output')
-	parser.add_argument('-rf', '--refId', help='Name of the genome reference used')
-	parser.add_argument('-qr', '--qryId', help='Name of the query sample')
-	parser.add_argument('-chr', '--chrNme', help='Chromosome name to be plotted')
-	parser.add_argument('-w', '--windSize', help='Windows size to count variations within')
-	parser.add_argument('-vf', '--varfltr', help='Filter variations above this threshold to compute GMM model')
-	parser.add_argument('-nc', '--comNumber', help='Number of componenets for the GMM model')
-	parser.add_argument('-st', '--StitchVarNum', help='Consecutive "outliers" in windows to stitch')
-	parser.add_argument('-gm','--output', help='GMM model output file grouped by IBS (1) and non-IBS (0)')
-	parser.add_argument('-ho','--hap_output', help='GMM model output file grouped by IBS (1) and non-IBS (0)')
-	
-	args = parser.parse_args()
-	return args
-
-args = parse_arguments()
-
-transform_gmm(args).to_csv(args.output, index=False,  sep='\t', compression="gzip")
-stitch_haploBlocks(args).to_csv(args.hap_output, index=False,  sep='\t', compression="gzip")
-
diff --git a/IBSplot/README.md b/IBSplot/README.md
deleted file mode 100755
index c6a5875..0000000
--- a/IBSplot/README.md
+++ /dev/null
@@ -1,143 +0,0 @@
-# IBSplot
-
-## Count variations using different windows sizes:
-
-	Use the script: var_byWind.py
-	Or adjust the bash script to your folder paths and arguments: run_var_byWind.sh
-
-### Run the code in the command line to generate a 500,000 bp windows file from the 50,000 bp original file:
-```sh
-	python var_byWind.py \
-		-db kmeribs-Wheat_Jagger-Flame.tsv.gz.gz \
-		-cl chr_sizes_jagger.genome.txt \
-		-w 500000 \
-		-o kmeribs-Wheat_Jagger-Flame_500000.tsv.gz
-```
-### Arguments description:
-```
-	-db [ output ".tsv.gz" file generated by the IBSpy analysis ]
-	-cl [ tab separated no header file with column one; chromosome names and column two; chromosome lengths]
-	-w [ size of the windows to count variations ]
-	-o [ file output ".tsv" with three columns "seqname", "start", "variations"; this can be used by the HapGMMbuild.py script exactly as described above ]
-```
-
-```
-	chromosome lengths example of the reference jagger: chr_sizes_jagger.genome.txt
-
-chr1A	596211899
-chr1B	705338699
-chr1D	493450010
-chr2A	804285258
-chr2B	810734643
-chr2D	673981989
-chr3A	754496630
-chr3B	855759449
-chr3D	614042580
-chr4A	743847818
-chr4B	673340788
-chr4D	515668560
-chr5A	713320806
-chr5B	703598484
-chr5D	570159854
-chr6A	625793224
-chr6B	721110502
-chr6D	459355444
-chr7A	745201001
-chr7B	749284433
-chr7D	643344672
-chrUn	379824255
-```
-
-
-## Build haplotypes using the 500,000 bp windows using the file created above:
-
-	Use the python script: HapGMMbuild.py
-	Or adjust the bash script to your folder paths and arguments: run_HapGMMbuild.sh
-
-### Run the code in the command line as the follow example for 500,000 bp windows and 1800 variation count filter using Jagger and Flame samples:
-NOTE: you will need to adjust the windows and filters according to the input data of the reference and query samples.
-
-```sh
-	python HapGMMbuild.py \
-		-db kmeribs-Wheat_Jagger-Flame_500000.tsv.gz \
-		-rf jagger \
-		-qr Flame \
-		-w 500000 \
-		-vf 1800 \
-		-nc 3 \
-		-chr chr1A \
-		-st 3 \
-		-gm gmm_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz \
-		-ho chr1A_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz
-
-```
-### Arguments description:
-```
-	-db [ output ".tsv.gz" file generated by the var_byWind.py above script ]
-	-rf [ genome reference used in the IBSpy (Jagger in this example) analysis ]
-	-qr [ name of the query or test sample used in the IBSpy analysis (Flame in this example) ]
-	-w [ size of the windows to count variations ]
-	-vf [ filter variations above this threshold to compute GMM model ]
-	-nc [ number of components used in the GMM model, we are mainly using 3, but this can be changed according to the samples used ]
-	-chr [ chromosome name to plot; must match the name in the "-cl" argument file, column one ]
-	-st [consecutive "outliers" in windows to stitch]
-	-gm [ output ".tsv.gz" file categorized in the column "gauss_mx_IBS" as a IBS == 1 and non-IBS == 0 ]
-	-ho [ haplotypes stitched using the GMM model output file grouped by IBS (1) and non-IBS (0) ]
-```
-
-
-## Generate histograms to visualize variation count distributions and categories generated by the GMM model:
-
-	Use the script: plot_HapGMM.py
-	Or adjust the bash script to your folder paths and arguments: run_histo_HapGMM.sh
-
-```sh
-	python histo_GMM.py \
-		-db kmeribs-Wheat_Jagger-Flame_500000.tsv.gz \
-		-gm gmm_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz \
-		-rf jagger \
-		-qr Flame \
-		-w 500000 \
-		-vf 1800 \
-		-do hist_Jagger_vs_Flame_500000_wd_1800_flt.pdf \
-		-go hist_Jagger_vs_Flame_500000_wd_whole_genome.pdf
-```
-### Arguments description:
-```
-Same as "HapGMMbuild.py", but with the additional following arguments: 
-	-do [ Output histogram of the variations from the whole genome in ".pdf" format ]
-	-ho [ Output histogram of the variations from the whole genome by IBS and non-IBS in ".pdf" format ]
-```
-
-
-## Plot final haploblocks:
-
-	Use the script: plot_HapGMM.py
-	Or adjust the bash script to your folder paths and arguments: run_plot_HapGMM.sh
-	
-### Run the code in the command line as follow to generate haplotypes using 500,000 bp windows and filter for 1800 count using jagger and flame samples:
-```sh
-	python plot_HapGMM.py \
-		-ho chr1A_Jagger_vs_Flame_500000_wd_1800_flt.tsv \
-		-cl chr_sizes_jagger.genome.txt \
-		-rf jagger \
-		-qr Flame \
-		-chr chr1A \
-		-w 500000 \
-		-vf 1800 \
-		-o hap_blocks_chr1A_Jagger_vs_Flame_500000_wd_1800_flt.jpg
-```
-
-### Arguments description:
-```
-	-ho [ file generated by the HapGMMbuild.py described above with the -gm command, which will have the "gauss_mx_IBS" column with 0/1 ]
-	-cl [ tab separated no header file with column one; chromosome names and column two; chromosome lengths (described above)]
-	-rf [ genome reference used in the IBSpy analysis ]
-	-qr [ name of the query or test sample used in the IBSpy analysis ]
-	-chr [ chromosome name to plot; must match the name in the "-cl" argument file, column one]
-	-w [ size of the windows to count variations; must match the windows used to compute the GMM model ]
-	-vf [ filter variations above this threshold to compute GMM model; must match the windows used to compute the GMM model ]
-	-o [ output scatter plot and bars with haplotypes in ".jpg" format ]
-
-```
-end...
\ No newline at end of file
diff --git a/IBSplot/examples/output/Jagger_vs_Flame/chr1A_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz b/IBSplot/examples/output/Jagger_vs_Flame/chr1A_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz
deleted file mode 100755
index e999876..0000000
Binary files a/IBSplot/examples/output/Jagger_vs_Flame/chr1A_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz and /dev/null differ
diff --git a/IBSplot/examples/output/Jagger_vs_Flame/gmm_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz b/IBSplot/examples/output/Jagger_vs_Flame/gmm_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz
deleted file mode 100755
index f0d2e3a..0000000
Binary files a/IBSplot/examples/output/Jagger_vs_Flame/gmm_Jagger_vs_Flame_500000_wd_1800_flt.tsv.gz and /dev/null differ
diff --git a/IBSplot/examples/output/Jagger_vs_Flame/hap_blocks_chr1A_Jagger_vs_Flame_500000_wd_1800_flt.jpg b/IBSplot/examples/output/Jagger_vs_Flame/hap_blocks_chr1A_Jagger_vs_Flame_500000_wd_1800_flt.jpg
deleted file mode 100755
index 59c2450..0000000
Binary files a/IBSplot/examples/output/Jagger_vs_Flame/hap_blocks_chr1A_Jagger_vs_Flame_500000_wd_1800_flt.jpg and /dev/null differ
diff --git a/IBSplot/examples/output/Jagger_vs_Flame/hist_Jagger_vs_Flame_500000_wd_1800_flt.pdf b/IBSplot/examples/output/Jagger_vs_Flame/hist_Jagger_vs_Flame_500000_wd_1800_flt.pdf
deleted file mode 100755
index 81a3265..0000000
Binary files a/IBSplot/examples/output/Jagger_vs_Flame/hist_Jagger_vs_Flame_500000_wd_1800_flt.pdf and /dev/null differ
diff --git a/IBSplot/examples/output/Jagger_vs_Flame/hist_Jagger_vs_Flame_500000_wd_whole_genome.pdf b/IBSplot/examples/output/Jagger_vs_Flame/hist_Jagger_vs_Flame_500000_wd_whole_genome.pdf
deleted file mode 100755
index 2f7b88d..0000000
Binary files a/IBSplot/examples/output/Jagger_vs_Flame/hist_Jagger_vs_Flame_500000_wd_whole_genome.pdf and /dev/null differ
diff --git a/IBSplot/examples/output/Jagger_vs_Flame/kmeribs-Wheat_Jagger-Flame_500000.tsv.gz b/IBSplot/examples/output/Jagger_vs_Flame/kmeribs-Wheat_Jagger-Flame_500000.tsv.gz
deleted file mode 100755
index 16dc0c6..0000000
Binary files a/IBSplot/examples/output/Jagger_vs_Flame/kmeribs-Wheat_Jagger-Flame_500000.tsv.gz and /dev/null differ
diff --git a/IBSplot/examples/run_HapGMMbuild.sh b/IBSplot/examples/run_HapGMMbuild.sh
deleted file mode 100755
index 0cf1da8..0000000
--- a/IBSplot/examples/run_HapGMMbuild.sh
+++ /dev/null
@@ -1,24 +0,0 @@
-#!/bin/bash
-
-REF=Jagger
-QUERY=Flame
-WINDOWS=500000
-FILTER=1800
-CHROM=ch1A
-
-scrip_dir='path/to/scripts'
-in_dir='path/to/database'
-out_dir='path/to/output_folder'
-
-
-python $scrip_dir/HapGMMbuild.v2.py \
-	-db $in_dir/kmeribs-Wheat_${REF}-${QUERY}_${WINDOWS}.tsv.gz \
-	-rf ${REF} \
-	-qr ${QUERY} \
-	-w ${WINDOWS} \
-	-vf ${FILTER} \
-	-nc 3 \
-	-chr $CHROM \
-	-st 3 \
-	-gm $out_dir/gmm_${REF}_vs_${QUERY}_${WINDOWS}_wd_${FILTER}_flt.tsv.gz \
-	-ho $out_dir/${CHROM}_${REF}_vs_${QUERY}_${WINDOWS}_wd_${FILTER}_flt.tsv.gz 
\ No newline at end of file
diff --git a/IBSplot/examples/run_histo_HapGMM.sh b/IBSplot/examples/run_histo_HapGMM.sh
deleted file mode 100755
index dedda5a..0000000
--- a/IBSplot/examples/run_histo_HapGMM.sh
+++ /dev/null
@@ -1,22 +0,0 @@
-#!/bin/bash
-
-
-REF=Jagger
-QUERY=Flame
-WINDOWS=500000
-FILTER=1800
-CHROM=chr1A
-
-scrip_dir='path/to/scripts'
-in_dir='path/to/database'
-out_dir='path/to/output_folder'
-
-python $scrip_dir/histo_GMM.py \
-	-db $in_dir/kmeribs-Wheat_${REF}-${QUERY}_${WINDOWS}.tsv.gz \
-	-gm $in_dir/gmm_${REF}_vs_${QUERY}_${WINDOWS}_wd_${FILTER}_flt.tsv.gz \
-	-rf ${REF} \
-	-qr ${QUERY} \
-	-w ${WINDOWS} \
-	-vf ${FILTER} \
-	-do $out_dir/hist_${REF}_vs_${QUERY}_${WINDOWS}_wd_${FILTER}_flt.pdf \
-	-go $out_dir/hist_${REF}_vs_${QUERY}_${WINDOWS}_wd_whole_genome.pdf
\ No newline at end of file
diff --git a/IBSplot/examples/run_plot_HapGMM.sh b/IBSplot/examples/run_plot_HapGMM.sh
deleted file mode 100755
index bd447b5..0000000
--- a/IBSplot/examples/run_plot_HapGMM.sh
+++ /dev/null
@@ -1,22 +0,0 @@
-#!/bin/bash
-
-REF=Jagger
-QUERY=Flame
-WINDOWS=500000
-FILTER=1800
-CHROM=chr1A
-
-scrip_dir='path/to/scripts'
-in_dir='path/to/database'
-out_dir='path/to/output_folder'
-chr_dir='path/to/genome_sizes_folder'
-
-python $scrip_dir/plot_HapGMM.v2.py \
-	-ho $in_dir/${CHROM}_${REF}_vs_${QUERY}_${WINDOWS}_wd_${FILTER}_flt.tsv.gz \
-	-cl $chr_dir/chr_sizes_jagger.genome.txt \
-	-rf ${REF} \
-	-qr ${QUERY} \
-	-chr ${CHROM} \
-	-w ${WINDOWS} \
-	-vf ${FILTER} \
-	-o $out_dir/hap_blocks_${CHROM}_${REF}_vs_${QUERY}_${WINDOWS}_wd_${FILTER}_flt.jpg
\ No newline at end of file
diff --git a/IBSplot/examples/run_var_byWind.sh b/IBSplot/examples/run_var_byWind.sh
deleted file mode 100755
index 03aa1f0..0000000
--- a/IBSplot/examples/run_var_byWind.sh
+++ /dev/null
@@ -1,17 +0,0 @@
-#!/bin/bash
-
-REF=Jagger
-QUERY=Flame
-IN_WINDOWS=50000
-OUT_WINDOWS=500000
-
-scrip_dir='path/to/scripts'
-in_dir='path/to/database'
-out_dir='path/to/output_folder'
-chr_dir='path/to/genome_sizes_folder'
-
-python $scrip_dir/var_byWind.v2.py \
-	-db $in_dir/kmeribs-Wheat_${REF}-${QUERY}.tsv.gz.gz \
-	-cl $chr_dir/chr_sizes_Jagger.genome.txt \
-	-w ${OUT_WINDOWS} \
-	-o $out_dir/kmeribs-Wheat_${REF}-${QUERY}_${OUT_WINDOWS}.tsv.gz
\ No newline at end of file
diff --git a/IBSplot/histo_GMM.py b/IBSplot/histo_GMM.py
deleted file mode 100755
index 008a073..0000000
--- a/IBSplot/histo_GMM.py
+++ /dev/null
@@ -1,76 +0,0 @@
-import argparse
-import pandas as pd
-import matplotlib.pyplot as plt
-from matplotlib.ticker import (MultipleLocator, NullFormatter, ScalarFormatter)
-
-def read_db():
-	inDB = pd.read_csv(args.db_file, delimiter='\t')
-	return inDB
-
-def read_gmm_f():
-    in_gmm_f = pd.read_csv(args.gmm_File, delimiter='\t')
-    return in_gmm_f
-
-# data
-def plot_histo_db(args):
-	inDB = read_db()
-	X_a = inDB.loc[inDB['variations']]
-	X_a = inDB['variations'].astype(int)
-
-	#plot
-	fig = plt.figure(figsize=(10, 2))
-	plt.hist(X_a, log=False, alpha=1, bins=range(min(X_a), max(X_a)), color='darkred', label='All-regions')
-	plt.title(f'All variations [{args.refId} vs {args.qryId}], {args.windSize}_windows', fontweight="bold", fontsize=12)
-	plt.legend(prop={'size': 15})
-	plt.grid(True)
-	plt.ylabel('Frequency', fontweight="bold", fontsize=11)
-	plt.xlabel('Variations count', fontweight="bold", fontsize=11)
-	plt.xscale('symlog')
-	plt.xlim(0,5000)
-	fig.axes[0].xaxis.set_major_formatter(ScalarFormatter())
-	plt.xticks([0,1,2,3,5,10,20,40,100,300,1000,2000,5000], fontsize = 10)
-	return fig
-
-#### GMM variations count histogram after classification ####
-# data
-def plot_histo_gmm(args):
-	varDF = read_gmm_f()
-	X_b = varDF.loc[varDF['gauss_mx_IBS'] == 1]
-	X_b = X_b['variations'].astype(int)
-	X_c = varDF.loc[varDF['gauss_mx_IBS'] == 0]
-	X_c = X_c['variations'].astype(int)
-
-	# plot
-	fig = plt.figure(figsize=(10, 2))
-	plt.hist(X_b, alpha=0.5, bins=range(min(X_b), max(X_b)), color='darkblue', label='IBS')
-	plt.hist(X_c, alpha=0.5, bins=range(min(X_c), max(X_c)), color='darkgreen', label='non-IBS')
-	plt.title(f'GMM: filtered variations <= {args.varfltr} [{args.refId} vs {args.qryId}], {args.windSize}_windows', fontweight="bold", fontsize=12)
-	plt.legend(prop={'size': 15})
-	plt.grid(True)
-	plt.ylabel('Frequency', fontweight="bold", fontsize=11)
-	plt.xlabel('Variations count', fontweight="bold", fontsize=11)
-	plt.xscale('symlog')
-	plt.xlim(0,5000)
-	fig.axes[0].xaxis.set_major_formatter(ScalarFormatter())
-	plt.xticks([0,1,2,3,5,10,20,40,100,300,1000,2000,5000], fontsize=10)
-	return fig
-
-def parse_arguments():
-	parser = argparse.ArgumentParser()
-
-	parser.add_argument('-db', '--db_file', help='Tab separated file with variations genetared by IBSpy output')
-	parser.add_argument('-gm', '--gmm_File', help='GMM model output file grouped by IBS (1) and non-IBS (0)')
-	parser.add_argument('-rf', '--refId', help='Name of the genome reference')
-	parser.add_argument('-qr', '--qryId', help='Name of the query sample')
-	parser.add_argument('-w', '--windSize', help='Windows size to count variations within')
-	parser.add_argument('-vf', '--varfltr', help='Filter variations above this threshold to compute GMM model')
-	parser.add_argument('-do','--db_output', help='Output histogram of the variations from the whole genome in ".pdf" format')
-	parser.add_argument('-go','--gmm_output', help='Output histogram of the variations from the whole genome by IBS and non-IBS in ".pdf" format')
-
-	args = parser.parse_args()
-	return args
-
-args = parse_arguments()
-
-plot_histo_db(args).savefig(args.db_output, format='pdf', bbox_inches='tight')
-plot_histo_gmm(args).savefig(args.gmm_output, format='pdf', bbox_inches='tight')
diff --git a/IBSplot/plot_HapGMM.py b/IBSplot/plot_HapGMM.py
deleted file mode 100755
index 09d1116..0000000
--- a/IBSplot/plot_HapGMM.py
+++ /dev/null
@@ -1,102 +0,0 @@
-import argparse
-import pandas as pd
-import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib.ticker import (MultipleLocator, NullFormatter, ScalarFormatter)
-
-def read_hapBlocks_f():
-    hoFile = pd.read_csv(args.hoFile, delimiter='\t')
-    return hoFile
-def read_chrLen_f():
-    in_chrLen = pd.read_csv(args.chrLen_f, delimiter='\t', header=None)
-    return in_chrLen
-
-def plot_haploBlocks(args):
-    hoFile = read_hapBlocks_f()
-    IBS = hoFile[hoFile['hap_block'] == 1]
-    strtIbsCln = np.array(IBS['start'])
-    hapBlkIbsCln = np.array(IBS['hap_block'])
-    hapArrIbsCln = np.vstack((strtIbsCln, hapBlkIbsCln))
-
-    non_IBS = hoFile[hoFile['hap_block'] == 0].copy()
-    chr_st_non_IBS = np.array(non_IBS['start'])
-    hp_blk_non_IBS = np.array(non_IBS['hap_block'])
-    chr_hp_arr_non_IBS = np.vstack((chr_st_non_IBS, hp_blk_non_IBS))
-
-    #### GMM before stitching ####
-    IBS_org = hoFile[hoFile['gauss_mx_IBS'] == 1].copy()
-    strtIbsCln_org = np.array(IBS_org['start'])
-    IBS_org = np.array(IBS_org['gauss_mx_IBS'])
-    hapArrIbsCln_org = np.vstack((strtIbsCln_org, IBS_org))
-
-    non_IBS_org = hoFile[hoFile['gauss_mx_IBS'] == 0].copy()
-    chr_st_non_IBS_org = np.array(non_IBS_org['start'])
-    non_IBS_org = np.array(non_IBS_org['gauss_mx_IBS'])
-    chr_hp_arr_non_IBS_org = np.vstack((chr_st_non_IBS_org, non_IBS_org))
-
-    #### base subplots ####
-    fig, axs = plt.subplots(3, figsize=(20, 7), sharex=True, gridspec_kw={'height_ratios': [1, 1,3]})
-    fig.suptitle(f'{args.chrNme}: {args.refId} vs {args.qryId}, {args.windSize} bp windows, filter: {args.varfltr}', fontweight="bold", fontsize=25)
-
-    #### horizontal bar plot ####
-    ### GMM_stitched
-    axs[0].eventplot(hapArrIbsCln, colors='darkcyan', alpha=1, lineoffsets=[1,1], linelengths=1.5)
-    axs[0].eventplot(chr_hp_arr_non_IBS, colors='goldenrod', lineoffsets=[1,1], linelengths=1)
-
-    ### GMM_original (before stitching) ###
-    axs[1].eventplot(hapArrIbsCln_org, colors='darkcyan', alpha=1, lineoffsets=[1,1], linelengths=1.5)
-    axs[1].eventplot(chr_hp_arr_non_IBS_org, colors='goldenrod', lineoffsets=[1,1], linelengths=1)
-
-    ### scatter plot ###
-    IBS_block = hoFile[(hoFile['hap_block'] == 1)]
-    x = IBS_block['start']
-    y = IBS_block['variations']
-    axs[2].scatter(x,y, c="darkcyan", alpha=0.5, s=4, label='IBS')
-
-    non_IBS_block = hoFile[(hoFile['hap_block'] == 0)]
-    a = non_IBS_block['start']
-    b = non_IBS_block['variations']
-    axs[2].scatter(a,b, c="goldenrod", alpha=0.6, s=4, marker='>', label='non-IBS')
-    axs[2].set_yscale('symlog')
-    axs[2].set_ylabel('Variations count', fontweight="bold", fontsize=17)
-    axs[2].set_xlabel('Chromosome position [Mbp]', fontweight="bold", fontsize=17)
-    axs[2].yaxis.set_major_formatter(ScalarFormatter())
-
-    ### get chromosome length to plot X_axis limits ###
-    chrLen_f = read_chrLen_f()
-    chr_length = int(chrLen_f[chrLen_f[0] == args.chrNme][1].values) + 10000000
-
-    #### labels and axis limits ####
-    names = ['gmm_stitched', 'gmm_origin']
-    for axs, names in zip(axs, names):
-            axs.text(-0.01, 0.5, names, va='center', ha='right', fontweight= "bold", fontsize=17, transform=axs.transAxes)
-            plt.grid(True, color='gray', linestyle='--', linewidth=0.5)
-            
-            axs.set_xlim(0,chr_length)
-            x_labels = list(range(0,chr_length,20000000)) # one label every 20 Mb
-            list_labels = int(chr_length/1000000) # transform to Mb to reduce number size
-            x_ticklabels = list(range(0,list_labels,20)) # use tranformed labels every 20 Mb
-            plt.xticks(x_labels, x_ticklabels, rotation=45, fontsize=15)
-            axs.yaxis.set_major_formatter(ScalarFormatter())
-            plt.yticks([1,2,3,10,20,40,100,300,1000,2000,5000], fontsize = 15)
-            plt.legend(prop={'size': 10}, markerscale=4, borderaxespad=0.3, handletextpad=0.1)
-            axs.set_axis_off()
-    return fig
-
-def parse_arguments():
-    parser = argparse.ArgumentParser()
-
-    parser.add_argument('-ho', '--hoFile', help='Tab separated file with IBS and non_IBS categorized by GMM model')
-    parser.add_argument('-rf', '--refId', help='Name of the genome reference used')
-    parser.add_argument('-qr', '--qryId', help='Name of the query sample')
-    parser.add_argument('-chr', '--chrNme', help='Chromosome name to be plotted')
-    parser.add_argument('-cl', '--chrLen_f', help='Reference chromosome lenghts file')
-    parser.add_argument('-w', '--windSize', help='Windows size to count variations within')
-    parser.add_argument('-vf', '--varfltr', help='Filter variations above this threshold to compute GMM model')
-    parser.add_argument('-o','--outImg', help='Output scatter plot and bars with haplotypes in ".jpg" format')
-    
-    args = parser.parse_args()
-    return args
-
-args = parse_arguments()
-plot_haploBlocks(args).savefig(args.outImg, format='jpg', dpi=150, bbox_inches='tight')
diff --git a/IBSpy/IBSPy_results.py b/IBSpy/IBSPy_results.py
index 8ac234b..c205758 100755
--- a/IBSpy/IBSPy_results.py
+++ b/IBSpy/IBSPy_results.py
@@ -1,48 +1,104 @@
 import pandas as pd
+import numpy as np
+from sklearn.mixture import GaussianMixture
 
 class IBSpyResults:
-    def __init__(self, filename, chromosome_lengths, window_size):
-        self.filename = filename
+    # class variables go here
+    def __init__(self, filename, window_size, filter_counts):
+
         self.db = pd.read_csv(filename, delimiter='\t')
-        self.lengths =  pd.read_csv(chromosome_lengths, delimiter='\t', header=None)
         self.window_size = window_size
+        self.filter_counts = filter_counts
 
-
-#### file with chromosome lengths ####
     def count_by_windows(self):
-        windSize = int(self.window_size)
-        db_ByGenom = pd.DataFrame()
-        for index, row in self.lengths.iterrows():
-            ## get individual chromosome and size
-            chrID = self.lengths[self.lengths[0] == row[0]]
-            chrLen = chrID.iloc[0][1]
-            ## access db data 
-            db_DF = self.db[self.db['seqname'].str.contains(row[0])] # made a change to fit new files from IBSpy
-            db_DF = db_DF[['seqname','start','variations']]
-
-            ## empty DF to merge by chromosome
-            w_pos = 0
-            db_byChr = pd.DataFrame()
-            while w_pos <= chrLen:
-                db_DF_ByWind = db_DF[(db_DF['start'] > w_pos) & (db_DF['start'] <= w_pos + windSize)]
-                db_DF_ByWind = db_DF_ByWind.loc[:, ('seqname','start','variations')]
-                db_DF_ByWind['w_num'] = w_pos + windSize
-                w_pos += windSize
-                db_byChr = db_byChr.append(db_DF_ByWind)
+        windSize = self.window_size
+        db_DF = self.db[['seqname','start','end','variations']]
+        # get the longest chromosome
+        chrLen = db_DF['end'].max()
+        # windows to iterate
+        w_pos = 0
+        db_byChr = pd.DataFrame()
+        while w_pos <= chrLen:
+            db_DF_ByWind = db_DF[(db_DF['start'] >= w_pos) & (db_DF['start'] < w_pos + windSize)]
+            db_DF_ByWind = db_DF_ByWind.loc[:, ('seqname','start','variations')]
+            db_DF_ByWind['window'] = w_pos + windSize
+            w_pos += windSize
+            db_byChr = db_byChr.append(db_DF_ByWind)
 
-            db_byChr = db_byChr.groupby(['w_num']).sum().reset_index()
-            db_byChr = db_byChr[['w_num', 'variations']]
-            db_byChr['seqname'] = row[0]
-            db_byChr = db_byChr[['seqname', 'w_num', 'variations']]
-            db_byChr.rename(columns={'w_num':'start'}, inplace=True)
-            db_ByGenom = db_ByGenom.append(db_byChr)
-        self.tmp_table = db_ByGenom
-        return db_ByGenom
+        # calculate statistics by chromosome, windows, and variations
+        out_db = db_byChr.groupby(['seqname','window'])['variations'].sum().reset_index()
+        out_db['mean'] = db_byChr.groupby(['seqname','window'])['variations'].mean().values.round(2)
+        out_db['median'] = db_byChr.groupby(['seqname','window'])['variations'].median().values.round(2)
+        out_db['variance'] = db_byChr.groupby(['seqname','window'])['variations'].var().values.round(2)
+        out_db['std'] = db_byChr.groupby(['seqname','window'])['variations'].std().values.round(2)
+        return out_db
+    
+    def transform_counts_to_log(self):
+        by_windows_db = self.count_by_windows()
+        if self.filter_counts is not None:
+            applied_filter = by_windows_db[by_windows_db['variations'] <= self.filter_counts]
+        else:
+            applied_filter = by_windows_db
+        varDF = pd.DataFrame(applied_filter[['seqname','window','variations']])
+        varDF.reset_index(drop=True, inplace=True)
+        varArray = np.array(varDF['variations'])
+        log_varArray = np.log(varArray, where=(varArray != 0))
+        log_varArray = np.nan_to_num(log_varArray, nan=0.0)
+        log_varArray = log_varArray.reshape((len(log_varArray),1))
+        return log_varArray, varDF
 
-    def transform_gmm(self):
-        return None
+    def build_gmm_model(self, n_components, covariance_type):
+        log_varArray, varDF = self.transform_counts_to_log()
+        model = GaussianMixture(n_components=n_components, covariance_type=covariance_type)
+        model.fit(log_varArray)
+        varGMM_predict = model.predict(log_varArray)
+        varGMM_predict = varGMM_predict.reshape((len(varGMM_predict),1)) 
+        varGMM_predict = pd.DataFrame({'v_gmm': varGMM_predict[:, 0]})     
+        varDF['v_gmm'] = varGMM_predict['v_gmm']                                   
+        gmmIBS_DF = varDF.groupby(by=['v_gmm'])
+        # build binary
+        key_group = []
+        value_group = []
+        for key, group in gmmIBS_DF:
+            key_group.append(key)
+            value_group.append(group['variations'].median())
+        dic_grupby = dict(zip(key_group, value_group))
+        IBS_group = min(dic_grupby, key=dic_grupby.get)
+        varDF['v_gmm'] = np.where(varDF['v_gmm'] == IBS_group, 1, 0)
+        
+        self.tmp_varDF = varDF
+        return varDF
 
+    def stitch_gmm_haplotypes(self, n_components, covariance_type, stitch_number):
+        stitch_db = self.build_gmm_model(n_components, covariance_type)
+        gmmDta = stitch_db['v_gmm'].copy()
 
+        StitchVarNum = stitch_number
+        for i in range(len(gmmDta)):
+            if gmmDta[i] == 1:
+                if i < (len(gmmDta) - StitchVarNum):
+                    count = 0
+                    for n in range(1, (StitchVarNum + 1)):
+                        if gmmDta[i+n] == 0:
+                            count += 1
+                    if count == StitchVarNum:
+                        continue
+                    else:
+                        gmmDta[i+1] = 1
+        hapBlock = gmmDta
+        hapBlock = pd.Series(hapBlock)
+        stitch_db['vh_block'] = hapBlock.values
+        # put back filtered data with haploblocks
+        hapCntFile = self.count_by_windows()
+        stitch_db = pd.merge(hapCntFile, stitch_db, left_on=['seqname','window'], right_on=['seqname','window'], how='left')
+        stitch_db.loc[:,'v_gmm':'vh_block'] = np.where(stitch_db.loc[:, 'v_gmm':'vh_block'] == 1, 1, 0)
+        stitch_db.rename(columns={'seqname_x':'seqname', 'variations_x':'variations'}, inplace=True)
+        stitch_db = stitch_db.drop(['variations_y'], axis=1)
+        return stitch_db
+    
     def run_analysis(self):
-        self.count_by_windows()
-        self.transform_gmm()
\ No newline at end of file
+        counts     = self.count_by_windows()
+        # by_windows_db = self.count_by_windows()
+        # normalised = self.transform_counts_to_log(by_windows_db)
+#         self.build_gmm_model()
+#         self.stitch_gmm_haplotypes()
\ No newline at end of file
diff --git a/IBSpy/IBSpy_plot.py b/IBSpy/IBSpy_plot.py
new file mode 100644
index 0000000..cc60dee
--- /dev/null
+++ b/IBSpy/IBSpy_plot.py
@@ -0,0 +1,114 @@
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.ticker import ScalarFormatter
+
+class IBSpyPlots:
+	# class variables go here
+	figure_size = (20, 5)
+	X_hist_lim = (0,10000)
+	v_count_ticks = [0,1,2,5,10,20,50,100,200,400,1000,2000,5000,10000]
+	ax_font_size = 17
+	axx_tick_size = 10
+	axy_tick_size = 15
+
+	def __init__(self, gmm_results, window_size, filter_counts, r_name, q_name):
+		
+		self.gmm_results = gmm_results
+		self.window_size = window_size
+		self.filter_counts = filter_counts
+		self.r_name = r_name
+		self.q_name = q_name
+
+	# raw_data
+	def raw_histogram(self):
+		inDB = self.gmm_results
+		X_a = inDB['variations']
+
+		fig = plt.figure(figsize= (20, 4))
+		plt.hist(X_a, log=False, alpha=1, bins=range(min(X_a), max(X_a)), color='darkred', label='All-counts')
+		plt.title(f'Raw variations count, whole genome: {self.r_name} vs {self.q_name}, {self.window_size}_windows', fontweight="bold", fontsize=22)
+		plt.legend(prop={'size': 15})
+		plt.grid(True)
+		plt.ylabel('Frequency', fontweight="bold", fontsize=self.ax_font_size)
+		plt.xlabel('Variations count', fontweight="bold", fontsize=self.ax_font_size)
+		plt.xscale('symlog')
+		plt.xlim(self.X_hist_lim)
+		fig.axes[0].xaxis.set_major_formatter(ScalarFormatter())
+		plt.xticks(self.v_count_ticks, fontsize = self.axx_tick_size)
+		plt.yticks(fontsize = self.axy_tick_size)
+		return fig
+	
+	# GMM
+	def gmm_histogram(self):
+		inDB = self.gmm_results
+		X_b = inDB['variations'][inDB['v_gmm'] == 1]
+		X_c = inDB['variations'][inDB['v_gmm'] == 0]
+
+		fig = plt.figure(figsize=(20, 4))
+		plt.hist(X_b, alpha=0.5, bins=range(min(X_b), max(X_b)), color='darkblue', label='IBS')
+		plt.hist(X_c, alpha=0.5, bins=range(min(X_c), max(X_c)), color='darkgreen', label='non-IBS')
+		plt.title(f'GMM variations count <= {self.filter_counts}, whole genome: {self.r_name} vs {self.q_name}, {self.window_size}_windows', fontweight="bold", fontsize=22)
+		plt.legend(prop={'size': 15})
+		plt.grid(True)
+		plt.ylabel('Frequency', fontweight="bold", fontsize=self.ax_font_size)
+		plt.xlabel('Variations count', fontweight="bold", fontsize=self.ax_font_size)
+		plt.xscale('symlog')
+		plt.xlim(self.X_hist_lim)
+		fig.axes[0].xaxis.set_major_formatter(ScalarFormatter())
+		plt.xticks(self.v_count_ticks, fontsize=self.axx_tick_size)
+		plt.yticks(fontsize = self.axy_tick_size)
+		return fig
+
+	def scatterplot_by_chromosome(self, chromosome_id):
+		variations_db = self.gmm_results
+		variations_db = variations_db[variations_db['seqname'].str.contains(chromosome_id)]
+		
+		IBS = variations_db[variations_db['vh_block'] == 1]
+		IBS_window_position = np.array(IBS['window'])
+		IBS_haplo_block = np.array(IBS['vh_block'])
+		IBS_window_block_array = np.vstack((IBS_window_position, IBS_haplo_block))
+
+		non_IBS = variations_db[variations_db['vh_block'] == 0]
+		non_IBS_window_position = np.array(non_IBS['window'])
+		non_IBS_haplo_block = np.array(non_IBS['vh_block'])
+		non_IBS_window_block_array = np.vstack((non_IBS_window_position, non_IBS_haplo_block))
+
+		# base subplots
+		fig, axs = plt.subplots(2, figsize=(20, 4), sharex=True, gridspec_kw={'height_ratios': [1, 3]})
+		plt.subplots_adjust(hspace=False)
+		fig.suptitle(f'{chromosome_id}: {self.r_name} vs {self.q_name}, {self.window_size} bp windows, filter: {self.filter_counts}', fontweight="bold", fontsize=22)
+
+		# horizontal bar plot
+		axs[0].eventplot(IBS_window_block_array, colors='darkcyan', alpha=1, lineoffsets=[1,1], linelengths=1.5)
+		axs[0].eventplot(non_IBS_window_block_array, colors='goldenrod', lineoffsets=[1,1], linelengths=1)
+
+		# scatter plot
+		x = IBS['window']
+		y = IBS['variations']
+		axs[1].scatter(x,y, c="darkcyan", alpha=0.5, s=4, label='IBS')
+		a = non_IBS['window']
+		b = non_IBS['variations']
+		axs[1].scatter(a,b, c="goldenrod", alpha=0.6, s=4, marker='>', label='non-IBS')
+
+		axs[1].set_yscale('symlog')
+		axs[1].set_ylabel('Variations count', fontweight="bold", fontsize = self.ax_font_size)
+		axs[1].set_xlabel('Chromosome position [Mbp]', fontweight="bold", fontsize = self.ax_font_size)
+		axs[1].yaxis.set_major_formatter(ScalarFormatter())
+		axs[1].set_ylim(1,10000)
+
+		# get chromosome length to plot X_axis limits
+		chr_length = variations_db['window'].max() + 2000000
+		names = ['GMM']
+		for axs, names in zip(axs, names):
+			axs.text(-0.01, 0.5, names, va='center', ha='right', fontweight= "bold", fontsize = self.ax_font_size, transform=axs.transAxes)
+			plt.grid(True, color='gray', linestyle='--', linewidth=0.5)
+			axs.set_xlim(0,chr_length)
+			x_labels = list(range(0,chr_length,10000000)) # one label every 20 Mb
+			list_labels = int(chr_length/1000000) # transform to Mb to reduce number size
+			x_ticklabels = list(range(0,list_labels,10)) 
+			plt.xticks(x_labels, x_ticklabels, rotation=45, fontsize = self.axx_tick_size)
+			plt.yticks(self.v_count_ticks, fontsize = self.axy_tick_size)
+			plt.legend(prop={'size': 10}, markerscale=4, borderaxespad=0.3, handletextpad=0.1)
+			axs.set_axis_off()
+		return fig
diff --git a/IBSpy/IBSpy_plot_window_results.py b/IBSpy/IBSpy_plot_window_results.py
index 736d690..8bfdd17 100644
--- a/IBSpy/IBSpy_plot_window_results.py
+++ b/IBSpy/IBSpy_plot_window_results.py
@@ -1,21 +1,55 @@
 import argparse
-from IBSpy import IBSpyResults
+from IBSpy import IBSpyResults, IBSpyPlots
+from matplotlib.backends.backend_pdf import PdfPages
+import matplotlib.pyplot as plt
+
 def parse_arguments():
     parser = argparse.ArgumentParser()
 
-    parser.add_argument('-d', '--IBSpy_counts',  required=True, help='Tab separated file with variations genetared by IBSpy output')
-    parser.add_argument('-c', '--chromosome_length',  required=True, help='Reference chromosome lenghts file')
-    parser.add_argument('-w', '--window_size',  required=True, help='Windows size to count variations within')
-    parser.add_argument('-o', '--output', help='File with variations by windows only')
+    parser.add_argument('-i', '--IBSpy_counts', help='tvs file genetared by IBSpy output')
+    parser.add_argument('-w', '--window_size', type=int, help='Windows size to count variations within')
+    parser.add_argument('-f', '--filter_counts', default=None, type=int, help='Filter number of variaitons above this threshold to compute GMM model, default=None')
+    parser.add_argument('-n', '--n_components', default=3, type=int, help='Number of componenets for the GMM model, default=3')
+    parser.add_argument('-c', '--covariance_type', default='full', help='type of covariance used for GMM model, default="full"')
+    parser.add_argument('-s', '--stitch_number', default=3, type=int, help='Consecutive "outliers" in windows to stitch, default=3')
+    parser.add_argument('-o', '--output', help='tsv file with variations count by windows and summary statistics')
+
+    # plot arguments
+    parser.add_argument('-r', '--reference', help='genome reference name')
+    parser.add_argument('-q', '--query', help='query sample')
+    parser.add_argument('-p','--plot_output', help=' histograms and ascatter files in .PDF format')
 
     args = parser.parse_args()
     return args
 
 def main():
-	args = parse_arguments()
-	results = IBSpyResults(args.IBSpy_counts, args.chromosome_length, args.window_size)
-	results_pd = results.count_by_windows()
-	results_pd.to_csv(args.output, index=False,  sep='\t', compression="gzip")
+    args = parse_arguments()
+
+    # count variations by windows from IBS_results
+    IBS_results = IBSpyResults(args.IBSpy_counts, args.window_size, args.filter_counts)
+    IBS_results = IBS_results.stitch_gmm_haplotypes(args.n_components, args.covariance_type, args.stitch_number)
+    IBS_results.to_csv(args.output, index=False,  sep='\t', compression="gzip")
+
+
+    plots =  IBSpyPlots(IBS_results, args.window_size, args.filter_counts, args.reference, args.query)
+    with PdfPages(args.plot_output) as pdf:
+
+        # histograms
+        raw_hist = plots.raw_histogram()
+        pdf.savefig(bbox_inches='tight')
+
+        gmm_hist = plots.gmm_histogram()
+        pdf.savefig(bbox_inches='tight')
+        
+        plt.close(raw_hist)
+        plt.close(gmm_hist)
+
+        # scatter plot
+        chromosomes = IBS_results["seqname"].str.split('_').str[0].unique()
+        for chromosome in chromosomes:
+            scatter_by_chr = plots.scatterplot_by_chromosome(chromosome)
+            pdf.savefig(bbox_inches='tight')
+            plt.close(scatter_by_chr)
 
 if __name__ == '__main__':
 	main()
\ No newline at end of file
diff --git a/IBSpy/__init__.py b/IBSpy/__init__.py
index d460a64..c9d120f 100644
--- a/IBSpy/__init__.py
+++ b/IBSpy/__init__.py
@@ -3,6 +3,7 @@
 from .jellyfish_kmer_database import JellyfishSDB, JellyfishBuilder
 from .IBSpy_window_count import window_count, parse_arguments
 from .IBSpy_results import IBSpyResults
+from .IBSpy_plot import IBSpyPlots
 import argparse
 
 def main(): 
diff --git a/README.md b/README.md
index f43ebed..c9bc54e 100644
--- a/README.md
+++ b/README.md
@@ -33,7 +33,7 @@ IBSpy requires to have a kmer database from the sequencing files. Currently two
 
 ## Running IBSPy
 
-IBSpy has relatively few options, you can look at them with the ```---help``` command. 
+IBSpy has relatively few options, you can look at them with the ```--help``` command. 
 
 ```sh
 IBSPy --help
@@ -59,8 +59,60 @@ optional arguments:
 
 To generate the table with the number of observed kmers and variants run the following command, using the kmer database from kmerGWAS use the following command:
 
+
 ```sh
  IBSpy --output "kmer_windows_LineXXX.tsv.gz" -z --database kmers_with_strand  --reference arinaLrFor.fa --window_size 50000 --compress --database_format kmerGWAS
 ```
 
+## Running IBSplot
+
+Look at the IBSplot commands using ```--help```.
+
+```sh
+IBSPy --help
+usage: IBSplot [-h] [-i IBSPY_COUNTS] [-w WINDOW_SIZE] [-f FILTER_COUNTS]
+               [-n N_COMPONENTS] [-c COVARIANCE_TYPE] [-s STITCH_NUMBER]
+               [-o OUTPUT] [-r REFERENCE] [-q QUERY] [-p PLOT_OUTPUT]
+
+optional arguments:
+  -h, --help            show this help message and exit
+  -i IBSPY_COUNTS, --IBSpy_counts IBSPY_COUNTS
+                        tvs file genetared by IBSpy output
+  -w WINDOW_SIZE, --window_size WINDOW_SIZE
+                        Windows size to count variations within
+  -f FILTER_COUNTS, --filter_counts FILTER_COUNTS
+                        Filter number of variaitons above this threshold to
+                        compute GMM model, default=None
+  -n N_COMPONENTS, --n_components N_COMPONENTS
+                        Number of componenets for the GMM model, default=3
+  -c COVARIANCE_TYPE, --covariance_type COVARIANCE_TYPE
+                        type of covariance used for GMM model, default="full"
+  -s STITCH_NUMBER, --stitch_number STITCH_NUMBER
+                        Consecutive "outliers" in windows to stitch, default=3
+  -o OUTPUT, --output OUTPUT
+                        tsv file with variations count by windows and summary
+                        statistics
+  -r REFERENCE, --reference REFERENCE
+                        genome reference name
+  -q QUERY, --query QUERY
+                        query sample
+  -p PLOT_OUTPUT, --plot_output PLOT_OUTPUT
+                        histograms and ascatter files in .PDF format
+```
+
+IBSplot uses the output table generated by IBSpy described above (e.g., ```"kmer_windows_LineXXX.tsv.gz"```). It can be used to count variant assigning larger windows. In the example below it is using 400,000 bp windows to compute  a GMM model and generate the plots.
+
+To generate the table with variant count categorized by the GMM model as IBS or non-IBS and generate the plots, run the following command:
+The description of the GMM model is [here](https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html#sklearn.mixture.GaussianMixture)
+
+```sh
+# minimal arguments
+IBSplot --IBSpy_counts "kmeribs-Wheat_Jagger-Flame.tsv.gz" --window_size 400000 --output gmm_ibs.tsv.gz --reference Jagger --query Flame --plot_output gmm_plots.pdf
+```
+
+In addition, you can include some or all of the following commands to tune the GMM model parameters and define the best IBS and non-IBS according to the reference and query sample used:
+
+```sh
+IBSplot --filter_counts 1000 --n_components 3 --covariance_type 'full' --stitch_number 3
+```
 
diff --git a/requirements.txt b/requirements.txt
index 6fd8c19..67ab208 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,7 @@
 cython
 biopython
-pyfaidx
\ No newline at end of file
+pyfaidx
+sklearn
+numpy
+pandas
+matplotlib
\ No newline at end of file


Test patch:
diff --git a/IBSpy/test_bash_scrip b/IBSpy/test_bash_scrip
new file mode 100644
index 0000000..d1d40c8
--- /dev/null
+++ b/IBSpy/test_bash_scrip
@@ -0,0 +1,12 @@
+IBSplot -i test_kmeribs-Wheat_Jagger-Flame.tsv \
+-c chr_sizes_Jagger.genome.txt \
+-w 200000 \
+-f 500 \
+-n 3 \
+-v full \
+-s 3 \
+-o file_out_test.tsv \
+-r Jagger \
+-q Flame \
+-do my_first_file.jpg
+-go GMM_my_first_file.pdf
diff --git a/tests/data/chr_sizes_jagger.genome.txt b/tests/data/chr_sizes_jagger.genome.txt
deleted file mode 100755
index 1f859f4..0000000
--- a/tests/data/chr_sizes_jagger.genome.txt
+++ /dev/null
@@ -1,22 +0,0 @@
-chr1A	596211899
-chr1B	705338699
-chr1D	493450010
-chr2A	804285258
-chr2B	810734643
-chr2D	673981989
-chr3A	754496630
-chr3B	855759449
-chr3D	614042580
-chr4A	743847818
-chr4B	673340788
-chr4D	515668560
-chr5A	713320806
-chr5B	703598484
-chr5D	570159854
-chr6A	625793224
-chr6B	721110502
-chr6D	459355444
-chr7A	745201001
-chr7B	749284433
-chr7D	643344672
-chrUn	379824255
diff --git a/tests/data/test4B.jagger.temp.kmerGWAS_k31 b/tests/data/test4B.jagger.temp.kmerGWAS_k31
new file mode 100644
index 0000000..d8e711e
Binary files /dev/null and b/tests/data/test4B.jagger.temp.kmerGWAS_k31 differ
diff --git a/tests/plots_data/gmm_ibs.tsv.gz b/tests/plots_data/gmm_ibs.tsv.gz
new file mode 100644
index 0000000..abca1cb
Binary files /dev/null and b/tests/plots_data/gmm_ibs.tsv.gz differ
diff --git a/tests/plots_data/gmm_plots.pdf b/tests/plots_data/gmm_plots.pdf
new file mode 100644
index 0000000..afceb41
Binary files /dev/null and b/tests/plots_data/gmm_plots.pdf differ
diff --git a/tests/data/kmeribs-Wheat_Jagger-Flame.tsv.gz b/tests/plots_data/kmeribs-Wheat_Jagger-Flame.tsv.gz
similarity index 100%
rename from tests/data/kmeribs-Wheat_Jagger-Flame.tsv.gz
rename to tests/plots_data/kmeribs-Wheat_Jagger-Flame.tsv.gz
diff --git a/tests/data/test_chr_sizes_jagger.tsv b/tests/plots_data/test_chr_sizes_jagger.tsv
similarity index 100%
rename from tests/data/test_chr_sizes_jagger.tsv
rename to tests/plots_data/test_chr_sizes_jagger.tsv
diff --git a/tests/data/test_kmeribs-Wheat_Jagger-Flame.tsv b/tests/plots_data/test_kmeribs-Wheat_Jagger-Flame.tsv
similarity index 100%
rename from tests/data/test_kmeribs-Wheat_Jagger-Flame.tsv
rename to tests/plots_data/test_kmeribs-Wheat_Jagger-Flame.tsv
diff --git a/tests/plots_data/test_run_arguments.sh b/tests/plots_data/test_run_arguments.sh
new file mode 100644
index 0000000..cd2e139
--- /dev/null
+++ b/tests/plots_data/test_run_arguments.sh
@@ -0,0 +1,19 @@
+# minimal arguments
+IBSplot --IBSpy_counts "kmeribs-Wheat_Jagger-Flame.tsv.gz" \
+--window_size 400000 \
+--output gmm_ibs.tsv.gz \
+--reference Jagger \
+--query Flame \
+--plot_output gmm_plots.pdf
+
+# all arguments
+IBSplot --IBSpy_counts "kmeribs-Wheat_Jagger-Flame.tsv.gz" \
+--window_size 400000 \
+--output gmm_ibs.tsv.gz \
+--reference Jagger \
+--query Flame \
+--plot_output gmm_plots.pdf \
+--filter_counts 1000 \
+--n_components 3 \
+--covariance_type 'full' \
+--stitch_number 3
\ No newline at end of file
diff --git a/tests/test_results.py b/tests/test_results.py
index cadba53..97e6ae6 100644
--- a/tests/test_results.py
+++ b/tests/test_results.py
@@ -1,5 +1,7 @@
 import os
 import sys
+import numpy as np
+
 
 import logging
 import operator
@@ -7,41 +9,64 @@
 from IBSpy import IBSpyResults 
 
 class TestResults(unittest.TestCase):
-	logger = logging.getLogger("test")
-	logger.setLevel(logging.DEBUG)
-
-	def setUp(self):
-		db = "./tests/data/test_kmeribs-Wheat_Jagger-Flame.tsv" 
-		lens= "./tests/data/test_chr_sizes_jagger.tsv"
-		windows = 200000
-		self.results = IBSpyResults(db, lens, windows)
-
-	def test_windows(self):
-		pd = self.results.count_by_windows()
-
-		self.assertEqual(pd.iloc[0]['variations'], 243)
-		self.assertEqual(pd.iloc[1]['variations'], 180)
-		self.assertEqual(pd.iloc[2]['variations'],  48)
-		self.assertEqual(pd.iloc[3]['variations'], 103)
-		self.assertEqual(pd.iloc[4]['variations'],  66)
-		self.assertEqual(pd.iloc[5]['variations'],  41)
-		self.assertEqual(pd.iloc[6]['variations'],  64)
-
-		self.assertEqual(pd.iloc[0]['seqname'], "chr1A_WhJag")
-		self.assertEqual(pd.iloc[1]['seqname'], "chr1A_WhJag")
-		self.assertEqual(pd.iloc[2]['seqname'], "chr2A_WhJag")
-		self.assertEqual(pd.iloc[3]['seqname'], "chr2A_WhJag")
-		self.assertEqual(pd.iloc[4]['seqname'], "chr2A_WhJag")
-		self.assertEqual(pd.iloc[5]['seqname'], "chr1D_WhJag")
-		self.assertEqual(pd.iloc[6]['seqname'], "chr1D_WhJag")
-
-		self.assertEqual(pd.iloc[0]['start'], 200000)
-		self.assertEqual(pd.iloc[1]['start'], 400000)
-		self.assertEqual(pd.iloc[2]['start'], 400000)
-		self.assertEqual(pd.iloc[3]['start'], 600000)
-		self.assertEqual(pd.iloc[4]['start'], 800000)
-		self.assertEqual(pd.iloc[5]['start'], 800000)
-		self.assertEqual(pd.iloc[6]['start'], 1000000)
+    logger = logging.getLogger("test")
+    logger.setLevel(logging.DEBUG)
+
+    def setUp(self):
+        db = "./tests/data/test_kmeribs-Wheat_Jagger-Flame.tsv" 
+        lens= "./tests/data/test_chr_sizes_jagger.tsv"
+        windows = 200000
+        v_filter = 500
+        self.stitch_number = 3
+        self.n_components = 3
+        self.covariance_type = 'full'
+        self.results = IBSpyResults(db, lens, windows, v_filter)
+
+    def test_windows(self):
+        pd = self.results.count_by_windows()
+
+        self.assertEqual(pd.iloc[0]['variations'], 243)
+        self.assertEqual(pd.iloc[1]['variations'], 180)
+        self.assertEqual(pd.iloc[2]['variations'],  41)
+        self.assertEqual(pd.iloc[3]['variations'], 64)
+        self.assertEqual(pd.iloc[4]['variations'],  48)
+        self.assertEqual(pd.iloc[5]['variations'],  103)
+        self.assertEqual(pd.iloc[6]['variations'],  66)
+
+        self.assertEqual(pd.iloc[0]['seqname'], "chr1A_WhJag")
+        self.assertEqual(pd.iloc[1]['seqname'], "chr1A_WhJag")
+        self.assertEqual(pd.iloc[2]['seqname'], "chr1D_WhJag")
+        self.assertEqual(pd.iloc[3]['seqname'], "chr1D_WhJag")
+        self.assertEqual(pd.iloc[4]['seqname'], "chr2A_WhJag")
+        self.assertEqual(pd.iloc[5]['seqname'], "chr2A_WhJag")
+        self.assertEqual(pd.iloc[6]['seqname'], "chr2A_WhJag")
+
+        self.assertEqual(pd.iloc[0]['window'], 200000)
+        self.assertEqual(pd.iloc[1]['window'], 400000)
+        self.assertEqual(pd.iloc[2]['window'], 800000)
+        self.assertEqual(pd.iloc[3]['window'], 1000000)
+        self.assertEqual(pd.iloc[4]['window'], 400000)
+        self.assertEqual(pd.iloc[5]['window'], 600000)
+        self.assertEqual(pd.iloc[6]['window'], 800000)
+
+    def test_normalize_data(self):
+        log_test, pd = self.results.normalize_data()
+#         print(log_test.tolist())
+        self.assertEqual(log_test[0], np.array([5.493061443340548]))
+        self.assertEqual(log_test[6], np.array([4.189654742026425]))
+
+    def test_stitch_haploBlocks(self):
+        hap_pd = self.results.stitch_haploBlocks(self.n_components,self.covariance_type, self.stitch_number)
+        print(hap_pd)
+        
+        self.assertEqual(hap_pd.iloc[0]['v_mean'], 48)
+        self.assertEqual(hap_pd.iloc[1]['v_mean'], 60)
+        self.assertEqual(hap_pd.iloc[2]['v_mean'], 20)
+        self.assertEqual(hap_pd.iloc[3]['v_mean'], 32)
+        self.assertEqual(hap_pd.iloc[4]['v_mean'], 48)
+        self.assertEqual(hap_pd.iloc[5]['v_mean'], 25)
+        self.assertEqual(hap_pd.iloc[6]['v_mean'], 33)
+
 
 		#print(pd.iloc[0]['variations'])