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

Issue:
Fixes on the table



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index 7b46396..eecd772 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -15,7 +15,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ['3.6', '3.7', '3.8', '3.9']
+        python-version: ['3.7', '3.8', '3.9']
 
     steps:
     - uses: actions/checkout@v2
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
new file mode 100644
index 0000000..8b5ff10
--- /dev/null
+++ b/IBSpy/IBSpy_plot_window_results.py
@@ -0,0 +1,54 @@
+import argparse
+from IBSpy import IBSpyResults, IBSpyPlots
+from matplotlib.backends.backend_pdf import PdfPages
+import matplotlib.pyplot as plt
+import numpy as np
+def parse_arguments():
+    parser = argparse.ArgumentParser()
+
+    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()
+
+    # count variations by windows from IBS_results
+    IBS_results = IBSpyResults(args.IBSpy_counts, args.window_size, args.filter_counts)
+    IBS_results = IBS_results.run_analysis(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/IBSpy_results.py b/IBSpy/IBSpy_results.py
new file mode 100755
index 0000000..90b9ab1
--- /dev/null
+++ b/IBSpy/IBSpy_results.py
@@ -0,0 +1,108 @@
+import pandas as pd
+import numpy as np
+from sklearn.mixture import GaussianMixture
+
+class IBSpyResults:
+    # class variables go here
+    def __init__(self, filename, window_size, filter_counts):
+
+        self.db = pd.read_csv(filename, delimiter='\t')
+        self.window_size = window_size
+        self.filter_counts = filter_counts
+
+    def count_by_windows(self):
+        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)
+
+        # 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, counts):
+        by_windows_db = counts 
+        #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 = varDF['variations'].to_numpy()
+        log_varArray = np.log(varArray, where=varArray>0)
+        log_varArray = log_varArray.reshape((len(log_varArray),1))
+        return log_varArray, varDF
+
+    def build_gmm_model(self, log_varArray, varDF, 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, model, stitch_number):
+        #model = self.build_gmm_model(n_components, covariance_type)
+        gmmDta = model['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)
+        model['vh_block'] = hapBlock.values
+        # put back filtered data with haploblocks
+        hapCntFile = self.count_by_windows()
+        model = pd.merge(hapCntFile, model, left_on=['seqname','window'], right_on=['seqname','window'], how='left')
+        model.loc[:,'v_gmm':'vh_block'] = np.where(model.loc[:, 'v_gmm':'vh_block'] == 1, 1, 0)
+        model.rename(columns={'seqname_x':'seqname', 'variations_x':'variations'}, inplace=True)
+        model = model.drop(['variations_y'], axis=1)
+        return model
+    
+    def run_analysis(self, n_components, covariance_type, stitch_number):
+        counts     = self.count_by_windows()
+        log_test, pd = self.transform_counts_to_log(counts)
+
+        model = self.build_gmm_model(log_test, pd, n_components, covariance_type) 
+        hap_pd = self.stitch_gmm_haplotypes(model, stitch_number)
+        return hap_pd
+        # 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/__init__.py b/IBSpy/__init__.py
index 248e3a8..c9d120f 100644
--- a/IBSpy/__init__.py
+++ b/IBSpy/__init__.py
@@ -2,6 +2,8 @@
 from .kmerGWAS_kmer import KmerGWASDB, KmerGWASDBBuilder
 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/IBSpy/kmerGWAS/kmer_general.c b/IBSpy/kmerGWAS/kmer_general.c
index 2d34cd8..5a391be 100644
--- a/IBSpy/kmerGWAS/kmer_general.c
+++ b/IBSpy/kmerGWAS/kmer_general.c
@@ -52,7 +52,6 @@ kmerGWAS_kmer kmer_shift_and_insert(kmerGWAS_kmer kmer, char base, uint8_t kmer_
 	kmerGWAS_kmer wall  = (kmerGWAS_kmer)  KMERGWAS_MASK << kmer_size * 2;
 	kmer2 |= (kmerGWAS_kmer) n;
 	kmer2 &= ~wall;
-	//fprintf(stderr, "shifting %llu %llu %c %llu %llu\n", kmer, kmer2, base, n, wall);
 	return kmer2 ;
 }
 
@@ -64,10 +63,6 @@ kmerGWAS_kmer kmer2bits(char * kmer, uint8_t kmer_size) {
 		b = kmer_shift_and_insert(b, kmer[i], kmer_size);
 	}
 	b &= KMERGWAS_ORIENTATION_MASK;
-	// for(uint64_t i=0; i < kmer_size; i++) {
-	// 	Nucleotide n = char_to_binary_nucleotide(kmer[kmer_size - i - 1]);
-	// 	b |= ((uint64_t) n << (i*2));
-	// }
 	
 	bt = kmer_reverse_complement(b,  kmer_size);
 	if (bt < b)
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
diff --git a/setup.py b/setup.py
index f055c5d..e5ca3f8 100644
--- a/setup.py
+++ b/setup.py
@@ -22,7 +22,7 @@ def unit_tests():
 
 setuptools.setup(
     name="IBSpy", # Replace with your own username
-    version="0.1.0",
+    version="0.2.0-dev1",
     author="Ricardo H. Ramirez-Gonzalez",
     author_email="ricardo.ramirez-gonzalez@jic.ac.uk",
     description="A package to detect IBS regions",
@@ -44,7 +44,8 @@ def unit_tests():
     entry_points={  # Optional
         'console_scripts': [
             'IBSpy_wcount=IBSpy.IBSpy_window_count:main',
-            'IBSpy=IBSpy:main'
+            'IBSpy=IBSpy:main',
+            'IBSplot=IBSpy.IBSpy_plot_window_results:main'
         ],
     }
 )


Test patch:
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/expected_output/gmm_ibs.tsv.gz b/tests/plots_data/expected_output/gmm_ibs.tsv.gz
new file mode 100644
index 0000000..a93ed26
Binary files /dev/null and b/tests/plots_data/expected_output/gmm_ibs.tsv.gz differ
diff --git a/tests/plots_data/expected_output/gmm_plots.pdf b/tests/plots_data/expected_output/gmm_plots.pdf
new file mode 100644
index 0000000..ce89044
Binary files /dev/null and b/tests/plots_data/expected_output/gmm_plots.pdf differ
diff --git a/tests/plots_data/kmeribs-Wheat_Jagger-Flame.tsv.gz b/tests/plots_data/kmeribs-Wheat_Jagger-Flame.tsv.gz
new file mode 100755
index 0000000..e68454a
Binary files /dev/null and b/tests/plots_data/kmeribs-Wheat_Jagger-Flame.tsv.gz differ
diff --git a/tests/plots_data/test_chr_sizes_jagger.tsv b/tests/plots_data/test_chr_sizes_jagger.tsv
new file mode 100644
index 0000000..62bc739
--- /dev/null
+++ b/tests/plots_data/test_chr_sizes_jagger.tsv
@@ -0,0 +1,4 @@
+chr1A_WhJag	399783
+chr2A_WhJag	699597
+chr1D_WhJag	949442
+
diff --git a/tests/plots_data/test_kmeribs-Wheat_Jagger-Flame.tsv b/tests/plots_data/test_kmeribs-Wheat_Jagger-Flame.tsv
new file mode 100644
index 0000000..887b229
--- /dev/null
+++ b/tests/plots_data/test_kmeribs-Wheat_Jagger-Flame.tsv
@@ -0,0 +1,20 @@
+seqname	start	end	total_kmers	observed_kmers	variations	kmer_distance
+chr1A_WhJag	0	50000	45545	44640	51	736
+chr1A_WhJag	49969	99969	46404	46216	13	195
+chr1A_WhJag	99938	149938	45833	44106	97	1204
+chr1A_WhJag	149907	199907	47777	47315	27	342
+chr1A_WhJag	199876	249876	47765	47000	55	872
+chr1A_WhJag	249845	299845	45625	45115	34	500
+chr1A_WhJag	299814	349814	46707	46033	35	374
+chr1A_WhJag	349783	399783	46518	44008	111	1202
+chr2A_WhJag	399752	449752	48371	47430	48	550
+chr2A_WhJag	449721	499721	48554	48218	23	341
+chr2A_WhJag	499690	549690	49422	49033	28	437
+chr2A_WhJag	549659	599659	49704	49314	31	531
+chr2A_WhJag	599628	649628	49337	49058	21	336
+chr2A_WhJag	649597	699597	49338	48927	27	416
+chr2A_WhJag	699566	749566	49054	48446	39	547
+chr1D_WhJag	749535	799535	48527	48211	22	442
+chr1D_WhJag	799504	849504	49970	49683	19	272
+chr1D_WhJag	849473	899473	49544	49127	29	436
+chr1D_WhJag	899442	949442	49844	49473	35	647
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
new file mode 100644
index 0000000..02ed64d
--- /dev/null
+++ b/tests/test_results.py
@@ -0,0 +1,91 @@
+import os
+import sys
+import numpy as np
+
+
+import logging
+import operator
+import unittest
+from IBSpy import IBSpyResults 
+
+class TestResults(unittest.TestCase):
+    logger = logging.getLogger("test")
+    logger.setLevel(logging.DEBUG)
+
+    def setUp(self):
+        db = "./tests/plots_data/test_kmeribs-Wheat_Jagger-Flame.tsv" 
+        windows = 200000
+        filter_counts = 500
+        self.stitch_number = 3
+        self.n_components = 3
+        self.covariance_type = 'full'
+        self.results = IBSpyResults(db, windows, filter_counts)
+        # print(self.results.count_by_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'],  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_transform_counts_to_log(self):
+        counts = self.results.count_by_windows()
+        log_test, pd = self.results.transform_counts_to_log(counts)
+#         print(log_test.tolist())
+        self.assertEqual( round( log_test[0][0], 2),  round( 5.493061443340548, 2 ) )
+        self.assertEqual( round( log_test[6][0], 2),  round( 4.189654742026425, 2 ) )
+
+    def test_build_model(self):
+        counts = self.results.count_by_windows()
+        log_test, pd = self.results.transform_counts_to_log(counts)
+        model = self.results.build_gmm_model(log_test, pd,self.n_components, self.covariance_type) 
+        print("....\n")
+        print(model)
+        # self.assertEqual(model.iloc[0]['v_gmm'],  0)
+        # self.assertEqual(model.iloc[1]['v_gmm'],  0)
+        # self.assertEqual(model.iloc[2]['v_gmm'],  1)
+        # self.assertEqual(model.iloc[3]['v_gmm'],  1)
+        # self.assertEqual(model.iloc[4]['v_gmm'],  1)
+        # self.assertEqual(model.iloc[5]['v_gmm'],  0)
+        # self.assertEqual(model.iloc[6]['v_gmm'],  0)
+        #TODO: Finish ths test
+
+    def test_stitch_gmm_haplotypes(self):
+        counts = self.results.count_by_windows()
+        log_test, pd = self.results.transform_counts_to_log(counts)
+        model = self.results.build_gmm_model(log_test, pd,self.n_components, self.covariance_type) 
+        hap_pd = self.results.stitch_gmm_haplotypes(model, self.stitch_number)
+        # print(hap_pd)
+        
+        self.assertEqual(hap_pd.iloc[0]['mean'], 48.60)
+        self.assertEqual(hap_pd.iloc[1]['mean'], 60.00)
+        self.assertEqual(hap_pd.iloc[2]['mean'], 20.50)
+        self.assertEqual(hap_pd.iloc[3]['mean'], 32.00)
+        self.assertEqual(hap_pd.iloc[4]['mean'], 48.00)
+        self.assertEqual(hap_pd.iloc[5]['mean'], 25.75)
+        self.assertEqual(hap_pd.iloc[6]['mean'], 33.00)
+
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file