Repository: akikuno/DAJIN2 (language: python)

Issue:
DAJIN2 reports extremely low-frequency alleles  (e.g. 0.005%)
## Describe the Bug

Current DAJIN2 aims to detect alleles greater than 0.5% to differentiate them from Nanopore sequencing errors. In other words, alleles less than or equal to 0.5% are considered sequencing errors.

However, the mutation info and HTML reports show extremely low-frequency alleles, such as 0.005%.

## Expected Behavior

Extremely low-frequency alleles, such as 0.005%, should be classified under the most similar major allele.

## Actual Behavior

None

## Steps/Code to Reproduce

None

## Operating System

WSL2

## Python version

3.10

## DAJIN2 version

0.4.6

## Additional Context

None

An error occurs if there is an "_" in the FASTA header name.
## Describe the bug

If an underscore `_` (e.g., >1_hoge) is included in the header name of the input FASTA file for DAJIN2, the following error occurs:

```python
2024-05-29 17:42:33, INFO, 🏃 Start running DAJIN2 version 0.4.6
2024-05-29 17:42:33, INFO, example_single/control is now processing...
2024-05-29 17:42:33, INFO, Preprocess example_single/control...
2024-05-29 17:42:57, INFO, Output BAM files of example_single/control...
2024-05-29 17:42:57, INFO, 🍵 example_single/control is finished!
2024-05-29 17:42:57, INFO, example_single/sample is now processing...
2024-05-29 17:42:57, INFO, Preprocess example_single/sample...
2024-05-29 17:43:30, INFO, Classify example_single/sample...
2024-05-29 17:43:33, INFO, Clustering example_single/sample...
2024-05-29 17:43:41, INFO, Consensus calling of example_single/sample...
2024-05-29 17:43:41, ERROR, Catch an Exception. Traceback:
Traceback (most recent call last):
  File "/home/kuno/miniconda/envs/env-dajin2/bin/DAJIN2", line 10, in <module>
    sys.exit(execute())
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/main.py", line 236, in execute
    execute_single_mode(arguments)
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/main.py", line 46, in execute_single_mode
    core.execute_sample(arguments)
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/core/core.py", line 187, in execute_sample
    consensus.cache_mutation_loci(ARGS, clust_subset_sample)
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/core/consensus/mutation_extractor.py", line 101, in cache_mutation_loci
    cache_normalized_indels(ARGS, path_midsv_sample)
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/core/consensus/mutation_extractor.py", line 73, in cache_normalized_indels
    sequence = ARGS.fasta_alleles[allele]
KeyError: '1'
```
## Solutions

The cause of the error is the frequent use of `split("_")` on the path without considering the use of underscores in the header name. In DAJIN2, various annotations are added to the header name using `_` as the delimiter. If the user-specified FASTA header name contains `_`, the expected splits are misaligned.

To handle cases where "_" is included, appropriate splitting should be performed. Specifically, it is recommended to remove the FASTA header name before splitting.

The following script contains hard-coded instances of the above issue, which need to be corrected.

```python
DAJIN2/src/DAJIN2/core/consensus/mutation_extractor.py:    allele, label, *_ = path_midsv_sample.stem.split("_")
DAJIN2/src/DAJIN2/core/consensus/mutation_extractor.py:        allele, label, *_ = path_indels_normalized_sample.stem.split("_")
DAJIN2/src/DAJIN2/core/consensus/similarity_searcher.py:    allele, label, *_ = Path(path_midsv_sample).stem.split("_")
DAJIN2/src/DAJIN2/core/preprocess/midsv_caller.py:        preset = path.stem.split("_")[0]
DAJIN2/src/DAJIN2/core/preprocess/midsv_caller.py:        preset = path.stem.split("_")[0]
DAJIN2/src/DAJIN2/core/report/sequence_exporter.py:    allele = header.split("_")[1]
DAJIN2/src/DAJIN2/utils/report_generator.py:        label, allele, type_, *_ = reads["NAME"].split("_")
```


## Steps/Code to Reproduce


## Operating System

WLS2

## Python version

3.10

## DAJIN2 version

0.4.6

## Additional context

Thank you @geedrn for reporting the issue!!

An error occurs if there is an "_" in the FASTA header name.
## Describe the bug

If an underscore `_` (e.g., >1_hoge) is included in the header name of the input FASTA file for DAJIN2, the following error occurs:

```python
2024-05-29 17:42:33, INFO, 🏃 Start running DAJIN2 version 0.4.6
2024-05-29 17:42:33, INFO, example_single/control is now processing...
2024-05-29 17:42:33, INFO, Preprocess example_single/control...
2024-05-29 17:42:57, INFO, Output BAM files of example_single/control...
2024-05-29 17:42:57, INFO, 🍵 example_single/control is finished!
2024-05-29 17:42:57, INFO, example_single/sample is now processing...
2024-05-29 17:42:57, INFO, Preprocess example_single/sample...
2024-05-29 17:43:30, INFO, Classify example_single/sample...
2024-05-29 17:43:33, INFO, Clustering example_single/sample...
2024-05-29 17:43:41, INFO, Consensus calling of example_single/sample...
2024-05-29 17:43:41, ERROR, Catch an Exception. Traceback:
Traceback (most recent call last):
  File "/home/kuno/miniconda/envs/env-dajin2/bin/DAJIN2", line 10, in <module>
    sys.exit(execute())
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/main.py", line 236, in execute
    execute_single_mode(arguments)
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/main.py", line 46, in execute_single_mode
    core.execute_sample(arguments)
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/core/core.py", line 187, in execute_sample
    consensus.cache_mutation_loci(ARGS, clust_subset_sample)
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/core/consensus/mutation_extractor.py", line 101, in cache_mutation_loci
    cache_normalized_indels(ARGS, path_midsv_sample)
  File "/home/kuno/miniconda/envs/env-dajin2/lib/python3.10/site-packages/DAJIN2/core/consensus/mutation_extractor.py", line 73, in cache_normalized_indels
    sequence = ARGS.fasta_alleles[allele]
KeyError: '1'
```
## Solutions

The cause of the error is the frequent use of `split("_")` on the path without considering the use of underscores in the header name. In DAJIN2, various annotations are added to the header name using `_` as the delimiter. If the user-specified FASTA header name contains `_`, the expected splits are misaligned.

To handle cases where "_" is included, appropriate splitting should be performed. Specifically, it is recommended to remove the FASTA header name before splitting.

The following script contains hard-coded instances of the above issue, which need to be corrected.

```python
DAJIN2/src/DAJIN2/core/consensus/mutation_extractor.py:    allele, label, *_ = path_midsv_sample.stem.split("_")
DAJIN2/src/DAJIN2/core/consensus/mutation_extractor.py:        allele, label, *_ = path_indels_normalized_sample.stem.split("_")
DAJIN2/src/DAJIN2/core/consensus/similarity_searcher.py:    allele, label, *_ = Path(path_midsv_sample).stem.split("_")
DAJIN2/src/DAJIN2/core/preprocess/midsv_caller.py:        preset = path.stem.split("_")[0]
DAJIN2/src/DAJIN2/core/preprocess/midsv_caller.py:        preset = path.stem.split("_")[0]
DAJIN2/src/DAJIN2/core/report/sequence_exporter.py:    allele = header.split("_")[1]
DAJIN2/src/DAJIN2/utils/report_generator.py:        label, allele, type_, *_ = reads["NAME"].split("_")
```


## Steps/Code to Reproduce


## Operating System

WLS2

## Python version

3.10

## DAJIN2 version

0.4.6

## Additional context

Thank you @geedrn for reporting the issue!!



Discussion / hints:

Fix patch:
diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
deleted file mode 100644
index 95b2f423..00000000
--- a/.github/ISSUE_TEMPLATE/bug_report.md
+++ /dev/null
@@ -1,44 +0,0 @@
----
-name: 🐛 Bug Report
-about: Create a report to help us improve
-labels: bug
----
-
-## Describe the bug
-
-## Expected Behavior
-
-## Actual Behavior
-
-## Steps/Code to Reproduce
-
-## Operating System
-
-> [!IMPORTANT]
-> DAJIN2 supports only Linux and macOS. If you are using Windows, please use WSL2.
-
-## Python version
-
-You can track the version of python by running the following command in the terminal:
-
-```bash
-python --version
-```
-
-> [!IMPORTANT]
-> DAJIN2 supports Python 3.8 to 3.10, not Python 3.11 yet due to [Bioconda issue] (https://github.com/bioconda/bioconda-recipes/issues/37805).
-
-## DAJIN2 version
-
-You can track the version of DAJIN2 by running the following command in the terminal:
-
-```bash
-DAJIN2 --version
-```
-
-Or, it is recorded in the log file when DAJIN2 is executed.
-
-
-## Additional context
-
-Add any other context or screenshots about the bug here.
diff --git a/.github/ISSUE_TEMPLATE/bug_report.yml b/.github/ISSUE_TEMPLATE/bug_report.yml
new file mode 100644
index 00000000..0f44257c
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/bug_report.yml
@@ -0,0 +1,65 @@
+name: 🐞 Bug Report
+description: Create a report to help us improve
+labels: ['bug']
+
+body:
+  - type: textarea
+    id: description
+    attributes:
+      label: '📋 Description'
+      description: A clear and concise description of the bug.
+    validations:
+      required: true
+
+  - type: textarea
+    attributes:
+      label: Current Behavior
+      description: |
+        Optional: A concise description of what you're experiencing.
+    validations:
+      required: false
+
+  - type: textarea
+    attributes:
+      label: Expected Behavior
+      description: |
+        Optional: A concise description of what you expected to happen.
+    validations:
+      required: false
+
+  - type: textarea
+    attributes:
+      label: Steps To Reproduce
+      description: |
+        Optional: Steps to reproduce the behavior.
+      placeholder: |
+        1. In this environment...
+        1. With this config...
+        1. Run '...'
+        1. See error...
+    validations:
+      required: false
+
+  - type: textarea
+    id: environment
+    attributes:
+      label: '🔍 Environment'
+      description: |
+          Optional: The environment information.
+          - **OS**: WSL (Ubuntu 22.04)
+          - **DAJIN2**: 0.4.6
+      value: |
+          - OS:
+          - DAJIN2:
+      render: markdown
+    validations:
+      required: false
+
+  - type: textarea
+    id: anything_else
+    attributes:
+      label: Anything else?
+      description: |
+        Optional: Add any other contexts, links, or screenshots about the bug here.
+    validations:
+      required: false
diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md
deleted file mode 100644
index 58d911ac..00000000
--- a/.github/ISSUE_TEMPLATE/feature_request.md
+++ /dev/null
@@ -1,25 +0,0 @@
----
-name: 💡 Feature Request
-about: Suggest an idea for this project
-labels: request
----
-
-## Describe the features you want.
-
-
-## Describe the solution you'd like if you have any.
-
-
-## DAJIN2 version
-
-You can track the version of DAJIN2 by running the following command in the terminal:
-
-```bash
-DAJIN2 --version
-```
-
-Or, it is recorded in the log file when DAJIN2 is executed.
-
-## Additional context
-
-Add any other context or screenshots about the feature request here.
diff --git a/.github/ISSUE_TEMPLATE/feature_request.yml b/.github/ISSUE_TEMPLATE/feature_request.yml
new file mode 100644
index 00000000..e70c8b98
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/feature_request.yml
@@ -0,0 +1,37 @@
+name: 💡 Feature Request
+description: Suggest an idea for this project
+labels: ['request']
+
+body:
+  - type: textarea
+    id: description
+    attributes:
+      label: '📋 Description'
+      description: A clear and concise description of the request.
+    validations:
+      required: true
+
+  - type: textarea
+    id: environment
+    attributes:
+      label: '🔍 Environment'
+      description: |
+          Optional: The environment information.
+          - **OS**: WSL (Ubuntu 22.04)
+          - **DAJIN2**: 0.4.6
+      value: |
+          - OS:
+          - DAJIN2:
+      render: markdown
+    validations:
+      required: false
+
+  - type: textarea
+    id: anything_else
+    attributes:
+      label: Anything else?
+      description: |
+        Optional: Add any other contexts, links, or screenshots about the feature request here.
+    validations:
+      required: false
+
diff --git a/.github/workflows/deploy_pypi.yml b/.github/workflows/deploy_pypi.yml
index e038225f..20fe8ea9 100644
--- a/.github/workflows/deploy_pypi.yml
+++ b/.github/workflows/deploy_pypi.yml
@@ -1,28 +1,28 @@
-name: deploy to pypi
+name: Upload Python Package to PyPI
 
 on:
   release:
-    types: [created]
+    types: [published]
 
 jobs:
   deploy:
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v3
-    - name: Set up Python
-      uses: actions/setup-python@v4
-      with:
-        python-version: '3.x'
-    - name: Install dependencies
-      run: |
-        python -m pip install --upgrade pip
-        pip install setuptools wheel twine
-    - name: Build
-      run: |
-        python3 setup.py sdist --formats=zip
-    - name: Publish a Python distribution to PyPI
-      if: startsWith(github.event.ref, 'refs/tags') || github.event_name == 'release'
-      uses: pypa/gh-action-pypi-publish@release/v1
-      with:
-        user: __token__
-        password: ${{ secrets.PYPI_API_TOKEN }}
+      - uses: actions/checkout@v4
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: '3.x'
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install build
+      - name: Build package
+        run: python -m build
+      - name: Publish package
+        uses: pypa/gh-action-pypi-publish@release/v1
+        with:
+          password: ${{ secrets.PYPI_API_TOKEN }}
+
+# Reference:
+# https://docs.github.com/actions/automating-builds-and-tests/building-and-testing-python#publishing-to-package-registries
diff --git a/MANIFEST.in b/MANIFEST.in
deleted file mode 100644
index 7a976c51..00000000
--- a/MANIFEST.in
+++ /dev/null
@@ -1,6 +0,0 @@
-include requirements.txt
-include src/DAJIN2/template_igvjs.html
-
-graft src/DAJIN2/templates
-graft src/DAJIN2/static
-graft src/DAJIN2/utils
diff --git a/README.md b/README.md
index b2c1e280..ddff0255 100644
--- a/README.md
+++ b/README.md
@@ -127,7 +127,7 @@ Options:
 -c, --control             Specify the path to the directory containing control FASTQ files.
 -a, --allele              Specify the path to the FASTA file.
 -n, --name (Optional)     Set the output directory name. Default: 'Results'.
--g, --genome (Optional)   Specify the reference genome ID (e.g., hg38, mm39). Default: '' (empty string).
+-g, --genome (Optional)   Specify the reference UCSC genome ID (e.g., hg38, mm39). Default: '' (empty string).
 -t, --threads (Optional)  Set the number of threads. Default: 1.
 -h, --help                Display this help message and exit.
 -v, --version             Display the version number and exit.
diff --git a/docs/RELEASE.md b/docs/RELEASE.md
index f8e0a6a5..25180051 100644
--- a/docs/RELEASE.md
+++ b/docs/RELEASE.md
@@ -6,26 +6,85 @@
 ## 🐛 Bug Fixes
 ## 🔧 Maintenance
 ## ⛔️ Deprecated
-[Commit Detail](https://github.com/akikuno/DAJIN2/commit/xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx)
+[[Commit Detail](https://github.com/akikuno/DAJIN2/commit/xxxxx)]
 -->
 
 <!-- 💡 ToDo
-- consensus.similarity_searcher.onehot_by_mutationsでは、コントロールアレルに対して、何度も同一のものを生成している可能性がある。
-  - 一箇所にまとめられるようにする
-- リード数が多すぎる場合には、平均的なPhread scoreが高いものを優先して解析する、といったsubsetを行っても良い気が
-
-- deletions_to_fasta.pyを加える
 
 - VCF、PDFを出力する
 - 逆位アレルでの検証を加える
 - nCATSがほしい…
-- Docker imageにして、darwin-arm64でも動くようにする
 - Flaskではなく、streamlitでGUIを作る
 
  -->
 
 <!-- ############################################################# # -->
 
+# v0.5.0 (2024-MM-DD)
+
+## 📝 Documentation
+
++ Update the issue template from md to yml and modify it to make it easier for users to fill out each item.  [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/ee2c1784e3cb0e72fd09b7c7df577082c19c1a88)]
+
+
+## 💥 Breaking
+
++ Extremely low-frequency alleles (less than 0.05%) are considered Nanopore sequence errors and are not clustered #36.
+  + Configure `clustering.extract_labels` so that alleles with a low number of reads (0.05% or fewer or 5 reads or fewer) are not clustered. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/e671b5c84b4cf522faf51823e36fe075b049efcf)]
+  + Change `clustering.clustering` to stop if the minimum value of the elements in the cluster is 0.5% or less. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/74609d5048a4ad8d7004886bf411b2ed4be7fa4b)]
+  + Add `consensus.remove_minor_alleles` to remove minor alleles with fewer than 5 reads or less than 0.5% [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/70f675d6a8ea90e9fca51639ddb2b4609e0f4c80)]
+
+
++ Save subsetted fastq of a control sample if the read number is too large (> 10,000 reads). The control will have a maximum of 10,000 reads to avoid excessive computational load. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/d21827f8bbeec326fa2aa4f28feadd6fdecaf554)]
+
++ If the read length is 500 bases or less, change the mappy preset to `sr`. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/6e56804ad40780e200f4e9c9ea23294b95443aba)]
+
++ Update `extract_best_preset` to prioritize `map-ont` and remove `splice` preset if inversion is observed. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/aa7f6925d6ef4a80a1ba0bbf2b75d8e549ae9863)]
+
+
+Update the algorithms of `cssplits_hander.reallocate_insertion_within_deletion` to automate change point detection by incorporating temporal changes. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/7ed8ac8404d18b86c163c71ded6dd1ba784bce79)]
+
+
+
+## 🔧 Maintenance
+
++ Update `deploy_pypi.yml` to use the latest version of Actions. Refer to [the latest official YAML for guidance](https://docs.github.com/actions/automating-builds-and-tests/building-and-testing-python#publishing-to-package-registries). [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/1a54b40146acd21eee30a3a373c44b419d170ad4)]
+
+
++ Integrate `requirements.txt` and `MANIFEST.in` into `pyproject.toml` by replacing `setup.py` [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/12f255c3a280098f0310755c51e966031c724932)]
+
++ Modify to record the execution command of DAJIN2 in the log file [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/38c97a725f6dd3f00162325bf504142f8f8d6594)]
+
++ Add a test to check if the version in `test_version.sh` matches the version in `pyproject.toml` and `utils.config` [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/a06cb4593ef1a11b3c9826f7ca5532a1bf83f67f)]
+
+
++ Rename `consensus.subset_clust` to `consensus.downsample_by_label` to clarify the function's purpose. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/f6e3f0bc2982996a7dbbc4126a80a7dedd076430)]
+
+
++ Update `extract_unique_insertions` to merge highly similar extracted insertion sequences. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/50fe99f42bcd0bae85bcd0eb4ee371a65f38ea14)]
+  + Fix `extract_unique_insertions`: There is a bug where removing the key twice in fasta_insertions_unique caused the index and key to become misaligned in enumerate(distances) if i != key. Therefore, the removal of keys from fasta_insertions_unique is now done all at once at the end. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/162f248b4deee8c35512b84ec428baec65fd8466)]
+
+
++ Add control characters for `fastx_handler.sanitize_filename` as forbidden chars. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/6b74fce0caa0580c4629a132406206c27a66274d)]
+
+
++ Changed the naming convention for the temporary directory: `<sample_name>/<process_content>/<allele_name>/(<label_name>)/file_name`. Example: `flox/consensus/control/1/mutation_loci.pickle`. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/54fee2f48564c6a29fd5c4151126ba4246e9547c)]
+
++ Move `sanitze_name` function from `utils.fastx_handler`to `utils.io` [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/a78bd5c0ad8f26bafe369da69607faf9a467c039)]
+
+
+## 🐛 Bug Fixes
+
++ Removed `sam_handler.remove_overlapped_reads` to prevent unnecessary trimming of reads. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/a8991edc0620412c384760d0862e34cc4ea6c0f1)]
+
++ Fix `preprocess.insertions_to_fasta.remove_minor_groups` to delete the keys (insertion loci) when insertions are removed and result in an empty dict. This prevents errors when accessing non-existent keys in `subset_insertions`. [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/ae8d887282035552c8fbe5c587e43844d5199952)]
+
++ Fix the bug in `cssplits_handler.convert_cssplits_to_cstag` where the insertion cs tag is not merged with the next cs tag if they have the same operator (e.g., `+A|+A|=T, =T`: before: `+aa=T=T`, after: `+aa=TT`). [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/02d1b4c128004e02671e833136508328d699f53f)]
+
++ Modified the system to separate intermediate files using a directory structure instead of underscores (`_`), ensuring that no errors occur even if users use allele names containing underscores [[Commit Detail](https://github.com/akikuno/DAJIN2/commit/f70948315a114e0c182895ba4320233f26fc1025)]
+  + Thank you @geedrn for reporting the issue #39!
+
+
 
 <!-- ############################################################# # -->
 
@@ -95,8 +154,10 @@
 
 + Fix `preprocess.knockin_handler` to correctly identify the flox knock-in sites as deletions not present in the control.  [Commit Detail](https://github.com/akikuno/DAJIN2/commit/d4d267c99f8c51d3a3f88f67882bead66685f710)
 
-+ Bug fix of `reallocate_insertion_within_deletion` [Commit Detail](https://github.com/akikuno/DAJIN2/commit/7cafabf3a89b75af73f1a40da16e6390d47b48d5)
-  - In the script that considers the region between two deletions as an insertion sequence, the size of the other deletion was not taken into account. Even if there was a single base deletion, the entire sequence between the deletions was considered as an insertion sequence. Therefore, the region between two deletions is now defined only if the size of both deletions is equal to or greater than the specified threshold (default = 3).
++ Bug fix and update `reallocate_insertion_within_deletion` [Commit Detail](https://github.com/akikuno/DAJIN2/commit/2f356546999f645a8cb8d33a1fc2f64bc6742113)
+  - In the script that considers the region between two deletions as an insertion sequence, the size of the other deletion was not taken into account. Even if there was a single base deletion, the entire sequence between the deletions was considered as an insertion sequence. 
+  - Therefore, the region between two deletions is now defined as (1) identifying bins where deletions are enriched within appropriate bins (500 bp) continuously, and (2) extracting the precise break points from the start and end of these bins, implementing an algorithm to extract the large deletion region.
+
 
 </details>
 
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 00000000..1afb3ad4
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,53 @@
+[build-system]
+requires = ["poetry-core>=1.0.0"]
+build-backend = "poetry.core.masonry.api"
+
+[tool.poetry]
+name = "DAJIN2"
+version = "0.5.0"
+description = "One-step genotyping tools for targeted long-read sequencing"
+authors = ["Akihiro Kuno <akuno@md.tsukuba.ac.jp>"]
+readme = "README.md"
+homepage = "https://github.com/akikuno/DAJIN2"
+license = "MIT"
+packages = [
+    { include = "DAJIN2", from = "src" }
+]
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "License :: OSI Approved :: MIT License",
+    "Operating System :: POSIX",
+    "Operating System :: MacOS",
+    "Intended Audience :: Science/Research",
+    "Topic :: Scientific/Engineering :: Bio-Informatics",
+]
+include = [
+    "src/DAJIN2/template_igvjs.html", 
+    "src/DAJIN2/templates/**", 
+    "src/DAJIN2/static/**", 
+    "src/DAJIN2/utils/**",
+]
+
+[tool.poetry.dependencies]
+python = "^3.8"
+numpy = ">=1.24.0"
+scipy = ">=1.10.0"
+pandas = ">=1.0.0"
+openpyxl = ">=3.1.0"
+rapidfuzz = ">=3.6.0"
+scikit-learn = ">=1.3.0"
+mappy = ">=2.24"
+pysam = ">=0.21.0"
+Flask = ">=2.2.0"
+waitress = ">=2.1.0"
+Jinja2 = ">=3.1.0"
+plotly = ">=5.19.0"
+kaleido = ">=0.2.0"
+cstag = ">=1.0.0"
+midsv = ">=0.11.0"
+wslPath = ">=0.4.1"
+ruptures = ">=1.1.8"
+
+[tool.poetry.scripts]
+DAJIN2 = "DAJIN2.main:execute"
+
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 34726a6a..00000000
--- a/requirements.txt
+++ /dev/null
@@ -1,22 +0,0 @@
-certifi
-
-numpy >= 1.24.0
-scipy >= 1.10.0
-pandas >= 1.0.0
-openpyxl >= 3.1.0
-rapidfuzz >=3.6.0
-scikit-learn >= 1.3.0
-
-mappy >= 2.24
-pysam >= 0.21.0
-
-Flask >= 2.2.0
-waitress >= 2.1.0
-Jinja2 >= 3.1.0
-
-plotly >= 5.19.0
-kaleido >= 0.2.0
-
-cstag >= 1.0.0
-midsv >= 0.11.0
-wslPath >=0.4.1
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 7c9175a0..00000000
--- a/setup.py
+++ /dev/null
@@ -1,35 +0,0 @@
-import setuptools
-from pathlib import Path
-
-this_directory = Path(__file__).parent
-long_description = (this_directory / "README.md").read_text()
-
-with open("requirements.txt") as requirements_file:
-    install_requirements = requirements_file.read().splitlines()
-
-setuptools.setup(
-    name="DAJIN2",
-    version="0.4.6",
-    author="Akihiro Kuno",
-    author_email="akuno@md.tsukuba.ac.jp",
-    description="One-step genotyping tools for targeted long-read sequencing",
-    long_description=long_description,
-    long_description_content_type="text/markdown",
-    url="https://github.com/akikuno/DAJIN2",
-    install_requires=install_requirements,
-    packages=setuptools.find_packages(
-        where="src",
-    ),
-    package_dir={"": "src"},
-    entry_points={"console_scripts": ["DAJIN2=DAJIN2.main:execute"]},
-    include_package_data=True,
-    classifiers=[
-        "Environment :: Console",
-        "Programming Language :: Python :: 3",
-        "License :: OSI Approved :: MIT License",
-        "Operating System :: POSIX",
-        "Operating System :: MacOS",
-        "Intended Audience :: Science/Research",
-        "Topic :: Scientific/Engineering :: Bio-Informatics",
-    ],
-)
diff --git a/src/DAJIN2/core/classification/classifier.py b/src/DAJIN2/core/classification/classifier.py
index 52a78331..a3074d17 100644
--- a/src/DAJIN2/core/classification/classifier.py
+++ b/src/DAJIN2/core/classification/classifier.py
@@ -45,7 +45,7 @@ def extract_alleles_with_max_score(score_of_each_alleles: list[dict]) -> list[di
 def classify_alleles(TEMPDIR: Path, FASTA_ALLELES: dict, SAMPLE_NAME: str) -> list[dict]:
     score_of_each_alleles = []
     for allele in FASTA_ALLELES:
-        path_midsv = Path(TEMPDIR, SAMPLE_NAME, "midsv", f"{allele}.json")
+        path_midsv = Path(TEMPDIR, SAMPLE_NAME, "midsv", allele, f"{SAMPLE_NAME}.jsonl")
         score_of_each_alleles.extend(score_allele(path_midsv, allele))
 
     score_of_each_alleles_merged = merge_minor_alleles(score_of_each_alleles)
diff --git a/src/DAJIN2/core/clustering/clustering.py b/src/DAJIN2/core/clustering/clustering.py
index 3a9f1c73..c50646c2 100644
--- a/src/DAJIN2/core/clustering/clustering.py
+++ b/src/DAJIN2/core/clustering/clustering.py
@@ -23,13 +23,14 @@
 
 
 def count_number_of_clusters(labels_control: list[int], coverage_control: int) -> int:
-    """Reads < 1% in the control are considered clustering errors and are not counted"""
-    return sum(1 for reads in Counter(labels_control).values() if reads / coverage_control > 0.01)
+    """If there is less than 1% lead within a cluster, they are considered clustering errors and are not counted"""
+    return sum(1 for control_reads in Counter(labels_control).values() if control_reads / coverage_control > 0.01)
 
 
-def optimize_labels(X: spmatrix, coverage_sample: int, coverage_control: int) -> list[int]:
-    labels_previous = list(range(coverage_sample))
-    for i in range(1, coverage_sample):
+def optimize_labels(X: spmatrix, coverage_sample: int, coverage_control: int, min_cluster_size: int) -> list[int]:
+    labels_previous = [1 for _ in range(coverage_sample)]
+
+    for i in range(2, coverage_sample):
         np.random.seed(seed=1)
         labels_all = BisectingKMeans(n_clusters=i, random_state=1).fit_predict(X).tolist()
 
@@ -45,11 +46,14 @@ def optimize_labels(X: spmatrix, coverage_sample: int, coverage_control: int) ->
         Return the number of clusters when:
             - the number of clusters in control is split into more than one.
             - the mutual information between the current and previous labels is high enough (= similar).
+            - the minimum cluster size is reached.
         To reduce the allele number, previous labels are returned.
         """
-        if num_labels_control >= 2 or rand_index >= 0.95:
+        if num_labels_control >= 2 or rand_index >= 0.95 or min_cluster_size > min(Counter(labels_current).values()):
             return labels_previous
+
         labels_previous = labels_current
+
     return labels_previous
 
 
@@ -58,21 +62,31 @@ def get_label_most_common(labels: list[int]) -> int:
 
 
 def return_labels(
-    path_score_sample: Path, path_score_control: Path, path_sample: Path, strand_bias_in_control: bool
+    path_score_sample: Path,
+    path_score_control: Path,
+    path_sample: Path,
+    strand_bias_in_control: bool,
+    min_cluster_size: int,
 ) -> list[int]:
     np.random.seed(seed=1)
+
     score_control = list(io.read_jsonl(path_score_control))
     X_control = csr_matrix(score_control)
+
     """Subset to 1000 reads of controls in the most common cluster to remove outliers and reduce computation time"""
     labels_control = BisectingKMeans(n_clusters=2, random_state=1).fit_predict(X_control)
     label_most_common = get_label_most_common(labels_control)
     scores_control_subset = subset_scores(labels_control, io.read_jsonl(path_score_control), label_most_common, 1000)
+
     scores_sample = list(io.read_jsonl(path_score_sample))
     X = csr_matrix(list(chain(scores_sample, scores_control_subset)))
+
     coverage_sample = io.count_newlines(path_score_sample)
     coverage_control = len(scores_control_subset)
-    labels = optimize_labels(X, coverage_sample, coverage_control)
+    labels = optimize_labels(X, coverage_sample, coverage_control, min_cluster_size)
+
     """Re-allocate clusters with strand bias to clusters without strand bias"""
     if strand_bias_in_control is False:
         labels = remove_biased_clusters(path_sample, path_score_sample, labels)
+
     return labels
diff --git a/src/DAJIN2/core/clustering/label_extractor.py b/src/DAJIN2/core/clustering/label_extractor.py
index 5b8435d8..d74482fc 100644
--- a/src/DAJIN2/core/clustering/label_extractor.py
+++ b/src/DAJIN2/core/clustering/label_extractor.py
@@ -15,31 +15,36 @@
 def extract_labels(classif_sample, TEMPDIR, SAMPLE_NAME, CONTROL_NAME) -> list[dict[str]]:
     labels_result = []
     max_label = 1
-    strand_bias = is_strand_bias(Path(TEMPDIR, CONTROL_NAME, "midsv", "control.json"))
+
+    strand_bias = is_strand_bias(Path(TEMPDIR, CONTROL_NAME, "midsv", "control", f"{CONTROL_NAME}.jsonl"))
     classif_sample.sort(key=lambda x: x["ALLELE"])
+    min_samples = max(5, int(len(classif_sample) * 0.5 // 100))  # 0.5% of the samples
+
     for allele, group in groupby(classif_sample, key=lambda x: x["ALLELE"]):
         """Cache data to temporary files"""
-        if Path(TEMPDIR, CONTROL_NAME, "midsv", f"{allele}.json").exists():
-            path_control = Path(TEMPDIR, CONTROL_NAME, "midsv", f"{allele}.json")
-        else:
-            path_control = Path(TEMPDIR, CONTROL_NAME, "midsv", f"{allele}_{SAMPLE_NAME}.json")
+
+        # Insertion/Inversion allele
+        path_control = Path(TEMPDIR, CONTROL_NAME, "midsv", allele, f"{SAMPLE_NAME}.jsonl")
+        if not path_control.exists():
+            path_control = Path(TEMPDIR, CONTROL_NAME, "midsv", allele, f"{CONTROL_NAME}.jsonl")
 
         unique_id = str(uuid.uuid4())
         path_sample = Path(TEMPDIR, SAMPLE_NAME, "clustering", f"tmp_{allele}_{unique_id}.jsonl")
         io.write_jsonl(data=group, file_path=path_sample)
 
         """Load mutation_loci and knockin_loci."""
-        path_mutation_loci = Path(TEMPDIR, SAMPLE_NAME, "mutation_loci", f"{allele}.pickle")
+        path_mutation_loci = Path(TEMPDIR, SAMPLE_NAME, "mutation_loci", allele, "mutation_loci.pickle")
+        path_knockin_loci = Path(TEMPDIR, SAMPLE_NAME, "knockin_loci", allele, "knockin.pickle")
+
         mutation_loci: list[set[str]] = io.load_pickle(path_mutation_loci)
-        path_knockin_loci = Path(TEMPDIR, SAMPLE_NAME, "knockin_loci", f"{allele}.pickle")
         knockin_loci: set[int] = io.load_pickle(path_knockin_loci) if path_knockin_loci.exists() else set()
 
         """Skip clustering when the number of reads is too small or there is no mutation."""
         read_numbers = io.count_newlines(path_sample)
         is_no_mutation = all(m == set() for m in mutation_loci)
-        if read_numbers <= 5 or is_no_mutation:
-            max_label += 1
+        if read_numbers < max(5, min_samples) or is_no_mutation:
             labels_result.extend([max_label] * read_numbers)
+            max_label += 1
             continue
 
         """Calculate scores to temporary files."""
@@ -54,7 +59,7 @@ def extract_labels(classif_sample, TEMPDIR, SAMPLE_NAME, CONTROL_NAME) -> list[d
         io.write_jsonl(data=scores_control, file_path=path_score_control)
 
         """Extract labels."""
-        labels = return_labels(path_score_sample, path_score_control, path_sample, strand_bias)
+        labels = return_labels(path_score_sample, path_score_control, path_sample, strand_bias, min_samples)
         labels_result += relabel_with_consective_order(labels, start=max_label)
 
         """Remove temporary files."""
diff --git a/src/DAJIN2/core/consensus/__init__.py b/src/DAJIN2/core/consensus/__init__.py
index a22db208..5a8496ac 100644
--- a/src/DAJIN2/core/consensus/__init__.py
+++ b/src/DAJIN2/core/consensus/__init__.py
@@ -2,5 +2,5 @@
 from DAJIN2.core.consensus.name_handler import call_allele_name
 from DAJIN2.core.consensus.name_handler import update_key_by_allele_name
 from DAJIN2.core.consensus.name_handler import add_key_by_allele_name
-from DAJIN2.core.consensus.clust_formatter import subset_clust
+from DAJIN2.core.consensus.clust_formatter import remove_minor_alleles, downsample_by_label
 from DAJIN2.core.consensus.mutation_extractor import cache_mutation_loci
diff --git a/src/DAJIN2/core/consensus/clust_formatter.py b/src/DAJIN2/core/consensus/clust_formatter.py
index 8d2e9cd0..5c7c619f 100644
--- a/src/DAJIN2/core/consensus/clust_formatter.py
+++ b/src/DAJIN2/core/consensus/clust_formatter.py
@@ -1,11 +1,31 @@
 from __future__ import annotations
 
+import random
 from itertools import groupby
+from collections import defaultdict
 
 
-def subset_clust(clust_sample: list[dict], num: int = 1000) -> list[dict]:
-    clust_subset_sample = []
+def remove_minor_alleles(clust_sample: list[dict]) -> list[dict]:
+    """Remove minor alleles from the clustering sample."""
+    min_sample_size = max(5, int(len(clust_sample) * 0.5 // 100))
+
+    counts = defaultdict(int)
+    for c in clust_sample:
+        counts[c["LABEL"]] += 1
+
+    minor_labels = {label for label, count in counts.items() if count < min_sample_size}
+
+    return [c for c in clust_sample if c["LABEL"] not in minor_labels]
+
+
+def downsample_by_label(clust_sample: list[dict], num: int = 1000) -> list[dict]:
+    random.seed(1)
+
+    clust_downsampled = []
     clust_sample.sort(key=lambda x: x["LABEL"])
     for _, group in groupby(clust_sample, key=lambda x: x["LABEL"]):
-        clust_subset_sample.extend(list(group)[:num])
-    return clust_subset_sample
+        samples = list(group)
+        num_samples = min(len(samples), num)
+        clust_downsampled.extend(random.sample(samples, num_samples))
+
+    return clust_downsampled
diff --git a/src/DAJIN2/core/consensus/consensus.py b/src/DAJIN2/core/consensus/consensus.py
index ddb8e0ed..28eaba0f 100644
--- a/src/DAJIN2/core/consensus/consensus.py
+++ b/src/DAJIN2/core/consensus/consensus.py
@@ -98,7 +98,6 @@ class ConsensusKey:
 
 
 def call_consensus(tempdir: Path, sample_name: str, clust_sample: list[dict]) -> tuple[dict[list], dict[str]]:
-    path_consensus = Path(tempdir, sample_name, "consensus")
 
     clust_sample.sort(key=lambda x: [x["ALLELE"], x["LABEL"]])
 
@@ -108,10 +107,11 @@ def call_consensus(tempdir: Path, sample_name: str, clust_sample: list[dict]) ->
     for (allele, label), group in groupby(clust_sample, key=lambda x: [x["ALLELE"], x["LABEL"]]):
         clust = list(group)
 
-        mutation_loci = io.load_pickle(Path(path_consensus, f"{allele}_{label}_mutation_loci.pickle"))
+        path_consensus = Path(tempdir, sample_name, "consensus", allele, str(label))
+        cons_mutation_loci = io.load_pickle(Path(path_consensus, "mutation_loci.pickle"))
 
         cssplits = [cs["CSSPLIT"].split(",") for cs in clust]
-        cons_percentage = call_percentage(cssplits, mutation_loci)
+        cons_percentage = call_percentage(cssplits, cons_mutation_loci)
 
         key = ConsensusKey(allele, label, clust[0]["PERCENT"])
         cons_percentages[key] = cons_percentage
diff --git a/src/DAJIN2/core/consensus/mutation_extractor.py b/src/DAJIN2/core/consensus/mutation_extractor.py
index 44303e04..c81d9adf 100644
--- a/src/DAJIN2/core/consensus/mutation_extractor.py
+++ b/src/DAJIN2/core/consensus/mutation_extractor.py
@@ -36,11 +36,11 @@ def get_thresholds(path_indels_normalized_sample, path_indels_normalized_control
     return thresholds
 
 
-def extract_path_control(ARGS, allele: str) -> Path:
+def extract_path_control(tempdir: Path, control_name: str, sample_name: str, allele: str) -> Path:
     # Define potential file paths
     paths = [
-        Path(ARGS.tempdir, ARGS.control_name, "midsv", f"{allele}.json"),
-        Path(ARGS.tempdir, ARGS.control_name, "midsv", f"{allele}_{ARGS.sample_name}.json"),
+        Path(tempdir, control_name, "midsv", allele, f"{control_name}.jsonl"),
+        Path(tempdir, control_name, "midsv", allele, f"{sample_name}.jsonl"),
     ]
     # Return the first path that exists
     for path in paths:
@@ -48,64 +48,66 @@ def extract_path_control(ARGS, allele: str) -> Path:
             return path
 
 
-def extract_path_n_filtered_control(ARGS, path_midsv_control: Path) -> Path:
+def extract_path_n_filtered_control(
+    tempdir: Path, control_name: str, sample_name: str, path_control: Path, allele: str
+) -> Path:
     """
     Filter out reads with N enriched in the control.
     """
-    path_output = Path(ARGS.tempdir, ARGS.control_name, "consensus", f"{path_midsv_control.stem}_n_filtered.jsonl")
+    path_output = Path(tempdir, control_name, "consensus", allele, f"{sample_name}_n_filtered.jsonl")
     if path_output.exists():
         return path_output
 
-    midsv_control = io.read_jsonl(path_midsv_control)
+    path_output.parent.mkdir(parents=True, exist_ok=True)
+
+    midsv_control = io.read_jsonl(path_control)
     n_counts = np.array([sum(1 if cs == "N" else 0 for cs in c["CSSPLIT"].split(",")) for c in midsv_control])
 
     kmeans = MiniBatchKMeans(n_clusters=2, random_state=0).fit(n_counts.reshape(-1, 1))
     threshold = kmeans.cluster_centers_.mean()
     labels = np.where(n_counts <= threshold, True, False)
-    midsv_filtered_control = (m for m, label in zip(io.read_jsonl(path_midsv_control), labels) if label)
+    midsv_filtered_control = (m for m, label in zip(io.read_jsonl(path_control), labels) if label)
     io.write_jsonl(midsv_filtered_control, path_output)
 
     return path_output
 
 
-def cache_normalized_indels(ARGS, path_midsv_sample: Path) -> None:
-    allele, label, *_ = path_midsv_sample.stem.split("_")
+def cache_normalized_indels(ARGS, path_consensus_sample: Path, allele: str) -> None:
     sequence = ARGS.fasta_alleles[allele]
 
-    path_midsv_control = extract_path_control(ARGS, allele)
-    path_midsv_n_filtered_control = extract_path_n_filtered_control(ARGS, path_midsv_control)
-
-    cache_selected_control_by_similarity(
-        ARGS, path_midsv_n_filtered_control, path_midsv_sample, path_midsv_sample.parent
+    path_control = extract_path_control(ARGS.tempdir, ARGS.control_name, ARGS.sample_name, allele)
+    path_control_filtered = extract_path_n_filtered_control(
+        ARGS.tempdir, ARGS.control_name, ARGS.sample_name, path_control, allele
     )
 
-    path_midsv_similar_control = Path(path_midsv_sample.parent, f"{allele}_{label}_control.jsonl")
+    cache_selected_control_by_similarity(ARGS, path_control_filtered, path_consensus_sample, allele)
 
-    _, indels_normalized_sample = summarize_indels(path_midsv_sample, sequence)
+    path_midsv_similar_control = Path(path_consensus_sample.parent, "control.jsonl")
+
+    _, indels_normalized_sample = summarize_indels(path_consensus_sample, sequence)
     _, indels_normalized_control = summarize_indels(path_midsv_similar_control, sequence)
 
-    path_consensus = Path(ARGS.tempdir, ARGS.sample_name, "consensus")
-    io.save_pickle(indels_normalized_sample, Path(path_consensus, f"{allele}_{label}_normalized_sample.pickle"))
-    io.save_pickle(indels_normalized_control, Path(path_consensus, f"{allele}_{label}_normalized_control.pickle"))
+    io.save_pickle(indels_normalized_sample, Path(path_consensus_sample.parent, "normalized_sample.pickle"))
+    io.save_pickle(indels_normalized_control, Path(path_consensus_sample.parent, "normalized_control.pickle"))
 
 
 def cache_mutation_loci(ARGS, clust_sample: list[dict]) -> None:
     # Separate clusters by label and cache them
     clust_sample.sort(key=lambda x: [x["ALLELE"], x["LABEL"]])
-    path_consensus = Path(ARGS.tempdir, ARGS.sample_name, "consensus")
     for (allele, label), group in groupby(clust_sample, key=lambda x: [x["ALLELE"], x["LABEL"]]):
-        io.write_jsonl(group, Path(path_consensus, f"{allele}_{label}_sample.jsonl"))
+        path_consensus = Path(ARGS.tempdir, ARGS.sample_name, "consensus", allele, str(label))
+        path_consensus.mkdir(parents=True, exist_ok=True)
+
+        path_consensus_sample = Path(path_consensus, f"{ARGS.sample_name}_sample.jsonl")
+        io.write_jsonl(group, path_consensus_sample)
 
-    # Cache normalized indels counts
-    for path_midsv_sample in path_consensus.glob("*_sample.jsonl"):
-        cache_normalized_indels(ARGS, path_midsv_sample)
+        cache_normalized_indels(ARGS, path_consensus_sample, allele)
 
-    # Extract and cache mutation loci
-    for path_indels_normalized_sample in path_consensus.glob("*_normalized_sample.pickle"):
-        allele, label, *_ = path_indels_normalized_sample.stem.split("_")
-        path_indels_normalized_control = Path(path_consensus, f"{allele}_{label}_normalized_control.pickle")
+        # Extract and cache mutation loci
+        path_indels_normalized_sample = Path(path_consensus, "normalized_sample.pickle")
+        path_indels_normalized_control = Path(path_consensus, "normalized_control.pickle")
         sequence = ARGS.fasta_alleles[allele]
-        path_knockin = Path(ARGS.tempdir, ARGS.sample_name, "knockin_loci", f"{allele}.pickle")
+        path_knockin = Path(ARGS.tempdir, ARGS.sample_name, "knockin_loci", allele, "knockin.pickle")
 
         thresholds = get_thresholds(path_indels_normalized_sample, path_indels_normalized_control)
 
@@ -120,4 +122,4 @@ def cache_mutation_loci(ARGS, clust_sample: list[dict]) -> None:
             is_consensus,
         )
 
-        io.save_pickle(mutation_loci, Path(path_consensus, f"{allele}_{label}_mutation_loci.pickle"))
+        io.save_pickle(mutation_loci, Path(path_consensus, "mutation_loci.pickle"))
diff --git a/src/DAJIN2/core/consensus/similarity_searcher.py b/src/DAJIN2/core/consensus/similarity_searcher.py
index 0649a21a..698cbf63 100644
--- a/src/DAJIN2/core/consensus/similarity_searcher.py
+++ b/src/DAJIN2/core/consensus/similarity_searcher.py
@@ -72,21 +72,21 @@ def identify_normal_reads(
 ###########################################################
 
 
-def filter_control(ARGS, path_midsv_control: Path, path_midsv_sample: Path) -> list[bool]:
+def filter_control(ARGS, path_consensus_control: Path, path_consensus_sample: Path, allele: str) -> list[bool]:
     """
     find similar control reads compared to sample reads
     """
-    cssplits = (m["CSSPLIT"].split(",") for m in io.read_jsonl(path_midsv_sample))
+    cssplits = (m["CSSPLIT"].split(",") for m in io.read_jsonl(path_consensus_sample))
     coverage_match = np.array([sum(1 for cs in cssplit if cs.startswith("=")) for cssplit in zip(*cssplits)])
-    mut_onehot_sample = onehot_by_mutations(io.read_jsonl(path_midsv_sample))
+    mut_onehot_sample = onehot_by_mutations(io.read_jsonl(path_consensus_sample))
 
     path_mut_onehot_control = Path(
-        ARGS.tempdir, ARGS.control_name, "consensus", f"{path_midsv_control.stem}_onehot.pickle"
+        ARGS.tempdir, ARGS.control_name, "consensus", allele, f"{ARGS.sample_name}_onehot.pickle"
     )
     if path_mut_onehot_control.exists():
         mut_onehot_control = io.load_pickle(path_mut_onehot_control)
     else:
-        mut_onehot_control = onehot_by_mutations(io.read_jsonl(path_midsv_control))
+        mut_onehot_control = onehot_by_mutations(io.read_jsonl(path_consensus_control))
         io.save_pickle(mut_onehot_control, path_mut_onehot_control)
 
     mut_percentage_sample = calculate_percentage(mut_onehot_sample, coverage_match)
@@ -99,11 +99,10 @@ def filter_control(ARGS, path_midsv_control: Path, path_midsv_sample: Path) -> l
 
 
 def cache_selected_control_by_similarity(
-    ARGS, path_midsv_control: Path, path_midsv_sample: Path, path_output: Path
+    ARGS, path_consensus_control: Path, path_consensus_sample: Path, allele: str
 ) -> None:
-    normal_reads_flags = filter_control(ARGS, path_midsv_control, path_midsv_sample)
-    midsv_control = io.read_jsonl(path_midsv_control)
+    normal_reads_flags = filter_control(ARGS, path_consensus_control, path_consensus_sample, allele)
+    midsv_control = io.read_jsonl(path_consensus_control)
     midsv_filtered = (m for m, flag in zip(midsv_control, normal_reads_flags) if flag is True)
 
-    allele, label, *_ = Path(path_midsv_sample).stem.split("_")
-    io.write_jsonl(midsv_filtered, Path(path_output, f"{allele}_{label}_control.jsonl"))
+    io.write_jsonl(midsv_filtered, Path(path_consensus_sample.parent, "control.jsonl"))
diff --git a/src/DAJIN2/core/core.py b/src/DAJIN2/core/core.py
index e9f650ac..5f6dead0 100644
--- a/src/DAJIN2/core/core.py
+++ b/src/DAJIN2/core/core.py
@@ -5,7 +5,7 @@
 
 from pathlib import Path
 
-from DAJIN2.utils import io, config, fastx_handler
+from DAJIN2.utils import io, fastx_handler
 from DAJIN2.core import classification, clustering, consensus, preprocess, report
 from DAJIN2.core.preprocess.input_formatter import FormattedInputs
 
@@ -19,8 +19,6 @@
 
 def execute_control(arguments: dict):
 
-    logger.info(f"\N{runner} Start running DAJIN2 version {config.DAJIN_VERSION}")
-
     logger.info(f"{arguments['control']} is now processing...")
 
     ###########################################################
@@ -36,21 +34,24 @@ def execute_control(arguments: dict):
     ###########################################################
     if Path(ARGS.tempdir, "report", "BAM", ARGS.control_name, f"{ARGS.control_name}.bam").exists():
         logger.info(f"{arguments['control']} is already preprocessed and reuse the results for the current run...")
-        return
+        return None
+
     logger.info(f"Preprocess {arguments['control']}...")
 
     ###########################################################
     # Merge fastq files
     ###########################################################
-    fastx_handler.save_concatenated_fastx(ARGS.tempdir, ARGS.path_control)
 
+    fastx_handler.save_concatenated_fastx(ARGS.tempdir, ARGS.path_control, ARGS.control_name)
+    # Save subsetted fastq if the read number is too large (> 10,000 reads)
+    fastx_handler.save_subsetted_fastx(
+        Path(ARGS.tempdir, ARGS.control_name, "fastq", f"{ARGS.control_name}.fastq.gz"), num_reads=10_000
+    )
     ###########################################################
     # Mapping
     ###########################################################
 
-    # ============================================================
     # Export fasta files as single-FASTA format
-    # ============================================================
     fastx_handler.export_fasta_files(ARGS.tempdir, ARGS.fasta_alleles, ARGS.control_name)
 
     # ============================================================
@@ -99,7 +100,7 @@ def execute_sample(arguments: dict):
     # Merge fastq files
     # ============================================================
 
-    fastx_handler.save_concatenated_fastx(ARGS.tempdir, ARGS.path_sample)
+    fastx_handler.save_concatenated_fastx(ARGS.tempdir, ARGS.path_sample, ARGS.sample_name)
 
     # ============================================================
     # Mapping with mappy
@@ -114,6 +115,7 @@ def execute_sample(arguments: dict):
     # ============================================================
     # MIDSV conversion
     # ============================================================
+
     preprocess.generate_midsv(ARGS, is_control=False, is_insertion=False)
 
     # ============================================================
@@ -181,18 +183,21 @@ def execute_sample(arguments: dict):
 
     logger.info(f"Consensus calling of {arguments['sample']}...")
 
+    # Remove minor alleles with fewer than 5 reads or less than 0.5%
+    clust_sample_removed = consensus.remove_minor_alleles(clust_sample)
+
     # Downsampling to 1000 reads in each LABEL
-    clust_subset_sample = consensus.subset_clust(clust_sample, 1000)
+    clust_downsampled = consensus.downsample_by_label(clust_sample_removed, 1000)
 
-    consensus.cache_mutation_loci(ARGS, clust_subset_sample)
+    consensus.cache_mutation_loci(ARGS, clust_downsampled)
 
-    cons_percentage, cons_sequence = consensus.call_consensus(ARGS.tempdir, ARGS.sample_name, clust_subset_sample)
+    cons_percentage, cons_sequence = consensus.call_consensus(ARGS.tempdir, ARGS.sample_name, clust_downsampled)
 
     allele_names = consensus.call_allele_name(cons_sequence, cons_percentage, ARGS.fasta_alleles)
     cons_percentage = consensus.update_key_by_allele_name(cons_percentage, allele_names)
     cons_sequence = consensus.update_key_by_allele_name(cons_sequence, allele_names)
 
-    RESULT_SAMPLE = consensus.add_key_by_allele_name(clust_sample, allele_names)
+    RESULT_SAMPLE = consensus.add_key_by_allele_name(clust_sample_removed, allele_names)
     RESULT_SAMPLE.sort(key=lambda x: x["LABEL"])
 
     io.save_pickle(cons_percentage, Path(ARGS.tempdir, ARGS.sample_name, "consensus", "cons_percentage.pickle"))
diff --git a/src/DAJIN2/core/preprocess/input_formatter.py b/src/DAJIN2/core/preprocess/input_formatter.py
index d2260f5f..12885085 100644
--- a/src/DAJIN2/core/preprocess/input_formatter.py
+++ b/src/DAJIN2/core/preprocess/input_formatter.py
@@ -19,9 +19,9 @@ def parse_arguments(arguments: dict) -> tuple:
         )
 
     return (
-        arguments["sample"],
-        arguments["control"],
-        arguments["allele"],
+        Path(arguments["sample"]),
+        Path(arguments["control"]),
+        Path(arguments["allele"]),
         arguments["name"],
         arguments["threads"],
         genome_urls,
@@ -73,9 +73,9 @@ def get_genome_coordinates(genome_urls: dict, fasta_alleles: dict, is_cache_geno
 
 @dataclass(frozen=True)
 class FormattedInputs:
-    path_sample: str
-    path_control: str
-    path_allele: str
+    path_sample: Path
+    path_control: Path
+    path_allele: Path
     sample_name: str
     control_name: str
     fasta_alleles: dict[str, str]
@@ -99,8 +99,8 @@ def format_inputs(arguments: dict) -> FormattedInputs:
         path_sample,
         path_control,
         path_allele,
-        sample_name,
-        control_name,
+        io.sanitize_name(sample_name),
+        io.sanitize_name(control_name),
         fasta_alleles,
         tempdir,
         genome_coordinates,
diff --git a/src/DAJIN2/core/preprocess/insertions_to_fasta.py b/src/DAJIN2/core/preprocess/insertions_to_fasta.py
index aa2ffdf4..dbe95c03 100644
--- a/src/DAJIN2/core/preprocess/insertions_to_fasta.py
+++ b/src/DAJIN2/core/preprocess/insertions_to_fasta.py
@@ -32,16 +32,14 @@ def remove_non_alphabets(cssplits: str) -> str:
 ###########################################################
 
 
-def clustering_insertions(cssplits_insertion: list[str], n_decoy: int = 1000) -> list[int]:
+def clustering_insertions(cssplits_insertion: list[str]) -> list[int]:
     seq_all = [remove_non_alphabets(seq) for seq in cssplits_insertion]
     query = seq_all[0]
     _, distances, _ = zip(*process.extract_iter(query, seq_all, scorer=DamerauLevenshtein.normalized_distance))
 
-    # By adding upper (1) and lower (0) limits, we prevent errors where minor differences are clustered (e.g., 0.1 and 0.2 becoming separate clusters).
-
     insertion_lengths = [[len(c) for c in cs.split(",")] for cs in cssplits_insertion]
 
-    scores = [s + [d] for s, d in zip(insertion_lengths, distances)]
+    scores = [length + [distance] for length, distance in zip(insertion_lengths, distances)]
 
     return MeanShift(bin_seeding=True).fit_predict(np.array(scores)).tolist()
 
@@ -164,12 +162,14 @@ def get_merged_insertion(insertion: dict[str, int], labels: np.ndarray) -> dict[
     return insertion_merged
 
 
-def remove_minor_groups(insertions_merged: dict[tuple[int], dict[tuple[str], int]], coverage: int, threshold: float = 0.5) -> dict[tuple[int], dict[tuple[str], int]]:
+def remove_minor_groups(insertions_merged: dict[tuple[int], dict[tuple[str], int]], coverage: int, percentage: float = 0.5) -> dict[tuple[int], dict[tuple[str], int]]:
+    """Remove minor groups with less than {percentage} % coverage or less than 5 reads."""
+    threshold = max(5, int(coverage * percentage // 100))
     for _, ins in insertions_merged.items():
         # Create a list of elements to delete
         to_delete = []
         for seq, count in ins.items():
-            if count < coverage * threshold // 100:
+            if count < threshold:
                 to_delete.append(seq)
 
         # Delete the collected elements
@@ -180,7 +180,7 @@ def remove_minor_groups(insertions_merged: dict[tuple[int], dict[tuple[str], int
 
 
 def merge_similar_insertions(
-    insertions: dict[int, dict[str, int]], mutation_loci: list[set[str]], coverage: int, threshold: float = 0.5
+    insertions: dict[int, dict[str, int]], mutation_loci: list[set[str]], coverage: int, percentage: float = 0.5
 ) -> dict[tuple[int], dict[tuple[str], int]]:
     index_grouped = group_index_by_consecutive_insertions(mutation_loci)
     insertions_grouped = group_insertions(insertions, index_grouped)
@@ -190,10 +190,11 @@ def merge_similar_insertions(
             seq, count = next(iter(insertion.items()))
             insertions_merged[idx] = {tuple([seq]): count}
             continue
-        labels = clustering_insertions(insertion, n_decoy=1000)
+
+        labels = clustering_insertions(insertion)
         insertions_merged[idx] = get_merged_insertion(insertion, labels)
 
-    return remove_minor_groups(insertions_merged, coverage, threshold)
+    return remove_minor_groups(insertions_merged, coverage, percentage)
 
 
 ###########################################################
@@ -294,10 +295,11 @@ def filter_minor_label(
     labels: list[int],
     insertions_scores_sequences: list[tuple[list[int], str]],
     coverage: int,
-    threshold: float = 0.5,
+    percentage: float = 0.5,
 ) -> tuple[list[int], list[str]]:
+    threshold = max(5, int(coverage * percentage // 100))
     labels_, counts_ = np.unique(labels, return_counts=True)
-    minor_labels = {label for label, count in zip(labels_, counts_) if count < coverage * threshold // 100}
+    minor_labels = {label for label, count in zip(labels_, counts_) if count < threshold}
     index_minor_labels = {i for i, label in enumerate(labels) if label in minor_labels}
     labels_filtered = [label for i, label in enumerate(labels) if i not in index_minor_labels]
     score_seq_filterd = [
@@ -449,6 +451,7 @@ def extract_index_of_insertions(
                 max_idx = idx
         if max_idx != -1:
             index_of_insertions.append(max_idx)
+
     return index_of_insertions
 
 
@@ -464,6 +467,7 @@ def generate_cstag(
                 continue
             list_sequence[idx] = convert_cssplits_to_cstag([seq]) + "="
         cstag_insertions[label] = "cs:Z:=" + "".join(list_sequence)
+
     return cstag_insertions
 
 
@@ -471,21 +475,33 @@ def generate_fasta(cstag_insertions: dict[str, str]) -> dict[str, str]:
     fasta_insertions = dict()
     for label, cs_tag in cstag_insertions.items():
         fasta_insertions[label] = cstag.to_sequence(cs_tag)
+
     return fasta_insertions
 
 
-def extract_unique_insertions(FASTA_ALLELES: dict[str, str], fasta_insertions: dict[str, str]) -> dict[str, str]:
+def extract_unique_insertions(fasta_insertions: dict[str, str], FASTA_ALLELES: dict[str, str]) -> dict[str, str]:
     """
     Extract unique insertion alleles if they are dissimilar to the FASTA_ALLELES input by the user.
     "Unique insertion alleles" are defined as sequences that have a difference of more than 10 bases compared to the sequences in FASTA_ALLELES
     """
-    to_keep = []
-    for query_key, query_seq in fasta_insertions.items():
-        _, distances, _ = zip(*process.extract_iter(query_seq, FASTA_ALLELES.values(), scorer=DamerauLevenshtein.distance))
-        if all(d > 10 for d in distances):
-            to_keep.append(query_key)
+    fasta_insertions_unique = fasta_insertions.copy()
+
+    # Remove insertion alleles that are similar to the FASTA_ALLELES input by the user
+    to_delete = set()
+    for key, seq in fasta_insertions_unique.items():
+        _, distances, _ = zip(*process.extract_iter(seq, FASTA_ALLELES.values(), scorer=DamerauLevenshtein.distance))
+        if any(d < 10 for d in distances):
+            to_delete |= {key}
+
+    # Remove insertion alleles that are similar to each other
+    for key, seq in fasta_insertions_unique.items():
+        if key in to_delete:
+            continue
+        _, distances, _ = zip(*process.extract_iter(seq, fasta_insertions_unique.values(), scorer=DamerauLevenshtein.distance))
+        similar_index = {i if d < 10 else None for i, d in enumerate(distances) if i != key}
+        to_delete |= similar_index
 
-    return {key: fasta_insertions[key] for key in to_keep if key in fasta_insertions}
+    return {k: v for k, v in fasta_insertions_unique.items() if k not in to_delete}
 
 
 def update_labels(d: dict, FASTA_ALLELES: dict) -> dict:
@@ -504,6 +520,7 @@ def update_labels(d: dict, FASTA_ALLELES: dict) -> dict:
         if user_defined_alleles.isdisjoint(set(d_updated)):
             break
         digits_up += 1
+
     return d_updated
 
 
@@ -532,9 +549,9 @@ def remove_temporal_files(TEMPDIR: Path, SAMPLE_NAME: str) -> None:
 
 
 def generate_insertion_fasta(TEMPDIR, SAMPLE_NAME, CONTROL_NAME, FASTA_ALLELES) -> None:
-    PATH_SAMPLE = Path(TEMPDIR, SAMPLE_NAME, "midsv", "control.json")
-    PATH_CONTROL = Path(TEMPDIR, CONTROL_NAME, "midsv", "control.json")
-    MUTATION_LOCI = io.load_pickle(Path(TEMPDIR, SAMPLE_NAME, "mutation_loci", "control.pickle"))
+    PATH_SAMPLE = Path(TEMPDIR, SAMPLE_NAME, "midsv", "control", f"{SAMPLE_NAME}.jsonl")
+    PATH_CONTROL = Path(TEMPDIR, CONTROL_NAME, "midsv", "control", f"{CONTROL_NAME}.jsonl")
+    MUTATION_LOCI = io.load_pickle(Path(TEMPDIR, SAMPLE_NAME, "mutation_loci", "control", "mutation_loci.pickle"))
 
     insertions = extract_insertions(PATH_SAMPLE, PATH_CONTROL, MUTATION_LOCI)
     if insertions == dict():
@@ -543,7 +560,7 @@ def generate_insertion_fasta(TEMPDIR, SAMPLE_NAME, CONTROL_NAME, FASTA_ALLELES)
 
     coverage: int = io.count_newlines(PATH_SAMPLE)
 
-    insertions_merged = merge_similar_insertions(insertions, MUTATION_LOCI, coverage, threshold=0.5)
+    insertions_merged = merge_similar_insertions(insertions, MUTATION_LOCI, coverage, percentage=0.5)
     if all(True if v == {} else False for v in insertions_merged.values()):
         """"If all the insertion alleles were minor, return None"""
         return None
@@ -551,9 +568,9 @@ def generate_insertion_fasta(TEMPDIR, SAMPLE_NAME, CONTROL_NAME, FASTA_ALLELES)
     # Clustering similar insertion alleles
     insertions_scores_sequences = extract_score_and_sequence(PATH_SAMPLE, insertions_merged)
     cssplits_insertion = [cssplit for _, cssplit in insertions_scores_sequences]
-    labels = clustering_insertions(cssplits_insertion, n_decoy=1000)
+    labels = clustering_insertions(cssplits_insertion)
     labels_filtered, insertion_scores_sequences_filtered = filter_minor_label(
-        labels, insertions_scores_sequences, coverage, threshold=0.5
+        labels, insertions_scores_sequences, coverage, percentage=0.5
     )
 
     # Consensus calling
@@ -573,7 +590,7 @@ def generate_insertion_fasta(TEMPDIR, SAMPLE_NAME, CONTROL_NAME, FASTA_ALLELES)
     index_of_insertions = extract_index_of_insertions(insertions, insertions_merged)
     cstag_insertions = generate_cstag(consensus_of_insertions, index_of_insertions, FASTA_ALLELES["control"])
     fasta_insertions = generate_fasta(cstag_insertions)
-    fasta_insertions_unique = extract_unique_insertions(FASTA_ALLELES, fasta_insertions)
+    fasta_insertions_unique = extract_unique_insertions(fasta_insertions, FASTA_ALLELES)
 
     if fasta_insertions_unique == dict():
         remove_temporal_files(TEMPDIR, SAMPLE_NAME)
diff --git a/src/DAJIN2/core/preprocess/knockin_handler.py b/src/DAJIN2/core/preprocess/knockin_handler.py
index 13b08a2c..d3b86a36 100644
--- a/src/DAJIN2/core/preprocess/knockin_handler.py
+++ b/src/DAJIN2/core/preprocess/knockin_handler.py
@@ -35,7 +35,8 @@ def extract_knockin_loci(TEMPDIR: str | Path, SAMPLE_NAME: str) -> None:
             continue
 
         allele = knockin_fasta.stem
-        path_output = Path(TEMPDIR, SAMPLE_NAME, "knockin_loci", f"{allele}.pickle")
+        path_output = Path(TEMPDIR, SAMPLE_NAME, "knockin_loci", allele, "knockin.pickle")
+        path_output.parent.mkdir(parents=True, exist_ok=True)
         if path_output.exists():
             continue
 
diff --git a/src/DAJIN2/core/preprocess/mapping.py b/src/DAJIN2/core/preprocess/mapping.py
index 40c070d5..c74b9828 100644
--- a/src/DAJIN2/core/preprocess/mapping.py
+++ b/src/DAJIN2/core/preprocess/mapping.py
@@ -6,7 +6,7 @@
 from pathlib import Path
 from typing import Generator
 
-from DAJIN2.utils.dna_handler import revcomp
+from DAJIN2.utils import dna_handler
 
 
 def to_sam(
@@ -55,7 +55,7 @@ def to_sam(
 
             # Handle reverse complement for negative strand
             if hit.strand == -1:
-                query_seq = revcomp(query_seq)
+                query_seq = dna_handler.revcomp(query_seq)
                 if query_qual:
                     query_qual = query_qual[::-1]
 
@@ -109,17 +109,26 @@ def generate_sam(
         path_fastq = Path(ARGS.tempdir, ARGS.sample_name, "fastq", f"{ARGS.sample_name}.fastq.gz")
         name = ARGS.sample_name
 
-    out_directory = Path(ARGS.tempdir, name, "sam")
-
     for path_fasta in paths_fasta:
         name_fasta = Path(path_fasta).stem
-        for preset in ["map-ont", "splice"]:
+        len_sequence = len(Path(path_fasta).read_text().split("\n")[1])
+        if len_sequence < 500:
+            presets = ["sr"]
+        else:
+            presets = ["map-ont", "splice"]
+
+        path_sam_directory = Path(ARGS.tempdir, name, "sam", name_fasta)
+        path_sam_directory.mkdir(parents=True, exist_ok=True)
+
+        for preset in presets:
             sam = to_sam(path_fasta, path_fastq, preset=preset, threads=ARGS.threads, options=mappy_options)
+
             if is_control and is_insertion:
-                path_sam = Path(out_directory, f"{preset}_{name_fasta}_{ARGS.sample_name}.sam")
+                path_sam_file = Path(path_sam_directory, f"{ARGS.sample_name}_{preset}.sam")
             else:
-                path_sam = Path(out_directory, f"{preset}_{name_fasta}.sam")
-            path_sam.write_text("\n".join(sam))
+                path_sam_file = Path(path_sam_directory, f"{preset}.sam")
+
+            path_sam_file.write_text("\n".join(sam))
 
 
 ########################################################################
diff --git a/src/DAJIN2/core/preprocess/midsv_caller.py b/src/DAJIN2/core/preprocess/midsv_caller.py
index 60b282a2..e0998ea2 100644
--- a/src/DAJIN2/core/preprocess/midsv_caller.py
+++ b/src/DAJIN2/core/preprocess/midsv_caller.py
@@ -4,9 +4,7 @@
 
 from pathlib import Path
 from typing import Generator
-from itertools import chain, groupby
-
-from collections import Counter
+from collections import Counter, defaultdict
 
 from DAJIN2.utils import io, sam_handler, cssplits_handler
 
@@ -24,51 +22,69 @@ def has_inversion_in_splice(CIGAR: str) -> bool:
     return False
 
 
-def extract_qname_of_map_ont(sam_ont: Generator[list[str]], sam_splice: Generator[list[str]]) -> set[str]:
-    """Extract qname of reads from `map-ont` when:
-    - no inversion signal in `splice` alignment (insertion + deletion)
-    - single read
-    - long alignment length
-    """
-    alignments_splice = {s[0]: s for s in sam_splice if not s[0].startswith("@")}
-    alignments_ont = sorted([s for s in sam_ont if not s[0].startswith("@")], key=lambda x: x[0])
-    qname_of_map_ont = set()
-    for qname_ont, group in groupby(alignments_ont, key=lambda x: x[0]):
-        group = list(group)
-
-        if qname_ont not in alignments_splice:
-            qname_of_map_ont.add(qname_ont)
+def extract_preset_and_cigar_by_qname(path_sam_files: list[Path]) -> dict[dict[str, str]]:
+    preset_cigar_by_qname = defaultdict(dict)
+    # Extract preset and CIGAR
+    for path in path_sam_files:
+        preset = path.stem
+        sam: list[list[str]] = io.read_sam(path)
+        for record in sam:
+            if record[0].startswith("@"):
+                continue
+            qname = record[0]
+            cigar = record[5]
+            if qname in preset_cigar_by_qname and preset in preset_cigar_by_qname[qname]:
+                cigar += preset_cigar_by_qname[qname][preset]
+            preset_cigar_by_qname[qname].update({preset: cigar})
+    return dict(preset_cigar_by_qname)
+
+
+def extract_best_preset(preset_cigar_by_qname: dict[str, dict[str, str]]) -> dict[str, str]:
+    best_preset = defaultdict(str)
+    for qname in preset_cigar_by_qname:
+        preset_cigar = preset_cigar_by_qname[qname]
+        # If there is an inversion in the splice preset, remove it
+        if "splice" in preset_cigar and has_inversion_in_splice(preset_cigar["splice"]):
+            preset_cigar.pop("splice")
+
+        if preset_cigar == {}:
             continue
 
-        cigar_splice = alignments_splice[qname_ont][5]
+        alignment_lengths = {
+            preset: sam_handler.calculate_alignment_length(cigar) for preset, cigar in preset_cigar.items()
+        }
 
-        # If preset=splice and inversion is present, `midsv.transform`` will not work, so use preset=map-ont.
-        if has_inversion_in_splice(cigar_splice):
-            qname_of_map_ont.add(qname_ont)
+        # If all alignment lengths are the same, prioritize map-ont
+        if len(set(alignment_lengths.values())) == 1 and "map-ont" in alignment_lengths:
+            best_preset[qname] = "map-ont"
             continue
 
-        # only accept single read or inversion reads
-        # TODO: 逆位のリードは単純に「３つのリードがある」という条件でいいのか？
-        if len(group) == 2 or len(group) >= 4:
-            continue
-        cigar_ont = group[0][5]
+        # Define a custom key function to prioritize map-ont
+        def custom_key(key: str) -> tuple[int, bool]:
+            return (alignment_lengths[key], key == "map-ont")
 
-        alignment_length_ont = sam_handler.calculate_alignment_length(cigar_ont)
-        alignment_length_splice = sam_handler.calculate_alignment_length(cigar_splice)
+        max_key = max(alignment_lengths, key=custom_key)
+        best_preset[qname] = max_key
 
-        if alignment_length_ont >= alignment_length_splice:
-            qname_of_map_ont.add(qname_ont)
-    return qname_of_map_ont
+    return dict(best_preset)
 
 
-def filter_sam_by_preset(sam: Generator[list[str]], qname_of_map_ont: set, preset: str = "map-ont") -> Generator:
-    for alignment in sam:
-        if alignment[0].startswith("@"):
-            yield alignment
-        elif preset == "map-ont" and alignment[0] in qname_of_map_ont:
-            yield alignment
-        elif preset != "map-ont" and alignment[0] not in qname_of_map_ont:
-            yield alignment
+def extract_best_alignment_length_from_sam(
+    path_sam_files: list[Path], best_preset: dict[str, str]
+) -> Generator[list[str]]:
+    flag_header = False
+    for path in path_sam_files:
+        preset = path.stem
+        sam = io.read_sam(path)
+        for record in sam:
+            if record[0].startswith("@"):
+                if not flag_header:
+                    yield record
+            else:
+                qname = record[0]
+                if best_preset.get(qname) == preset:
+                    yield record
+        flag_header = True
 
 
 def transform_to_midsv_format(sam: Generator[list[str]]) -> Generator[list[dict]]:
@@ -127,11 +143,6 @@ def filter_samples_by_n_proportion(midsv_sample: Generator[dict], threshold: int
 # convert_consecutive_indels_to_match
 ###########################################################
 
-"""
-Due to alignment errors, there can be instances where a true match is mistakenly replaced with "insertion following a deletion".
-For example, although it should be "=C,=T", it gets replaced by "-C,+C|=T". In such cases, a process is performed to revert it back to "=C,=T".
-"""
-
 
 def convert_consecutive_indels_to_match(cssplit: str) -> str:
     i = 0
@@ -176,17 +187,15 @@ def convert_consecutive_indels_to_match(cssplit: str) -> str:
 
 
 def convert_consecutive_indels(midsv_sample: Generator) -> Generator[list[dict]]:
+    """
+    Due to alignment errors, there can be instances where a true match is mistakenly replaced with "insertion following a deletion".
+    For example, although it should be "=C,=T", it gets replaced by "-C,+C|=T". In such cases, a process is performed to revert it back to "=C,=T".
+    """
     for m in midsv_sample:
         m["CSSPLIT"] = convert_consecutive_indels_to_match(m["CSSPLIT"])
         yield m
 
 
-def reallocate_insertions_within_deletion(midsv_sample: Generator) -> Generator[list[dict]]:
-    for m in midsv_sample:
-        m["CSSPLIT"] = cssplits_handler.reallocate_insertion_within_deletion(m["CSSPLIT"])
-        yield m
-
-
 ###########################################################
 # main
 ###########################################################
@@ -196,33 +205,31 @@ def generate_midsv(ARGS, is_control: bool = False, is_insertion: bool = False) -
     name = ARGS.control_name if is_control else ARGS.sample_name
 
     for allele, sequence in ARGS.fasta_alleles.items():
-        if Path(ARGS.tempdir, name, "midsv", f"{allele}.json").exists():
+        path_midsv_directory = Path(ARGS.tempdir, name, "midsv", allele)
+        path_midsv_directory.mkdir(parents=True, exist_ok=True)
+
+        if Path(path_midsv_directory, f"{name}.jsonl").exists():
             continue
 
         if is_control and is_insertion:
             """
             Set the destination for midsv as `barcode01/midsv/insertion1_barcode02.json` when control is barcode01, sample is barcode02, and the allele is insertion1.
             """
-            path_ont = Path(ARGS.tempdir, name, "sam", f"map-ont_{allele}_{ARGS.sample_name}.sam")
-            path_splice = Path(ARGS.tempdir, name, "sam", f"splice_{allele}_{ARGS.sample_name}.sam")
-            path_output_midsv = Path(ARGS.tempdir, name, "midsv", f"{allele}_{ARGS.sample_name}.json")
+            path_sam_files = list(Path(ARGS.tempdir, name, "sam", allele).glob(f"{ARGS.sample_name}_*.sam"))
+            path_midsv_output = Path(ARGS.tempdir, name, "midsv", allele, f"{ARGS.sample_name}.jsonl")
         else:
             """
             Set the destination for midsv as `barcode02/midsv/insertion1.json` when the sample is barcode02 and the allele is insertion1.
             """
-            path_ont = Path(ARGS.tempdir, name, "sam", f"map-ont_{allele}.sam")
-            path_splice = Path(ARGS.tempdir, name, "sam", f"splice_{allele}.sam")
-            path_output_midsv = Path(ARGS.tempdir, name, "midsv", f"{allele}.json")
-
-        sam_ont = sam_handler.remove_overlapped_reads(list(io.read_sam(path_ont)))
-        sam_splice = sam_handler.remove_overlapped_reads(list(io.read_sam(path_splice)))
-        qname_of_map_ont = extract_qname_of_map_ont(sam_ont, sam_splice)
-        sam_of_map_ont = filter_sam_by_preset(sam_ont, qname_of_map_ont, preset="map-ont")
-        sam_of_splice = filter_sam_by_preset(sam_splice, qname_of_map_ont, preset="splice")
-        midsv_chaind = transform_to_midsv_format(chain(sam_of_map_ont, sam_of_splice))
+            path_sam_files = list(Path(ARGS.tempdir, name, "sam", allele).glob(f"*.sam"))
+            path_midsv_output = Path(ARGS.tempdir, name, "midsv", allele, f"{name}.jsonl")
+
+        preset_cigar_by_qname = extract_preset_and_cigar_by_qname(path_sam_files)
+        best_preset = extract_best_preset(preset_cigar_by_qname)
+        sam_best_alignments = extract_best_alignment_length_from_sam(path_sam_files, best_preset)
+        midsv_chaind = transform_to_midsv_format(sam_best_alignments)
         midsv_sample = replace_internal_n_to_d(midsv_chaind, sequence)
         midsv_sample = convert_flag_to_strand(midsv_sample)
         midsv_sample = filter_samples_by_n_proportion(midsv_sample)
         midsv_sample = convert_consecutive_indels(midsv_sample)
-        midsv_sample = reallocate_insertions_within_deletion(midsv_sample)
-        midsv.write_jsonl(midsv_sample, path_output_midsv)
+        midsv.write_jsonl(midsv_sample, path_midsv_output)
diff --git a/src/DAJIN2/core/preprocess/mutation_extractor.py b/src/DAJIN2/core/preprocess/mutation_extractor.py
index 2549182f..cb5800d1 100644
--- a/src/DAJIN2/core/preprocess/mutation_extractor.py
+++ b/src/DAJIN2/core/preprocess/mutation_extractor.py
@@ -261,19 +261,24 @@ def transpose_mutation_loci(mutation_loci: dict[str, set[int]], sequence: str) -
 
 
 def cache_indels_count(ARGS, is_control: bool = False) -> None:
+    dirname = ARGS.control_name if is_control else ARGS.sample_name
     for allele, sequence in ARGS.fasta_alleles.items():
-        if is_control and Path(ARGS.tempdir, ARGS.control_name, "midsv", f"{allele}_{ARGS.sample_name}.json").exists():
-            prefix = f"{allele}_{ARGS.sample_name}"
+        path_mutation_loci = Path(ARGS.tempdir, dirname, "mutation_loci", allele)
+        path_mutation_loci.mkdir(parents=True, exist_ok=True)
+
+        if not is_control:
+            prefix = ARGS.sample_name
         else:
-            prefix = allele
+            path_insertion = Path(ARGS.tempdir, ARGS.control_name, "midsv", allele, f"{ARGS.sample_name}.jsonl")
+            if path_insertion.exists():
+                prefix = ARGS.sample_name
+            else:
+                prefix = ARGS.control_name
 
-        path_mutation_loci = Path(ARGS.tempdir, ARGS.control_name if is_control else ARGS.sample_name, "mutation_loci")
         if Path(path_mutation_loci, f"{prefix}_count.pickle").exists():
             continue
 
-        path_midsv = Path(
-            ARGS.tempdir, ARGS.control_name if is_control else ARGS.sample_name, "midsv", f"{prefix}.json"
-        )
+        path_midsv = Path(ARGS.tempdir, dirname, "midsv", allele, f"{prefix}.jsonl")
         indels_count, indels_normalized = summarize_indels(path_midsv, sequence)
         io.save_pickle(indels_count, Path(path_mutation_loci, f"{prefix}_count.pickle"))
         io.save_pickle(indels_normalized, Path(path_mutation_loci, f"{prefix}_normalized.pickle"))
@@ -318,24 +323,25 @@ def cache_mutation_loci(ARGS, is_control: bool = False) -> None:
     if is_control:
         return
 
-    path_mutation_sample = Path(ARGS.tempdir, ARGS.sample_name, "mutation_loci")
-    path_mutation_control = Path(ARGS.tempdir, ARGS.control_name, "mutation_loci")
-
     for allele, sequence in ARGS.fasta_alleles.items():
-        path_output = Path(path_mutation_sample, f"{allele}.pickle")
-        if path_output.exists():
+        path_mutation_sample = Path(ARGS.tempdir, ARGS.sample_name, "mutation_loci", allele)
+        path_mutation_control = Path(ARGS.tempdir, ARGS.control_name, "mutation_loci", allele)
+
+        path_output_mutation_loci = Path(path_mutation_sample, "mutation_loci.pickle")
+        if path_output_mutation_loci.exists():
             continue
 
-        file_name = f"{allele}_{ARGS.sample_name}_normalized.pickle"
+        file_name = f"{ARGS.sample_name}_normalized.pickle"
         if not Path(path_mutation_control, file_name).exists():
-            file_name = f"{allele}_normalized.pickle"
+            file_name = f"{ARGS.control_name}_normalized.pickle"
+
         path_indels_normalized_control = Path(path_mutation_control, file_name)
 
-        path_indels_normalized_sample = Path(path_mutation_sample, f"{allele}_normalized.pickle")
-        path_knockin = Path(ARGS.tempdir, ARGS.sample_name, "knockin_loci", f"{allele}.pickle")
+        path_indels_normalized_sample = Path(path_mutation_sample, f"{ARGS.sample_name}_normalized.pickle")
+        path_knockin = Path(ARGS.tempdir, ARGS.sample_name, "knockin_loci", allele, "knockin.pickle")
 
         mutation_loci: list[set[str]] = extract_mutation_loci(
             sequence, path_indels_normalized_sample, path_indels_normalized_control, path_knockin
         )
 
-        io.save_pickle(mutation_loci, path_output)
+        io.save_pickle(mutation_loci, path_output_mutation_loci)
diff --git a/src/DAJIN2/core/report/bam_exporter.py b/src/DAJIN2/core/report/bam_exporter.py
index e6b86fc0..0d0eae38 100644
--- a/src/DAJIN2/core/report/bam_exporter.py
+++ b/src/DAJIN2/core/report/bam_exporter.py
@@ -100,7 +100,6 @@ def write_sam_to_bam(sam: list[list[str]], path_sam: str | Path, path_bam: str |
 
 def update_sam(sam: list, GENOME_COODINATES: dict = {}) -> list:
     sam_records = sam.copy()
-    sam_records = sam_handler.remove_overlapped_reads(sam_records)
     sam_records = sam_handler.remove_microhomology(sam_records)
     if GENOME_COODINATES["genome"]:
         return recalculate_sam_coodinates_to_reference(sam_records, GENOME_COODINATES)
@@ -109,7 +108,10 @@ def update_sam(sam: list, GENOME_COODINATES: dict = {}) -> list:
 
 
 def export_to_bam(TEMPDIR, NAME, GENOME_COODINATES, THREADS, UUID, RESULT_SAMPLE=None, is_control=False) -> None:
-    path_sam_input = Path(TEMPDIR, NAME, "sam", "map-ont_control.sam")
+    path_sam_input = Path(TEMPDIR, NAME, "sam", "control", "map-ont.sam")
+    if not path_sam_input.exists():  # In the case of short-read.
+        path_sam_input = Path(TEMPDIR, NAME, "sam", "control", "sr.sam")
+
     sam_records = list(io.read_sam(path_sam_input))
 
     # Update sam
diff --git a/src/DAJIN2/core/report/sequence_exporter.py b/src/DAJIN2/core/report/sequence_exporter.py
index 3661c5b5..a0518b5b 100644
--- a/src/DAJIN2/core/report/sequence_exporter.py
+++ b/src/DAJIN2/core/report/sequence_exporter.py
@@ -1,9 +1,10 @@
 from __future__ import annotations
 
-import textwrap
 import cstag
+import textwrap
 from pathlib import Path
-from DAJIN2.utils.cssplits_handler import convert_cssplits_to_cstag
+
+from DAJIN2.utils.cssplits_handler import convert_cssplits_to_cstag, reallocate_insertion_within_deletion
 from DAJIN2.core.report.insertion_reflector import reflect_ref_insertion_to_query
 
 
@@ -16,8 +17,9 @@ def convert_to_fasta(header: str, sequence: str) -> str:
 
 
 def convert_to_html(TEMPDIR: Path, SAMPLE_NAME: str, header: str, cons_per: list[dict]) -> str:
-    cons_cssplit = [max(cons, key=cons.get) for cons in cons_per]
-    cons_cstag = convert_cssplits_to_cstag(cons_cssplit)
+    cons_cssplits = [max(cons, key=cons.get) for cons in cons_per]
+    cons_cssplits_reallocated = reallocate_insertion_within_deletion(cons_cssplits, bin_size=500, percentage=50)
+    cons_cstag = convert_cssplits_to_cstag(cons_cssplits_reallocated)
 
     allele = header.split("_")[1]
     if Path(TEMPDIR, SAMPLE_NAME, "cstag", f"{allele}.txt").exists():
diff --git a/src/DAJIN2/main.py b/src/DAJIN2/main.py
index 95f1bb85..8ae4080e 100644
--- a/src/DAJIN2/main.py
+++ b/src/DAJIN2/main.py
@@ -11,6 +11,7 @@
 
 import sys
 import shutil
+import logging
 import argparse
 from pathlib import Path
 from copy import deepcopy
@@ -21,12 +22,9 @@
 from DAJIN2.utils import io, config, report_generator, input_validator, multiprocess
 
 
-def generate_report(name: str) -> None:
+def generate_report(name: str, logger: logging.Logger) -> None:
     report_generator.report(name)
-    print(
-        f"\N{party popper} Finished! Open {config.DAJIN_RESULTS_DIR}/{name} to see the report.",
-        file=sys.stderr,
-    )
+    logger.info(f"\N{party popper} Finished! Open {config.DAJIN_RESULTS_DIR}/{name} to see the report.")
 
 
 ################################################################################
@@ -35,18 +33,26 @@ def generate_report(name: str) -> None:
 
 
 def execute_single_mode(arguments: dict[str]):
-    # set logging to export log to stderr and file
+    # Set logging to export log to stderr and file
     path_logfile = config.get_logfile()
-    config.set_logging(path_logfile)
+    logger = config.set_logging(path_logfile)
+    logger.info(f"\N{runner} Start running DAJIN2 version {config.DAJIN_VERSION}")
+    logger.info(f"\N{Personal Computer} {' '.join(sys.argv)}")
+
+    # Validate input files
     input_validator.validate_files(arguments["sample"], arguments["control"], arguments["allele"])
     if arguments.get("genome"):
         arguments.update(input_validator.validate_genome_and_fetch_urls(arguments["genome"]))
-    # run DAJIN2
+
+    # Run DAJIN2
     core.execute_control(arguments)
     core.execute_sample(arguments)
-    # finish
-    generate_report(arguments["name"])
+
+    # Finish call
+    generate_report(arguments["name"], logger)
     shutil.move(path_logfile, Path("DAJIN_Results", arguments["name"]))
+
+    # Remove temporary files
     if not arguments["debug"]:
         shutil.rmtree(Path("DAJIN_Results", ".tempdir", arguments["name"]))
 
@@ -137,14 +143,16 @@ def execute_batch_mode(arguments: dict[str]):
         # Set logging to export log to stderr and file
         config.reset_logging()
         path_logfile = config.get_logfile()
-        config.set_logging(path_logfile)
+        logger = config.set_logging(path_logfile)
+        logger.info(f"\N{runner} Start running DAJIN2 version {config.DAJIN_VERSION}")
+        logger.info(f"\N{Personal Computer} {' '.join(sys.argv)}")
 
-        # Start DAJIN2
+        # Run DAJIN2
         run_DAJIN2(groups, cache_urls_genome, is_control=True, num_workers=arguments["threads"])
         run_DAJIN2(groups, cache_urls_genome, is_control=False, num_workers=arguments["threads"])
 
         # Finish call
-        generate_report(name)
+        generate_report(name, logger)
         shutil.move(path_logfile, Path("DAJIN_Results", name))
         if not arguments["debug"]:
             shutil.rmtree(Path("DAJIN_Results", ".tempdir", name))
diff --git a/src/DAJIN2/utils/config.py b/src/DAJIN2/utils/config.py
index 1a9cd9ad..8553b67f 100644
--- a/src/DAJIN2/utils/config.py
+++ b/src/DAJIN2/utils/config.py
@@ -9,7 +9,7 @@
 
 from sklearn.exceptions import ConvergenceWarning
 
-DAJIN_VERSION = "0.4.6"
+DAJIN_VERSION = "0.5.0"
 DAJIN_RESULTS_DIR = Path("DAJIN_Results")
 TEMP_ROOT_DIR = Path(DAJIN_RESULTS_DIR, ".tempdir")
 
@@ -31,7 +31,7 @@ def get_logfile() -> Path:
     return Path(f"DAJIN2_log_{current_time}.txt")
 
 
-def set_logging(path_logfile: Path) -> None:
+def set_logging(path_logfile: Path) -> logging.Logger:
     log_format = "%(asctime)s, %(levelname)s, %(message)s"
     date_format = "%Y-%m-%d %H:%M:%S"
 
@@ -52,8 +52,10 @@ def handle_uncaught_exception(exc_type, exc_value, exc_traceback):
 
     sys.excepthook = handle_uncaught_exception
 
+    return logging.getLogger(__name__)
 
-def reset_logging():
+
+def reset_logging() -> None:
     """
     Reset the logging system to its default state.
     """
diff --git a/src/DAJIN2/utils/cssplits_handler.py b/src/DAJIN2/utils/cssplits_handler.py
index c66795de..2ebd27f2 100644
--- a/src/DAJIN2/utils/cssplits_handler.py
+++ b/src/DAJIN2/utils/cssplits_handler.py
@@ -1,6 +1,8 @@
 from __future__ import annotations
 import re
 import cstag
+import numpy as np
+import ruptures as rpt
 
 
 def find_n_boundaries(cssplits: list[str]) -> tuple[int, int]:
@@ -28,14 +30,14 @@ def find_n_boundaries(cssplits: list[str]) -> tuple[int, int]:
 ###########################################################
 
 
-def _reverse_cssplits(cssplits: list) -> list:
+def _reverse_cssplits(cssplits: list[str]) -> list[str]:
     for i, cs in enumerate(cssplits):
         if cs.startswith("+"):
             cssplits[i] = "+" + "|".join(cs.split("|")[::-1])
     return cssplits[::-1]
 
 
-def _realign_insertion(cssplits: list) -> list:
+def _realign_insertion(cssplits: list[str]) -> list[str]:
     for i, cs in enumerate(cssplits):
         if not cs.startswith("+"):
             continue
@@ -49,7 +51,7 @@ def _realign_insertion(cssplits: list) -> list:
     return cssplits
 
 
-def _complement_cssplit(cssplits: list) -> list:
+def _complement_cssplit(cssplits: list[str]) -> list[str]:
     comp = {"A": "T", "C": "G", "G": "C", "T": "A", "N": "N", "a": "t", "c": "g", "g": "c", "t": "a", "n": "n"}
     for i, cs in enumerate(cssplits):
         op = cs[0]
@@ -76,76 +78,79 @@ def revcomp_cssplits(cssplits: list[str]) -> list[str]:
 ###########################################################
 
 
-def add_match_operator_to_n(cssplits: list[str]) -> list[str]:
-    """Add "=" (match operator) to the sequences that start with "N"."""
-    return ["=" + seq if seq.startswith("N") or seq.startswith("n") else seq for seq in cssplits]
-
+def _add_match_operator_to_n(cssplits: list[str]) -> list[str]:
+    """Add "=" (match operator) to the sequences with "N"."""
+    cssplits_add_match_op = []
+    for cs in cssplits:
+        if cs.startswith("N") or cs.startswith("n"):
+            cssplits_add_match_op.append("=" + cs)
+        elif cs.startswith("+") and (cs[-1] == "N" or cs[-1] == "n"):
+            cs_ins = cs.split("|")
+            cs_last = "=" + cs_ins[-1]
+            cs = "|".join(cs_ins[:-1]) + "|" + cs_last
+            cssplits_add_match_op.append(cs)
+        else:
+            cssplits_add_match_op.append(cs)
+    return cssplits_add_match_op
 
-def format_insertion(cs: str) -> str:
-    """Reformat insertion sequence by consolidating variants."""
-    # Remove "+" and split by "|"
-    variants = cs.replace("+", "").split("|")
-    # Extract the last character from all but the last variant
-    consolidated_variants = "".join([variant[-1] for variant in variants[:-1]])
-    # Combine consolidated variants with the last variant
-    formatted_insertion = "+" + consolidated_variants + variants[-1]
 
-    return formatted_insertion
+def _split_cssplits_by_delimiter(cssplits: list[str]) -> list[str]:
+    cssplits_break = []
+    for cs in cssplits:
+        if cs.startswith("+"):
+            cssplits_break.extend(cs.split("|"))
+        else:
+            cssplits_break.append(cs)
+    return cssplits_break
 
 
-def concatenate_cssplits(cssplits: list[str]) -> str:
-    """Concatenate list of sequences based on certain variants."""
+def _combine_cssplits_by_prefix(cssplits: list[str]) -> list[str]:
     if not cssplits:
-        return ""
-    concatenated = []
-    prev = cssplits[0]
-    if prev[0] == "+":
-        concatenated.append(format_insertion(prev))
-    else:
-        concatenated.append(prev)
-    for prev, current in zip(cssplits, cssplits[1:]):
-        if prev[0] == current[0] and current[0] in {"=", "-"}:
-            concatenated.append(current[1:])
-        elif current[0] == "+":
-            concatenated.append(format_insertion(current))
+        return []
+
+    combined_cssplits = []
+    prev_prefix: str = cssplits[0][0]
+    current_cs: list[str] = [cssplits[0][1:]]
+
+    for item in cssplits[1:]:
+        current_prefix, cs = item[0], item[1:]
+        if prev_prefix == current_prefix:
+            if current_prefix == "*":
+                current_cs.append(item)
+            else:
+                current_cs.append(cs)
         else:
-            concatenated.append(current)
+            combined_cssplits.append(prev_prefix + "".join(current_cs))
+            prev_prefix = current_prefix
+            current_cs = [cs]
 
-    return "".join(concatenated)
+    combined_cssplits.append(prev_prefix + "".join(current_cs))
+    return combined_cssplits
 
 
-def standardize_case(sequence: str) -> str:
-    """Standardize the case of characters based on preceding variants."""
+def _standardize_case(cssplits: list[str]) -> list[str]:
+    """Standardize the case of characters based on mutation types."""
     transformations = {"*": str.lower, "-": str.lower, "+": str.lower, "~": str.lower, "=": str.upper}
 
     transformed = []
-    i = 0
-    while i < len(sequence):
-        char = sequence[i]
-        trans_func = transformations.get(char)
-        if trans_func:
-            start = i + 1
-            while start < len(sequence) and sequence[start].isalpha():
-                start += 1
-            transformed.append(char + trans_func(sequence[i + 1 : start]))
-            i = start - 1
-        else:
-            transformed.append(char)
-        i += 1
+    for cs in cssplits:
+        prefix = cs[0]
+        trans_func = transformations.get(prefix)
+        transformed.append(prefix + trans_func(cs[1:]))
 
-    return "".join(transformed)
+    return transformed
 
 
 def convert_cssplits_to_cstag(cssplits: list[str]) -> str:
-    cssplits = add_match_operator_to_n(cssplits)
-    cssplits_concatenated = concatenate_cssplits(cssplits)
-    return standardize_case(cssplits_concatenated)
-
-
-# convert position weight matrix (cons_pergentage) to sequence
+    cssplits_matched_op = _add_match_operator_to_n(cssplits)
+    cssplits_splitted = _split_cssplits_by_delimiter(cssplits_matched_op)
+    cssplits_combined = _combine_cssplits_by_prefix(cssplits_splitted)
+    return "".join(_standardize_case(cssplits_combined))
 
 
 def call_sequence(cons_percentage: list[dict[str, float]], sep: str = "") -> str:
+    """convert position weight matrix (cons_pergentage) to sequence"""
+
     consensus_sequence = []
     for cons_per in cons_percentage:
         cssplits = max(cons_per, key=cons_per.get)
@@ -156,63 +161,103 @@ def call_sequence(cons_percentage: list[dict[str, float]], sep: str = "") -> str
 
 
 ###########################################################
-# detect_insertion_within_deletion
+# reallocate_insertion_within_deletion
 ###########################################################
 
 
-def is_start_of_consecutive_deletions(cssplits: list[str], i: int, del_range: int = 3) -> bool:
-    """Determine if the current index is the start of a consective deletion."""
-    if i + del_range > len(cssplits):
-        return False
-
-    for j in range(del_range):
-        if not cssplits[i + j].startswith("-"):
-            return False
-
-    return True
-
-
-def get_range_of_inclusive_deletions(cssplits: list[str], del_range: int = 3, distance: int = 10) -> set[int]:
-    inclusive_deletion_range: set[int] = set()
-    n: int = len(cssplits)
-    i: int = 0
-    while i < n:
-        # Check for start of deletion
-        if cssplits[i].startswith("-"):
-            start: int = i
-            while i < n and cssplits[i].startswith("-"):
-                i += 1
-            # We have found the end of a deletion range
-            end: int = i
-            # Check if deletion length is enough to be considered
-
-            if end - start >= del_range:
-                # Now check for matching sequence within 'distance'
-                match_count: int = 0
-                while i < n and match_count <= distance:
-                    if cssplits[i].startswith("="):
-                        match_count += 1
-                    elif cssplits[i].startswith("-"):
-                        # If another deletion is found, reset match count and update range
-                        next_start: int = i
-                        while i < n and cssplits[i].startswith("-"):
-                            i += 1
-                            next_end: int = i
-                        if next_end - next_start >= del_range:
-                            inclusive_deletion_range.update(range(start, next_end))
-
-                            start = next_start
-                            end = next_end
-                            match_count = 0
-                            break
-                    i += 1
-
-        i += 1
-
-    return inclusive_deletion_range
-
-
-def adjust_cs_insertion(cs: str) -> str:
+def _extract_candidate_index_of_large_deletions(
+    cssplits: list[str], bin_size: int = 500, percentage: int = 50
+) -> list[dict[str, int]]:
+    deletion_counts = []
+    for start in range(0, len(cssplits), bin_size // 10):
+        end = start + bin_size
+        counts_deletion = "".join(cssplits[start:end]).count("-")
+        deletion_counts.append((start, counts_deletion))
+
+    range_of_large_deletions = []
+    on_the_peak = False
+    start = -1
+    end = -1
+    for i, count in deletion_counts:
+        del_percentage = count / bin_size * 100
+        if del_percentage > percentage and on_the_peak is False:
+            start = i
+            on_the_peak = True
+        if del_percentage <= percentage and on_the_peak is True:
+            end = i
+            on_the_peak = False
+        if start != -1 and end != -1:
+            range_of_large_deletions.append({"start": start, "end": end})
+            start = -1
+            end = -1
+
+    return range_of_large_deletions
+
+
+def _extract_break_points_of_large_deletions(
+    cssplits: list[str], range_of_large_deletions: list[dict[str, int]], bin_size: int = 500
+):
+    cssplits_match_bool = np.array([1 if cs.startswith("=") else 0 for cs in cssplits], dtype=bool)
+
+    break_points = []
+    for range_of_deletion in range_of_large_deletions:
+        break_point = {}
+        for key, index in range_of_deletion.items():
+            window_start = index
+            window_end = index + bin_size
+
+            # Find breakpoints(bps). First, predict approximate breakpoints using the signal of the window.
+            signal = cssplits_match_bool[window_start:window_end]
+            bps = rpt.Window(model="ar").fit(signal).predict(n_bkps=10)[:-1]
+
+            # Aling breakpoints to the original index
+            if bps:
+                approximate_bp: int = bps[0] if key == "start" else bps[-1]
+            else:
+                approximate_bp: int = 0 if key == "start" else bin_size
+
+            # Find the exact breakpoint because the approximate breakpoint may not be accurate
+            signal_subset = signal[max(0, approximate_bp - 50) : min(approximate_bp + 50, len(signal))]
+            bp = -1
+            for i, sig in enumerate(signal_subset):
+                if key == "start" and sig == False:
+                    bp = i + max(0, approximate_bp - 50) + window_start
+                    break
+                if key == "end" and sig == True:
+                    bp = i + max(0, approximate_bp - 50) + window_start - 1
+                    break
+            if bp == -1:
+                bp = approximate_bp
+
+            # Store the breakpoint
+            break_point[key] = bp
+
+        if break_point["start"] == break_point["end"]:
+            continue
+
+        break_points.append(break_point)
+
+    return break_points
+
+
+def _convert_break_points_to_index(break_points: list[dict[str, int]]) -> set[int]:
+    index_of_large_deletions = set()
+    for break_point in break_points:
+        start = break_point["start"]
+        end = break_point["end"]
+        index_of_large_deletions |= set(range(start, end + 1))
+
+    return index_of_large_deletions
+
+
+def _get_index_of_large_deletions(cssplits: list[str], bin_size: int = 500, percentage: int = 50) -> set[int]:
+    range_of_large_deletions = _extract_candidate_index_of_large_deletions(cssplits, bin_size, percentage)
+    break_points = _extract_break_points_of_large_deletions(cssplits, range_of_large_deletions, bin_size)
+
+    return _convert_break_points_to_index(break_points)
+
+
+def _adjust_cs_insertion(cs: str) -> str:
     if cs.startswith("-"):
         return None
     elif cs.startswith("+"):
@@ -224,18 +269,19 @@ def adjust_cs_insertion(cs: str) -> str:
     return cs_insertion
 
 
-def reallocate_insertion_within_deletion(cssplit: str, del_range: int = 3, distance: int = 10) -> str:
+def reallocate_insertion_within_deletion(cssplits: list[str], bin_size: int = 500, percentage: int = 50) -> list[str]:
     """
     Since the mapping in minimap2 is local alignment, insertion bases within large deletions may be partially mapped to the reference genome and not detected as insertion bases. Therefore, update cssplits to detect insertions within large deletions as insertions.
     """
-    cssplits = cssplit.split(",")
     cssplits_updated = cssplits.copy()
 
-    range_of_inclusive_deletion: set[int] = get_range_of_inclusive_deletions(cssplits, del_range, distance)
+    index_of_large_deletions: set[int] = _get_index_of_large_deletions(
+        cssplits, bin_size=bin_size, percentage=percentage
+    )
 
     insertion_within_deletion = []
     for i, cs in enumerate(cssplits):
-        if i in range_of_inclusive_deletion:
+        if i in index_of_large_deletions:
             if cs.startswith("-"):
                 continue
             if cs.startswith("N") or cs.startswith("n"):
@@ -248,11 +294,11 @@ def reallocate_insertion_within_deletion(cssplit: str, del_range: int = 3, dista
                 cs_new = f"-{cs[-1]}"
             cssplits_updated[i] = cs_new
 
-            insertion_within_deletion.append(adjust_cs_insertion(cs))
+            insertion_within_deletion.append(_adjust_cs_insertion(cs))
         else:
             if not insertion_within_deletion:
                 continue
             cssplits_updated[i] = "".join(insertion_within_deletion) + cs
             insertion_within_deletion = []
 
-    return ",".join(cssplits_updated)
+    return cssplits_updated
diff --git a/src/DAJIN2/utils/fastx_handler.py b/src/DAJIN2/utils/fastx_handler.py
index 9a819865..58455ac2 100644
--- a/src/DAJIN2/utils/fastx_handler.py
+++ b/src/DAJIN2/utils/fastx_handler.py
@@ -4,23 +4,10 @@
 import gzip
 from pathlib import Path
 
-import mappy
-
-
-#################################################
-# Helper function
-#################################################
-
-
-def sanitize_filename(path_file: Path | str) -> str:
-    """
-    Sanitize the path_file by replacing invalid characters on Windows OS with '-'
-    """
-    path_file = str(path_file).lstrip()
-    if not path_file:
-        raise ValueError("Provided FASTA/FASTQ is empty or consists only of whitespace")
-    return re.sub(r'[\\/:?.,\'"<>| ]', "-", path_file)
+import random
 
+import mappy
+from DAJIN2.utils.io import sanitize_name
 
 #################################################
 # Extract filename
@@ -30,7 +17,7 @@ def sanitize_filename(path_file: Path | str) -> str:
 def extract_filename(path_fasta: Path | str) -> str:
     filename = Path(path_fasta).name
     filename = re.sub(r"\..*$", "", filename)  # Remove file extension
-    return sanitize_filename(filename)
+    return sanitize_name(filename)
 
 
 #################################################
@@ -39,7 +26,7 @@ def extract_filename(path_fasta: Path | str) -> str:
 
 
 def dictionize_allele(path_fasta: str | Path) -> dict[str, str]:
-    return {sanitize_filename(name): seq.upper() for name, seq, _ in mappy.fastx_read(str(path_fasta))}
+    return {sanitize_name(name): seq.upper() for name, seq, _ in mappy.fastx_read(str(path_fasta))}
 
 
 #################################################
@@ -48,11 +35,12 @@ def dictionize_allele(path_fasta: str | Path) -> dict[str, str]:
 
 
 def export_fasta_files(TEMPDIR: Path, FASTA_ALLELES: dict, NAME: str) -> None:
-    """+ Save multiple FASTAs in separate single-FASTA format files."""
+    """Save multiple FASTAs in separate single-FASTA format files."""
     for identifier, sequence in FASTA_ALLELES.items():
+        identifier = sanitize_name(identifier)
         contents = "\n".join([">" + identifier, sequence]) + "\n"
-        output_fasta = Path(TEMPDIR, NAME, "fasta", f"{identifier}.fasta")
-        output_fasta.write_text(contents)
+        path_output_fasta = Path(TEMPDIR, NAME, "fasta", f"{identifier}.fasta")
+        path_output_fasta.write_text(contents)
 
 
 #################################################
@@ -60,15 +48,15 @@ def export_fasta_files(TEMPDIR: Path, FASTA_ALLELES: dict, NAME: str) -> None:
 #################################################
 
 
-def extract_extention(path_file: Path) -> str:
-    suffixes = path_file.suffixes
+def extract_extention(path_file: str | Path) -> str:
+    suffixes = Path(path_file).suffixes
     return "".join(suffixes)
 
 
-def is_gzip_file(path_file: Path) -> bool:
+def is_gzip_file(path_file: str | Path) -> bool:
     """Check if a file is a GZip compressed file."""
     try:
-        with path_file.open("rb") as f:
+        with Path(path_file).open("rb") as f:
             return f.read(2) == b"\x1f\x8b"
     except IOError:
         return False
@@ -86,9 +74,52 @@ def save_fastq_as_gzip(TEMPDIR: Path, path_fastx: list[Path], barcode: str) -> N
                     merged_file.write(f.read().encode())
 
 
-def save_concatenated_fastx(TEMPDIR: Path, directory: str) -> None:
+def save_concatenated_fastx(TEMPDIR: Path, path_directory: Path, name: str) -> None:
     fastx_suffix = {".fa", ".fq", ".fasta", ".fastq", ".fa.gz", ".fq.gz", ".fasta.gz", ".fastq.gz"}
-    path_directory = Path(directory)
-    barcode = path_directory.stem
     path_fastx = [path for path in path_directory.iterdir() if extract_extention(path) in fastx_suffix]
-    save_fastq_as_gzip(TEMPDIR, path_fastx, barcode)
+    save_fastq_as_gzip(TEMPDIR, path_fastx, name)
+
+
+#################################################
+# save_subsetted_fastx
+#################################################
+
+
+def read_lines(path_file) -> list[str]:
+    if is_gzip_file(path_file):
+        with gzip.open(path_file, "rt") as f:
+            return [line.strip() for line in f]
+    else:
+        with open(path_file, "r") as f:
+            return [line.strip() for line in f]
+
+
+def parse_fastq(fastq) -> list[dict]:
+    fastq_parsed = []
+    iterator = iter(fastq)
+    try:
+        while True:
+            header = next(iterator).strip()
+            sequence = next(iterator).strip()
+            annotate = next(iterator).strip()
+            quality = next(iterator).strip()
+            fastq_parsed.append({"header": header, "sequence": sequence, "annotate": annotate, "quality": quality})
+    except StopIteration:
+        pass
+    return fastq_parsed
+
+
+def save_subsetted_fastx(path_fastq: str | Path, num_reads: int = 10_000) -> None:
+    """If the number of control reads is too high, it unnecessarily slows down the computation speed. Therefore, we perform random sampling to reduce the number of reads to below 10,000."""
+
+    fastq: list[str] = read_lines(path_fastq)
+    if len(fastq) // 4 <= 10_000:
+        return None
+
+    fastq: list[dict] = parse_fastq(fastq)
+
+    random.seed(1)
+    fastq_subset = random.sample(fastq, num_reads)
+    with gzip.open(path_fastq, "wt") as f:
+        for read in fastq_subset:
+            f.write(f"{read['header']}\n{read['sequence']}\n{read['annotate']}\n{read['quality']}\n")
diff --git a/src/DAJIN2/utils/io.py b/src/DAJIN2/utils/io.py
index 66399bd4..cd52b613 100644
--- a/src/DAJIN2/utils/io.py
+++ b/src/DAJIN2/utils/io.py
@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import re
 import csv
 import json
 import pickle
@@ -166,6 +167,18 @@ def convert_to_posix(path: str) -> str:
     return path
 
 
+def sanitize_name(file_name: Path | str) -> str:
+    """
+    Sanitize the file name by replacing invalid characters on Windows OS with '-'
+    """
+    file_name = str(file_name).strip()
+    if not file_name:
+        raise ValueError("Provided name is empty or consists only of whitespace")
+    forbidden_chars = r'[<>:"/\\|?*\x00-\x1F .]'
+
+    return re.sub(forbidden_chars, "-", file_name)
+
+
 ###########################################################
 # Cache control hash
 ###########################################################


Test patch:
diff --git a/.github/workflows/cli-tests.yml b/.github/workflows/cli-tests.yml
new file mode 100644
index 00000000..19ac50b1
--- /dev/null
+++ b/.github/workflows/cli-tests.yml
@@ -0,0 +1,53 @@
+name: cli-tests
+
+on:
+  push:
+  schedule:
+    - cron: "0 2 * * 1" # The task runs at 2 a.m. every Monday.
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      max-parallel: 1
+      matrix:
+        os: [ubuntu-latest]
+        python-version: ['3.10']
+    name: Python ${{ matrix.python-version }} on ${{ matrix.os }}
+
+    defaults:
+      run:
+        shell: bash -l {0}
+
+    steps:
+      - uses: actions/checkout@v4
+      - uses: conda-incubator/setup-miniconda@v3
+        with:
+          python-version: ${{ matrix.python-version }}
+          miniconda-version: "latest"
+          activate-environment: test-env
+          channels: bioconda, conda-forge, defaults, anaconda
+
+      - name: Setup conda environment on ubuntu-latest🐧
+        if: runner.os == 'Linux'
+        run: |
+          conda create -y -n env-dajin2 python=${{ matrix.python-version }}
+
+      - name: Install dependencies
+        run: |
+          conda activate env-dajin2
+
+          # ==========================================
+          # Use pip because conda does not support pyproject.toml.
+          # ==========================================
+          pip install -U pip pytest
+          pip install -e .
+
+      ##########################################
+      # Execute tests
+      ##########################################
+
+      - name: Run tests version
+        run: |
+          conda activate env-dajin2
+          sh tests/test_version.sh
diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index ee5ac730..5c8f9852 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -43,16 +43,20 @@ jobs:
           conda activate env-dajin2
           conda config --env --set subdir osx-64
 
+      - name: Install weasyprint on macos-latest🍎
+        if: runner.os == 'macOS'
+        run: |
+          conda install -n env-dajin2 weasyprint
+
       - name: Install dependencies
         run: |
           conda activate env-dajin2
 
           # ==========================================
-          # Remove kaleido from requirements.txt because it has a different name in conda (python-kaleido).
+          # Use pip because conda does not support pyproject.toml.
           # ==========================================
-          grep -v kaleido requirements.txt > requirements_.txt 
-          conda install --yes --file requirements_.txt
-          pip install -U pip pytest kaleido
+          pip install -U pip pytest
+          pip install -e .
 
       - name: Execute tests
         env:
diff --git a/tests/src/consensus/test_clust_subsetter.py b/tests/src/consensus/test_clust_formatter.py
similarity index 59%
rename from tests/src/consensus/test_clust_subsetter.py
rename to tests/src/consensus/test_clust_formatter.py
index fb4d4f4d..eda65d84 100644
--- a/tests/src/consensus/test_clust_subsetter.py
+++ b/tests/src/consensus/test_clust_formatter.py
@@ -1,42 +1,42 @@
-from DAJIN2.core.consensus.clust_formatter import subset_clust
+from DAJIN2.core.consensus.clust_formatter import downsample_by_label
 
 
-def test_basic_subset_clust():
+def test_basic_downsample_by_label():
     sample = [{"LABEL": 1, "DATA": i} for i in range(1500)]
-    subset = subset_clust(sample)
+    subset = downsample_by_label(sample)
     assert len(subset) == 1000  # Default num value is 1000
     assert all(item["LABEL"] == 1 for item in subset)
 
 
-def test_subset_clust_with_custom_num():
+def test_downsample_by_label_with_custom_num():
     sample = [{"LABEL": 1, "DATA": i} for i in range(1500)]
-    subset = subset_clust(sample, num=500)
+    subset = downsample_by_label(sample, num=500)
     assert len(subset) == 500
 
 
-def test_subset_clust_with_multiple_labels():
+def test_downsample_by_label_with_multiple_labels():
     sample = [{"LABEL": 1, "DATA": i} for i in range(500)] + [{"LABEL": 2, "DATA": i} for i in range(500)]
-    subset = subset_clust(sample)
+    subset = downsample_by_label(sample)
     assert len(subset) == 1000
     assert len([item for item in subset if item["LABEL"] == 1]) == 500
     assert len([item for item in subset if item["LABEL"] == 2]) == 500
 
 
-def test_subset_clust_with_less_than_num():
+def test_downsample_by_label_with_less_than_num():
     sample = [{"LABEL": 1, "DATA": i} for i in range(400)]
-    subset = subset_clust(sample, num=500)
+    subset = downsample_by_label(sample, num=500)
     assert len(subset) == 400
 
 
-def test_subset_clust_with_multiple_labels_with_less_than_nu():
+def test_downsample_by_label_with_multiple_labels_with_less_than_nu():
     sample = [{"LABEL": 1, "DATA": i} for i in range(500)] + [{"LABEL": 2, "DATA": i} for i in range(500)]
-    subset = subset_clust(sample, num=200)
+    subset = downsample_by_label(sample, num=200)
     assert len(subset) == 400
     assert len([item for item in subset if item["LABEL"] == 1]) == 200
     assert len([item for item in subset if item["LABEL"] == 2]) == 200
 
 
-def test_subset_clust_with_empty_sample():
+def test_downsample_by_label_with_empty_sample():
     sample = []
-    subset = subset_clust(sample)
+    subset = downsample_by_label(sample)
     assert len(subset) == 0
diff --git a/tests/src/preprocess/test_midsv_caller.py b/tests/src/preprocess/test_midsv_caller.py
index a8dd507e..0b945e59 100644
--- a/tests/src/preprocess/test_midsv_caller.py
+++ b/tests/src/preprocess/test_midsv_caller.py
@@ -1,10 +1,8 @@
 from __future__ import annotations
 
 import pytest
-import midsv
 
 from DAJIN2.core.preprocess.midsv_caller import has_inversion_in_splice
-from DAJIN2.core.preprocess.midsv_caller import extract_qname_of_map_ont
 from DAJIN2.core.preprocess.midsv_caller import replace_internal_n_to_d
 from DAJIN2.core.preprocess.midsv_caller import convert_flag_to_strand
 from DAJIN2.core.preprocess.midsv_caller import convert_consecutive_indels_to_match
@@ -61,47 +59,6 @@ def test_has_inversion_in_splice_random_deletion():
 #     test_cigar = [s[5] for s in test if not s[0].startswith("@")]
 #     assert not has_inversion_in_splice(test_cigar[0])
 
-###########################################################
-# extract_qname_of_map_ont
-###########################################################
-
-
-def test_extract_qname_of_map_ont_simulation():
-    sam_ont = [["@header"], ["read1", "", "", "", "", "5M"]]
-    sam_splice = [["@header"], ["read1", "", "", "", "", "5M"]]
-    qname_of_map_ont = extract_qname_of_map_ont(iter(sam_ont), iter(sam_splice))
-    assert qname_of_map_ont == {"read1"}
-
-    # Large deletion
-    sam_ont = [["@header"], ["read1", "", "", "", "", "5M"], ["read1", "", "", "", "", "5M"]]
-    sam_splice = [["@header"], ["read1", "", "", "", "", "5M100D5M"]]
-    qname_of_map_ont = extract_qname_of_map_ont(iter(sam_ont), iter(sam_splice))
-    assert qname_of_map_ont == set()
-
-    # Inversion
-    sam_ont = [
-        ["@header"],
-        ["read1", "", "", "", "", "5M"],
-        ["read1", "", "", "", "", "100M"],
-        ["read1", "", "", "", "", "5M"],
-    ]
-    sam_splice = [["@header"], ["read1", "", "", "", "", "5M100I100N5M"]]
-    qname_of_map_ont = extract_qname_of_map_ont(iter(sam_ont), iter(sam_splice))
-    assert qname_of_map_ont == {"read1"}
-
-    # Inversion with single read in map-ont
-    sam_ont = [["@header"], ["read1", "", "", "", "", "5M10N5M"]]
-    sam_splice = [["@header"], ["read1", "", "", "", "", "5M100I100N5M"]]
-    qname_of_map_ont = extract_qname_of_map_ont(iter(sam_ont), iter(sam_splice))
-    assert qname_of_map_ont == {"read1"}
-
-
-def test_extract_qname_of_map_ont_real():
-    sam_ont = list(midsv.read_sam(Path("tests", "data", "preprocess", "midsv_caller", "stx2-ont_deletion.sam")))
-    sam_splice = list(midsv.read_sam(Path("tests", "data", "preprocess", "midsv_caller", "stx2-splice_deletion.sam")))
-    qname_of_map_ont = extract_qname_of_map_ont(iter(sam_ont), iter(sam_splice))
-    assert qname_of_map_ont == {"stx2-small-deletion"}
-
 
 ###########################################################
 # replace n to d
@@ -201,7 +158,7 @@ def test_convert_flag_to_strand(input_sample, expected_output):
         # no change
         ("=A,=C,=G,=T,=A", "=A,=C,=G,=T,=A"),
         # empty
-        ("", "")
+        ("", ""),
     ],
 )
 def test_convert_consecutive_indels_to_match(cons, expected):
diff --git a/tests/src/utils/test_cssplits_handler.py b/tests/src/utils/test_cssplits_handler.py
index c458264c..c8e9f3e5 100644
--- a/tests/src/utils/test_cssplits_handler.py
+++ b/tests/src/utils/test_cssplits_handler.py
@@ -1,18 +1,12 @@
 from __future__ import annotations
 
 import pytest
-from DAJIN2.utils.cssplits_handler import (
-    find_n_boundaries,
-    add_match_operator_to_n,
-    concatenate_cssplits,
-    standardize_case,
-    call_sequence,
-    is_start_of_consecutive_deletions,
-    # is_within_deletion,
-    get_range_of_inclusive_deletions,
-    adjust_cs_insertion,
-    reallocate_insertion_within_deletion,
-)
+from DAJIN2.utils import cssplits_handler
+
+
+###########################################################
+# find_n_boundaries
+###########################################################
 
 
 @pytest.mark.parametrize(
@@ -25,7 +19,7 @@
     ],
 )
 def test_find_n_boundaries(cssplits, expected):
-    assert find_n_boundaries(cssplits) == expected
+    assert cssplits_handler.find_n_boundaries(cssplits) == expected
 
 
 ###########################################################
@@ -33,76 +27,59 @@ def test_find_n_boundaries(cssplits, expected):
 ###########################################################
 
 
-def test_add_match_operator_to_n_empty_input():
-    assert add_match_operator_to_n([]) == []
-
-
-def test_add_match_operator_to_n_no_n_starting_strings():
-    cssplits = ["=A", "=C", "=G"]
-    assert add_match_operator_to_n(cssplits) == cssplits
-
-
-def test_add_match_operator_to_n_with_n_starting_strings():
-    cssplits = ["N", "=C", "N"]
-    expected_output = ["=N", "=C", "=N"]
-    assert add_match_operator_to_n(cssplits) == expected_output
-
-
-def test_concatenate_cssplits_empty_input():
-    assert concatenate_cssplits([]) == ""
-
-
-def test_concatenate_cssplits_single_element():
-    assert concatenate_cssplits(["=ACGT"]) == "=ACGT"
-
-
-def test_concatenate_cssplits_same_symbols():
-    cssplits = ["=AC", "=GT", "-CC", "-GG"]
-    assert concatenate_cssplits(cssplits) == "=ACGT-CCGG"
-
-
-def test_concatenate_cssplits_plus_symbol():
-    cssplits = ["+A|*GC|=T", "+A|+C|=G"]
-    assert concatenate_cssplits(cssplits) == "+AC=T+AC=G"
-
-
-def test_concatenate_cssplits_plus_symbol_torio():
-    cssplits = ["+A|*GC|=T", "+A|+C|=G", "+A|+C|=G"]
-    assert concatenate_cssplits(cssplits) == "+AC=T+AC=G+AC=G"
-
-
-def test_concatenate_cssplits_mixed():
-    cssplits = ["=AC", "=GT", "+A|+C|=G", "-CC", "-GG"]
-    assert concatenate_cssplits(cssplits) == "=ACGT+AC=G-CCGG"
-
-
-def test_standardize_case_empty_input():
-    assert standardize_case("") == ""
-
-
-def test_standardize_case_special_characters():
-    assert standardize_case("*AG-DEF+GHI~jkl=XYZ") == "*ag-def+ghi~jkl=XYZ"
+@pytest.mark.parametrize(
+    "cssplits, expected",
+    [
+        ([], []),
+        (["=A", "=C", "=G"], ["=A", "=C", "=G"]),
+        (["N", "=C", "N"], ["=N", "=C", "=N"]),
+        (["=A", "+G|+G|+G|N", "=T"], ["=A", "+G|+G|+G|=N", "=T"]),
+    ],
+)
+def test_add_match_operator_to_n(cssplits, expected):
+    assert cssplits_handler._add_match_operator_to_n(cssplits) == expected
 
 
-def test_standardize_case_no_special_characters():
-    assert standardize_case("=ABCDEF") == "=ABCDEF"
+# @pytest.mark.parametrize(
+#     "cssplits, expected",
+#     [
+#         ([], []),
+#         (["=A", "=C", "=G"], ["=A", "=C", "=G"]),
+#         (["+A|+A|*GC"], ["+A|+A|*GC"]),
+#         (["+A|*GC|=C"], ["+A|=C|=C"]),
+#     ],
+# )
+# def test_format_substitution_withtin_insertion(cssplits, expected):
+#     assert cssplits_handler._format_substitution_withtin_insertion(cssplits) == expected
 
 
-def test_standardize_case_mixed():
-    assert standardize_case("*abcDEF+ghi-JKL=XYZ123*AC") == "*abcdef+ghi-jkl=XYZ123*ac"
+@pytest.mark.parametrize(
+    "input_cssplits, expected",
+    [
+        ([], []),
+        (["=ACGT"], ["=ACGT"]),
+        (["=ACGT", "*GC", "=ACGT"], ["=ACGT", "*gc", "=ACGT"]),
+    ],
+)
+def test_standardize_case(input_cssplits, expected):
+    assert cssplits_handler._standardize_case(input_cssplits) == expected
 
 
 @pytest.mark.parametrize(
-    "input, expected",
+    "cssplits, expected",
     [
-        ("*abc", "*abc"),
-        ("-DEF", "-def"),
-        ("+GHI", "+ghi"),
-        ("=XYZ", "=XYZ"),
+        ([], ""),
+        (["=A", "=C", "=G", "=T"], "=ACGT"),
+        (["=A", "*GT", "=A"], "=A*gt=A"),
+        (["-C", "-G"], "-cg"),
+        (["*GC", "*TA"], "*gc*ta"),
+        (["+A|+A|=C", "=G"], "+aa=CG"),
+        (["+A|+C|=T", "+A|+C|=G"], "+ac=T+ac=G"),
+        (["=A", "=C", "+A|+C|=G", "-C", "-G", "*CG"], "=AC+ac=G-cg*cg"),
     ],
 )
-def test_standardize_case_parametrized(input, expected):
-    assert standardize_case(input) == expected
+def test_convert_cssplits_to_cstag(cssplits, expected):
+    assert cssplits_handler.convert_cssplits_to_cstag(cssplits) == expected
 
 
 ###########################################################
@@ -111,7 +88,7 @@ def test_standardize_case_parametrized(input, expected):
 
 
 @pytest.mark.parametrize(
-    "cons_percentage_by_key, expected_sequence",
+    "cons_percentage, expected_sequence",
     [
         ([{"=A": 1.0}, {"=T": 0.9, "-T": 0.1}], "AT"),  # match
         ([{"=A": 1.0}, {"-A": 0.9, "=A": 0.1}, {"=T": 1.0}], "AT"),  # deletion
@@ -124,75 +101,29 @@ def test_standardize_case_parametrized(input, expected):
         ([{"=A": 1.0}], "A"),
     ],
 )
-def test_call_sequence(cons_percentage_by_key, expected_sequence):
-    assert call_sequence(cons_percentage_by_key) == expected_sequence
+def test_call_sequence(cons_percentage, expected_sequence):
+    assert cssplits_handler.call_sequence(cons_percentage) == expected_sequence
 
 
 ###########################################################
-# get_range_of_inclusive_deletions
+# get_index_of_large_deletions
 ###########################################################
 
 
 @pytest.mark.parametrize(
     "cssplits, expected",
     [
-        (["-A", "-A", "-A", "=T", "-A", "-A", "-A"], {0, 1, 2, 3, 4, 5, 6}),
-        (
-            ["-A", "-A", "-A", "=T", "-A", "-A", "-A", "=G", "-A", "-A", "-A", "=T", "-A", "-A", "-A", "=G"],
-            {0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14},
-        ),
-        (
-            ["-A", "-A", "-A", "=T", "-A", "=T", "-A", "-A", "-A", "=G", "-G", "-G", "-G", "=G"],
-            {0, 1, 2, 3, 4, 5, 6, 7, 8},
-        ),
+        (["=T"] * 100 + ["-A"] * 300 + ["=T"] * 100, set(range(100, 400))),
     ],
 )
-def test_get_range_of_inclusive_deletions(cssplits, expected):
-    assert get_range_of_inclusive_deletions(cssplits) == expected
+def test_get_index_of_large_deletions(cssplits, expected):
+    assert cssplits_handler._get_index_of_large_deletions(cssplits) == expected
 
 
 ###########################################################
 # reallocate_insertion_within_deletion
 ###########################################################
 
-
-@pytest.mark.parametrize(
-    "cssplits, i, start_del_count, expected",
-    [
-        (0, ["-A", "-A", "-A"], 3, True),
-        (0, ["-A", "-A"], 3, False),
-        (1, ["-A", "-A", "-A"], 3, False),
-    ],
-)
-def test_is_start_of_consecutive_deletions(cssplits: list[str], i: int, start_del_count: int, expected: bool):
-    assert is_start_of_consecutive_deletions(i, cssplits, start_del_count) == expected
-
-
-# @pytest.mark.parametrize(
-#     "index, cssplits, del_range, distance, expected",
-#     [
-#         (0, ["-A", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "-T", "=G"], 1, 10, True),
-#         (11, ["-A", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "-T", "=G"], 1, 10, False),
-#         (0, ["-A", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "=C", "-T", "=G"], 1, 10, False),
-#         (0, ["-A", "=C", "=C", "=C", "=C", "=C", "*GC", "=C", "=C", "=C", "=C", "-T", "=G"], 1, 10, True),
-#         (0, ["-A", "=C", "=C", "=C", "=C", "=C", "-T", "=C", "=C", "=C", "=C", "-T", "=G"], 1, 10, True),
-#         (0, ["-A", "=C", "=C", "=C", "=C", "=C", "+A|=C", "=C", "=C", "=C", "=C", "-T", "=G"], 1, 10, True),
-#         (0, ["-A", "-A", "=C", "-T", "-A"], 2, 10, True),
-#     ],
-#     ids=[
-#         "-A has a 10-character match and is within the deletion cluster.",
-#         "-T is outside the deletion cluster.",
-#         "-A has a 11-character match and is outside the deletion cluster.",
-#         "-A is determined to be reset and within the deletion cluster due to a point mutation at the 6th character.",
-#         "-A is determined to be reset and within the deletion cluster due to a deletion at the 6th character.",
-#         "-A is determined to be reset and within the deletion cluster due to an insertion at the 6th character.",
-#         "-A is within the deletion cluster.",
-#     ],
-# )
-# def test_is_within_deletion(index: int, cssplits: list[str], del_range: int, distance: int, expected: bool):
-#     assert is_within_deletion(index, cssplits, del_range, distance) == expected
-
-
 @pytest.mark.parametrize(
     "cs, expected",
     [
@@ -204,27 +135,27 @@ def test_is_start_of_consecutive_deletions(cssplits: list[str], i: int, start_de
     ],
 )
 def test_adjust_cs_insertion(cs: str, expected: str):
-    assert adjust_cs_insertion(cs) == expected
+    assert cssplits_handler._adjust_cs_insertion(cs) == expected
 
 
-@pytest.mark.parametrize(
-    "input_str, expected_output",
-    [
-        ("-A,-A,-A,=C,=C,=C,-T,-T,-T,=G", "-A,-A,-A,-C,-C,-C,-T,-T,-T,+C|+C|+C|=G"),
-        ("-A,-A,-A,=C,=C,=C,=C,-T,-T,-T", "-A,-A,-A,=C,=C,=C,=C,-T,-T,-T"),
-        ("-A,-A,-A,N,=C,n,-T,-T,-T,=G", "-A,-A,-A,N,-C,n,-T,-T,-T,+N|+C|+n|=G"),
-        ("-A,-A,-A,=C,+T|+T|=C,=C,-T,-T,-T,=G", "-A,-A,-A,-C,-C,-C,-T,-T,-T,+C|+T|+T|+C|+C|=G"),
-        ("-A,-A,-A,=C,+T|+T|*CG,=C,-T,-T,-T,=G", "-A,-A,-A,-C,-C,-C,-T,-T,-T,+C|+T|+T|+G|+C|=G"),
-        ("-G,-G,-C,=A,=C,=C,*CA,=A,-T,-T,*AC", "-G,-G,-C,=A,=C,=C,*CA,=A,-T,-T,*AC"),
-    ],
-    ids=[
-        "insertion within deletion",
-        "4-character match",
-        "N and n",
-        "Insertion",
-        "Insertion followed by substitution",
-        "Should not be adjusted",
-    ],
-)
-def test_reallocate_insertion_within_deletion(input_str: str, expected_output: str):
-    assert reallocate_insertion_within_deletion(input_str, del_range=3, distance=3) == expected_output
+# @pytest.mark.parametrize(
+#     "input_str, expected_output",
+#     [
+#         ("-A,-A,-A,=C,=C,=C,-T,-T,-T,=G", "-A,-A,-A,-C,-C,-C,-T,-T,-T,+C|+C|+C|=G"),
+#         ("-A,-A,-A,=C,=C,=C,=C,-T,-T,-T", "-A,-A,-A,=C,=C,=C,=C,-T,-T,-T"),
+#         ("-A,-A,-A,N,=C,n,-T,-T,-T,=G", "-A,-A,-A,N,-C,n,-T,-T,-T,+N|+C|+n|=G"),
+#         ("-A,-A,-A,=C,+T|+T|=C,=C,-T,-T,-T,=G", "-A,-A,-A,-C,-C,-C,-T,-T,-T,+C|+T|+T|+C|+C|=G"),
+#         ("-A,-A,-A,=C,+T|+T|*CG,=C,-T,-T,-T,=G", "-A,-A,-A,-C,-C,-C,-T,-T,-T,+C|+T|+T|+G|+C|=G"),
+#         ("-G,-G,-C,=A,=C,=C,*CA,=A,-T,-T,*AC", "-G,-G,-C,=A,=C,=C,*CA,=A,-T,-T,*AC"),
+#     ],
+#     ids=[
+#         "insertion within deletion",
+#         "4-character match",
+#         "N and n",
+#         "Insertion",
+#         "Insertion followed by substitution",
+#         "Should not be adjusted",
+#     ],
+# )
+# def test_reallocate_insertion_within_deletion(input_str: str, expected_output: str):
+#     assert reallocate_insertion_within_deletion(input_str, del_range=3, distance=3) == expected_output
diff --git a/tests/src/utils/test_fastx_handler.py b/tests/src/utils/test_fastx_handler.py
index 8fe95f6b..7975d5c1 100644
--- a/tests/src/utils/test_fastx_handler.py
+++ b/tests/src/utils/test_fastx_handler.py
@@ -2,39 +2,8 @@
 
 import pytest
 
-from pathlib import Path
-
-
 from DAJIN2.utils import fastx_handler
 
-
-def test_sanitize_filename_with_valid_path():
-    assert fastx_handler.sanitize_filename("valid_name") == "valid_name"
-    assert fastx_handler.sanitize_filename(Path("valid/name")) == "valid-name"
-
-
-def test_sanitize_filename_with_invalid_characters():
-    assert fastx_handler.sanitize_filename("inva/lid:name?") == "inva-lid-name-"
-    assert fastx_handler.sanitize_filename(Path("/inva>lid|name.")) == "-inva-lid-name-"
-
-
-def test_sanitize_filename_with_whitespace():
-    assert fastx_handler.sanitize_filename("  leading_space") == "leading_space"
-    assert fastx_handler.sanitize_filename("trailing_space ") == "trailing_space-"
-
-
-def test_sanitize_filename_with_empty_string():
-    with pytest.raises(ValueError) as e:
-        fastx_handler.sanitize_filename(" ")
-    assert str(e.value) == "Provided FASTA/FASTQ is empty or consists only of whitespace"
-
-
-def test_sanitize_filename_with_empty_path():
-    with pytest.raises(ValueError) as e:
-        fastx_handler.sanitize_filename("")
-    assert str(e.value) == "Provided FASTA/FASTQ is empty or consists only of whitespace"
-
-
 ########################################################################
 # Extract basename
 ########################################################################
@@ -51,7 +20,7 @@ def test_extract_filename_fail():
     fastq_path = "tests/hoge/fuga/.fq.gz"
     with pytest.raises(ValueError) as e:
         fastx_handler.extract_filename(fastq_path)
-    assert str(e.value) == "Provided FASTA/FASTQ is empty or consists only of whitespace"
+    assert str(e.value) == "Provided name is empty or consists only of whitespace"
 
 
 def test_extract_filename_change_filename():
@@ -84,4 +53,4 @@ def test_dictionize_allele_empty():
     path_fasta = "tests/data/preprocess/format_input/fasta_empty.fa"
     with pytest.raises(ValueError) as e:
         fastx_handler.dictionize_allele(path_fasta)
-    assert str(e.value) == "Provided FASTA/FASTQ is empty or consists only of whitespace"
+    assert str(e.value) == "Provided name is empty or consists only of whitespace"
diff --git a/tests/src/utils/test_io.py b/tests/src/utils/test_io.py
index e4110ac2..425ea312 100644
--- a/tests/src/utils/test_io.py
+++ b/tests/src/utils/test_io.py
@@ -1,8 +1,39 @@
 from __future__ import annotations
 
-from src.DAJIN2.utils import io
 import os
 import json
+import pytest
+from pathlib import Path
+
+from src.DAJIN2.utils import io
+
+
+def test_sanitize_name_with_valid_path():
+    assert io.sanitize_name("valid_name") == "valid_name"
+    assert io.sanitize_name(Path("valid/name")) == "valid-name"
+
+
+def test_sanitize_name_with_invalid_characters():
+    assert io.sanitize_name("inva/lid:name?") == "inva-lid-name-"
+    assert io.sanitize_name(Path("/inva>lid|name.")) == "-inva-lid-name-"
+
+
+def test_sanitize_name_with_whitespace():
+    assert io.sanitize_name("  leading_space") == "leading_space"
+    assert io.sanitize_name("trailing space ") == "trailing-space"
+
+
+def test_sanitize_name_with_empty_string():
+    with pytest.raises(ValueError) as e:
+        io.sanitize_name(" ")
+    assert str(e.value) == "Provided name is empty or consists only of whitespace"
+
+
+def test_sanitize_name_with_empty_path():
+    with pytest.raises(ValueError) as e:
+        io.sanitize_name("")
+    assert str(e.value) == "Provided name is empty or consists only of whitespace"
+
 
 ########################################################################
 # Convert Path
diff --git a/tests/test_version.sh b/tests/test_version.sh
new file mode 100644
index 00000000..527b4afd
--- /dev/null
+++ b/tests/test_version.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+PIP_VERSION=$(echo $(DAJIN2 --version) | awk '{print $NF}')
+UTIL_VERSION=$(grep "DAJIN_VERSION" src/DAJIN2/utils/config.py | awk '{print $NF}' | tr -d '"')
+
+if [ -z "$PIP_VERSION" ]; then
+    echo "DAJIN2 is not installed"
+    exit 1
+fi
+
+if [ -z "$UTIL_VERSION" ]; then
+    echo "src/DAJIN2/utils/config.py is not found"
+    exit 1
+fi
+
+if [ "$PIP_VERSION" != "$UTIL_VERSION" ]; then
+    echo "Version mismatch: $PIP_VERSION != $UTIL_VERSION"
+    exit 1
+fi