Repository: bigbio/ibaqpy (language: python)

Issue:
update: todo



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index 19560d0..cf78424 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -40,4 +40,11 @@ jobs:
       run: |
         cd tests/
         python -m unittest
-
+    - name: Test commandline tool help
+      run: |
+        ibaqpy --help
+    - name: Test commandline normalization features
+      run: |
+        rm -rfv tests/PXD003947/PXD003947-peptides-norm.csv
+        rm -rfv tests/PXD003947/PXD003947-peptides-norm.parquet
+        ibaqpy features2peptides --parquet tests/PXD003947/PXD003947-feature.parquet --sdrf tests/PXD003947/PXD003947.sdrf.tsv --min_aa 7 --min_unique 2 --remove_ids data/contaminants_ids.tsv --remove_decoy_contaminants --remove_low_frequency_peptides --output tests/PXD003947/PXD003947-peptides-norm.csv --skip_normalization --nmethod median --pnmethod max_min --log2 --save_parquet
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 7415ce6..9f404e4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,3 +10,7 @@ venv
 /compute-all.sh
 /ibaqpy.egg-info/
 /ibaqpy_temp/
+/tests/PXD003947/IBAQ-QCprofile.pdf
+/tests/PXD003947/PXD003947-ibaq-norm.csv
+/tests/PXD003947/PXD003947-peptides-norm.csv
+/tests/PXD003947/PXD003947-peptides-norm.parquet
diff --git a/README.md b/README.md
index 35c0af7..7905cf4 100644
--- a/README.md
+++ b/README.md
@@ -6,20 +6,97 @@
 [![PyPI version](https://badge.fury.io/py/ibaqpy.svg)](https://badge.fury.io/py/ibaqpy)
 ![PyPI - Downloads](https://img.shields.io/pypi/dm/ibaqpy)
 
-iBAQ (intensity Based Absolute Quantification) determines the abundance of a protein by dividing the total precursor intensities by the number of theoretically observable peptides of the protein. The TPA (Total Protein Approach) value is determined by summing peptide intensities of each protein and then dividing by the molecular mass to determine the relative concentration of each protein. By using [ProteomicRuler](https://www.sciencedirect.com/science/article/pii/S1535947620337749), it is possible to calculate the protein copy number and absolute concentration. **ibaqpy** compute IBAQ values, TPA values, copy numbers and concentration for proteins starting from a feature parquet from [quantmsio](https://github.com/bigbio/quantms.io) and a [SDRF](https://github.com/bigbio/proteomics-sample-metadata) file. In addition, it supports the merging of iBAQ results from multiple datasets and the elimination of outliers and batch effects. This package provides multiple functions: 
+iBAQ (Intensity-Based Absolute Quantification) determines the abundance of a protein by dividing the total precursor intensities by the number of theoretically observable peptides of the protein [manuscript here](https://pubmed.ncbi.nlm.nih.gov/16219938/). ibaqpy is a Python package that computes iBAQ values starting from a feature parquet from [quantmsio](https://github.com/bigbio/quantms.io) and a [SDRF](https://github.com/bigbio/proteomics-sample-metadata) file. In addition, the package computes other ibaq values including rIBAQ, log2, and ppb.
 
-**Assemble features as peptides**: In features, the differences between different runs and different fractions in the sample were reduced, and then normalized after combining into peptides. See details in `commands/features2peptides`.
+ibaqpy also allows computing the TPA value (Total Protein Approach). TPA is determined by summing peptide intensities of each protein and then dividing by the molecular mass to determine the relative concentration of each protein. By using [ProteomicRuler](https://www.sciencedirect.com/science/article/pii/S1535947620337749), it is possible to calculate the protein copy number and absolute concentration. The OpenMS tool was used to calculate the theoretical molecular mass of each protein. Similar to the calculation of IBAQ, the TPA value of protein-group was the sum of its intensity divided by the sum of the theoretical molecular mass.
 
-**Assemble peptides as proteins**: In the sample, the protein was combined from its unique peptides, then the iBAQ value was calculated as the Intensity of the protein devided by the theoretical number of unique peptides cut by the enzyme. This command also normalized iBAQ as riBAQ. See details in `commands/peptides2proteins`.
+### Overview of ibaq-base values computation
 
-**Compute TPA**: Compute TPA values, protein copy numbers and concentration from the output file from script `commands/features2peptides`. See details in `commands/compute_tpa.py`.
+As mentioned before, ibaq values are calculated by dividing the total precursor intensities by the number of theoretically observable peptides of the protein. We use the following steps to calculate the iBAQ values:
 
-**Merge projects**: Merge ibaq results from multiple datasets. It consists of three steps: missing value imputation, outlier removal, and batch effect removal. See details in `commands/datasets_merger.py`.
+- _Observable peptides_, the protein sequence is digested in silico using a specific enzyme. The current version of this tool uses OpenMS method to load fasta file, and use [ProteaseDigestion](https://openms.de/current_doxygen/html/classOpenMS_1_1ProteaseDigestion.html) to enzyme digestion of protein sequences, and finally get the theoretical peptide number of each protein.
 
-**NOTE:** In all scripts and result files, *uniprot accession* is used as the protein identifier.
+- _Total precursor intensities_, the total intensity of a protein is calculated by summing the intensity of all peptides that belong to the protein. The intensity values are obtained from the feature parquet file in [quantms.io](https://github.com/bigbio/quantms.io).
+
+## TODO @PingZheng: Can you double-check that for protein group we multiply the intensity or divided? 
+
+> Note: If protein-group exists in the peptide intensity dataframe, the intensity of all proteins in the protein-group is summed based on the above steps, and then multiplied by the number of proteins in the protein-group.
+
+#### Other normalized IBAQ values  
+
+- `IbaqNorm` - normalize the ibaq values using the total ibaq of the sample `ibaq / sum(ibaq)`, the sum is applied for proteins in the same _sample + condition_.
+
+- `IbaqLog`  - The ibaq log is calculated as `10 + log10(IbaqNorm)`. This normalized ibaq value was developed [by ProteomicsDB Team](https://academic.oup.com/nar/article/46/D1/D1271/4584631).
+
+- `IbaqPpb` - The resulted IbaqNorm is multiplied by 100M `IbaqNorm * 100'000'000`. This method was developed originally [by PRIDE Team](https://www.nature.com/articles/s41597-021-00890-2).
+
+### From quantms to Ibaq values
 
 ![Ibaq](./data/ibaqpy.png "IBAQ")
 
+The output of quantms is converted into quantms.io feature file. A feature in quantms.io is the combination of the following columns: 
+
+- `ProteinName`: Protein name
+- `Peptidoform`: Peptide sequence including post-translation modifications `(e.g. .(Acetyl)ASPDWGYDDKN(Deamidated)GPEQWSK)`
+- `PrecursorCharge`: Precursor charge
+- `IsotopeLabelType`: Isotope label type
+- `Condition`: Condition label `(e.g. heart)`
+- `BioReplicate`: Biological replicate index `(e.g. 1)`
+- `Run`: Run index `(e.g. 1)`
+- `Fraction`: Fraction index `(e.g. 1)`
+- `Intensity`: Peptide intensity
+- `SampleID`: Sample ID `(e.g. PXD003947-Sample-3)`
+- `StudyID`: Study ID `(e.g. PXD003947)`. In most of the cases, the study ID is the same as the ProteomeXchange ID.
+
+In summary, each feature is the unique combination of a peptide sequence including modifications (peptidoform), precursor charge state, condition, biological replicate, run, fraction, isotopic label type, and a given intensity. In order to go from these features into protein ibaq values, the package does the following: 
+
+#### Data preprocessing
+
+## TODO @PingZheng: Can you confirm that the following steps are correct?
+
+In this section, ibaqpy will do: 
+- Remove lines where intensity or study condition is empty: This could happen in the following cases: 
+  - The DIA pipeline sometimes for some features releases intensities with value 0.
+  - The quantms.io do not contain feature information for some conditions. This extreme case could happen when not ID/Quant was found for a given condition during the analysis.
+- Parse the identifier of proteins and retain only unique peptides. 
+- Low-confidence proteins were removed according to the threshold of unique peptides: We use a thershold of 2 unique peptides to consider a protein as a high-confidence protein. This parameter is applied if not specified by the user, and the default value is 2. If users want to change this threshold, they can use the `--min_unique` parameter.
+- Filter decoy, contaminants, entrapment: Proteins with the following prefix are removed by default: `DECOY, CONTAMINANT, ENTRAPMENT` could be removed, by default, the filter is not applied. If users want to remove these proteins, they can use the `--remove_decoy_contaminants` parameter.
+- Filter user-specified proteins: The user can provide a list of protein identifiers to remove from the analysis using the `--remove_ids` parameter. The remove ids parameters will remove proteins from the analysis that could be potential to influence the intensity normalization. For example, ALBU_HUMAN could be over expressed in human tissues, and that is why we may want to remove it when analyzing tissue data.
+- Remove peptides with low frequency if `sample number > 1`: This parameter is applied always unless the user specifies the `--remove_low_frequency_peptides` parameter. The default value is 20% of the samples. If users want to change this threshold, they can use the `--remove_low_frequency_peptides` parameter.
+- Intensity transformation to log: The user can specify the `--log2` parameter to transform the peptide intensity values to log2 before normalization.
+
+> Note: At the moment, ibaqpy computes the ibaq values only based on unique peptides. Shared peptides are discarded. However, if a group of proteins share the same unique peptides (e.g., Pep1 -> Prot1;Prot2 and Pep2 -> Prot1;Prot2), the intensity of the proteins is summed and divided by the number of proteins in the group.
+
+#### Feature normalization
+
+## TODO @PingZheng: Can you confirm that the following steps are correct?
+
+In this section, ibaqpy corrects the intensity of each feature MS runs in the sample, eliminating the effect between runs (including technique repetitions). It will do:
+
+- When `MS runs > 1` in the sample, the `mean` of all average(`mean`, `median` or `iqr`) in each MS run is calculated(SampleMean)
+- The ratio between SampleMean and the average MS run is used as a reference to scale the original intensity
+
+#### Features to peptides
+
+A peptidoform is a combination of a `PeptideSequence(Modifications) + Charge + BioReplicate + Fraction` (among other features), and a peptide is a combination of a `PeptideSequence(Canonical) + BioReplicate`. ibaqpy will do:
+
+- Select peptidoforms with the highest intensity across different modifications, fractions, and technical replicates
+- Merge peptidoforms across different charges and combined into peptides. In order to merge peptidoforms, the package will applied the `sum` of the intensity values of the peptidoforms.
+
+#### Peptide Intensity Normalization
+
+Finally, the intensity of the peptide was normalized globally by `globalMedian`,`conditionMedian`. In the sample, the protein was combined with its unique peptides. 
+
+#### Computing Ibaq and TPA 
+
+The iBAQ value was calculated as the Intensity of the protein divided by the theoretical number of unique peptides cut by the enzyme. This command also normalized iBAQ as riBAQ. See details in `peptides2proteins`.
+
+Compute TPA values, protein copy numbers and concentration from the output file from script `commands/features2peptides`. 
+
+**Merge projects**: Merge ibaq results from multiple datasets. It consists of three steps: missing value imputation, outlier removal, and batch effect removal. 
+
+> Note: In all scripts and result files, *uniprot accession* is used as the protein identifier.
+
 ### How to install ibaqpy
 
 Ibaqpy is available in PyPI and can be installed using pip:
@@ -49,9 +126,9 @@ You can install the package from code:
 >$ python setup.py install
 ```
 
-### Collecting intensity files
+### Collecting intensity files from quantms.org 
 
-Absolute quantification files has been store in the following url:
+Absolute quantification files have been stored in the following url:
 
 ```
 http://ftp.pride.ebi.ac.uk/pub/databases/pride/resources/proteomes/absolute-expression/quantms-data/
@@ -61,10 +138,13 @@ Inside each project reanalysis folder, the folder proteomicslfq contains the mss
 
 E.g. http://ftp.pride.ebi.ac.uk/pub/databases/pride/resources/proteomes/absolute-expression/quantms-data/MSV000079033.1/MSV000079033.1-bd44c7e3-654c-444d-9e21-0f701d6dac94.feature.parquet
 
-### Assemble features as peptides
+### Major commands  
+
+
+#### Features to peptides
 
 ```asciidoc
-python commands/features2peptides.py -p tests/PXD003947/PXD003947-featrue.parquet -s tests/PXD003947/PXD003947.sdrf.tsv --remove_ids data/contaminants_ids.tsv --remove_decoy_contaminants --remove_low_frequency_peptides --output tests/PXD003947/PXD003947-peptides-norm.csv
+ibaqpy features2peptides -p tests/PXD003947/PXD003947-featrue.parquet -s tests/PXD003947/PXD003947.sdrf.tsv --remove_ids data/contaminants_ids.tsv --remove_decoy_contaminants --remove_low_frequency_peptides --output tests/PXD003947/PXD003947-peptides-norm.csv --log2 --save_parquet
 ```
 ```asciidoc
 Usage: features2peptides.py [OPTIONS]
@@ -87,74 +167,28 @@ Options:
                                   properties for normalization
   --skip_normalization            Skip normalization step
   --nmethod TEXT                  Normalization method used to normalize
-                                  feature intensities for tec
+                                  feature intensities for all samples
                                   (options: mean, median, iqr, none)
   --pnmethod TEXT                 Normalization method used to normalize
                                   peptides intensities for all samples
-                                  (options: globalMedian,conditionMedian,none)
+                                  (options: globalMedian,conditionMedian)
   --log2                          Transform to log2 the peptide intensity
                                   values before normalization
   --save_parquet                  Save normalized peptides to parquet
   --help                          Show this message and exit.
 ```
 
-Peptide normalization starts from the peptides dataframe extracted from feature parquet. It contains the following columns:
-
-- ProteinName: Protein name
-- PeptideSequence: Peptide sequence including post-translation modifications `(e.g. .(Acetyl)ASPDWGYDDKN(Deamidated)GPEQWSK)`
-- PrecursorCharge: Precursor charge
-- FragmentIon: Fragment ion
-- ProductCharge: Product charge
-- IsotopeLabelType: Isotope label type
-- Condition: Condition label `(e.g. heart)`
-- BioReplicate: Biological replicate index `(e.g. 1)`
-- Run: Run index `(e.g. 1)`
-- Fraction: Fraction index `(e.g. 1)`
-- Intensity: Peptide intensity
-- Reference: Name of the RAW file containing the peptide intensity `(e.g. Adult_Heart_Gel_Elite_54_f16)`
-- SampleID: Sample ID `(e.g. PXD003947-Sample-3)`
-- StudyID: Study ID `(e.g. PXD003947)`. In most of the cases the study ID is the same as the ProteomeXchange ID.
-
-#### 1. Data preprocessing
 
-In this section, ibaqpy will do: 
-- Remove lines where intensity or study condition is empty
-- Parse the identifier of proteins and retain only unique peptides
-- Low confidence proteins were removed according to the threshold of unique peptides (optional)
-- Filter decoy, contaminants, entrapments, and user-specified proteins (optional)
-- Remove peptides with low frequency if `sample number > 1` (optional)
-- Data logization (optional)
-
-#### 2. Feature normalization
-
-In this section, ibaqpy corrects the intensity of each MS runs in the sample, eliminating the effect between runs (including technique repetitions). It will do:
-
-- When `MS runs > 1` in the sample, the `mean` of all average(`mean`, `median` or `iqr`) in each MS run is calculated(SampleMean)
-- The ratio between SampleMean and the average MS run is used as reference to scale the original intensity
-
-#### 3. Assembly features to peptides
-
-A peptidoform is a combination of a `PeptideSequence(Modifications) + Charge + BioReplicate + Fraction`, and a peptide is a combination of a `PeptideSequence(Canonical) + BioReplicate`. ibaqpy will do:
-
-- Select peptidoforms with the highest intensity across different modifications, fractions, and technical repeats
-- Merge peptidoforms across different charges and combined into peptides
-
-#### 4. Peptide Intensity Normalization
-
-Finally, the intensity of the peptide was normalized globally by `globalMedian`,`conditionMedian`.
-
-
-### Compute IBAQ
-IBAQ is an absolute quantitative method based on strength that can be used to estimate the relative abundance of proteins in a sample. IBAQ value is the total intensity of a protein divided by the number of theoretical peptides.
+####  Compute IBAQ
 
 ```asciidoc
-python commands/peptides2proteins.py --fasta tests/PXD003947/Homo-sapiens-uniprot-reviewed-contaminants-decoy-202210.fasta --peptides tests/PXD003947/PXD003947-peptides-norm.csv --enzyme Trypsin --output tests/PXD003947/PXD003947-ibaq-norm.csv --normalize --verbose
+ibaqpy peptides2proteins --fasta tests/PXD003947/Homo-sapiens-uniprot-reviewed-contaminants-decoy-202210.fasta --peptides tests/PXD003947/PXD003947-peptides-norm.csv --enzyme Trypsin --output tests/PXD003947/PXD003947-ibaq-norm.csv --normalize --verbose
 ``` 
 
 The command provides an additional `flag` for normalize IBAQ values.
 
 ```asciidoc
-Usage: peptides2proteins.py [OPTIONS]
+Usage: peptides2proteins [OPTIONS]
 
 Options:
   -f, --fasta TEXT     Protein database to compute IBAQ values
@@ -174,29 +208,17 @@ Options:
   --help               Show this message and exit.
 ```
 
-#### 1. Performs the Enzymatic Digestion
-The current version of this tool uses OpenMS method to load fasta file, and use [ProteaseDigestion](https://openms.de/current_doxygen/html/classOpenMS_1_1ProteaseDigestion.html) to enzyme digestion of protein sequences, and finally get the theoretical peptide number of each protein.
-
-#### 2. Calculate the IBAQ Value
-First, peptide intensity dataframe was grouped according to protein name, sample name and condition. The protein intensity of each group was summed. Due to the experimental type, the same protein may exhibit missing peptides in different samples, resulting in variations in the number of peptides detected for the protein across different samples. To handle this difference, normalization within the same group can be achieved by using the formula `sum(peptides) / n`(n represents the number of detected peptide segments). Finally, the sum of the intensity of the protein is divided by the number of theoretical peptides.
-
-If protein-group exists in the peptide intensity dataframe, the intensity of all proteins in the protein-group is summed based on the above steps, and then multiplied by the number of proteins in the protein-group.
-
-#### 3. IBAQ Normalization  
-Normalize the ibaq values using the total ibaq of the sample. The resulted ibaq values are then multiplied by 100'000'000 (PRIDE database noramalization), for the ibaq ppb and log10 shifted by 10 (ProteomicsDB).
-
-### Compute TPA
-The total protein approach (TPA) is a label- and standard-free method for absolute protein quantitation of proteins using large-scale proteomic data. In the current version of the tool, the TPA value is the total intensity of the protein divided by its theoretical molecular mass.
+#### Compute TPA
 
-[ProteomicRuler](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4256500/) is a method for protein copy number and concentration estimation that does not require the use of isotope labeled standards. It uses the mass spectrum signal of histones as a "proteomic ruler" because it is proportional to the amount of DNA in the sample, which in turn depends on cell count. Thus, this approach can add an absolute scale to the mass spectrometry readout and allow estimates of the copy number of individual proteins in each cell.
+# TODO @PingZheng: Can you confirm that this command works?
 
 ```asciidoc
-python commands/compute_tpa.py --fasta Homo-sapiens-uniprot-reviewed-contaminants-decoy-202210.fasta --organism 'human' --peptides PXD003947-peptides.csv --ruler --ploidy 2 --cpc 200 --output PXD003947-tpa.tsv --verbose
+python compute_tpa --fasta Homo-sapiens-uniprot-reviewed-contaminants-decoy-202210.fasta --organism 'human' --peptides PXD003947-peptides.csv --ruler --ploidy 2 --cpc 200 --output PXD003947-tpa.tsv --verbose
 ```
 
 ```asciidoc
-python compute_tpa.py --help
-Usage: compute_tpa.py [OPTIONS]
+python compute_tpa --help
+Usage: compute_tpa [OPTIONS]
 
   Compute the protein copy numbers and concentrations according to a file output of peptides with the
   format described in peptide_normalization.py.
@@ -226,73 +248,17 @@ Options:
   --help                Show this message and exit.
 ```
 
-#### 1. Calculate the TPA Value
-The OpenMS tool was used to calculate the theoretical molecular mass of each protein. Similar to the calculation of IBAQ, the TPA value of protein-group was the sum of its intensity divided by the sum of the theoretical molecular mass.
+### Calculate the Cellular Protein Copy Number and Concentration
 
-#### 2. Calculate the Cellular Protein Copy Number and Concentration
 The protein copy calculation follows the following formula:
+
 ```
 protein copies per cell = protein MS-signal *  (avogadro / molecular mass) * (DNA mass / histone MS-signal)
 ```
-For cellular protein copy number calculation, the uniprot accession of histones were obtained from species first, and the molecular mass of DNA was calculated. Then the dataframe was grouped according to different conditions, and the copy number, molar number and mass of proteins were calculated.
-
-In the calculation of protein concentration, the volume is calculated according to the cell protein concentration first, and then the protein mass is divided by the volume to calculate the intracellular protein concentration.
-
-### Datasets Merge
-There are batch effects in protein identification and quantitative results between different studies, which may be caused by differences in experimental techniques, conditional methods, data analysis, etc. Here we provide a method to apply batch effect correction.  First to impute ibaq data, then remove outliers using `hdbscan`, and apply batch effect correction using `pycombat`.
-
-
-```asciidoc
-python commands/datasets_merge.py datasets_merge --data_folder ../ibaqpy_test/ --output datasets-merge.csv --verbose
-``` 
-
-```asciidoc
-python datasets_merge.py --help
-Usage: datasets_merge.py [OPTIONS]
-
-  Merge ibaq results from compute_ibaq.py.
-
-  :param data_folder: Data dolfer contains SDRFs and ibaq CSVs.
-  :param output: Output file after batch effect removal.
-  :param covariate: Indicators included in covariate consideration when datasets are merged.
-  :param organism: Organism to keep in input data.
-  :param covariate_to_keep: Keep tissue parts from metadata, e.g. 'LV,RV,LA,RA'.
-  :param non_missing_percent_to_keep: non-missing values in each group.
-  :param n_components: Number of principal components to be computed.
-  :param min_cluster: The minimum size of clusters.
-  :param min_sample_num: The minimum number of samples in a neighborhood for a point to be considered as a core point.
-  :param n_iter: Number of iterations to be performed.
-  :param verbose/quiet: Output debug information.
-
-Options:
-  Options:
-  -d, --data_folder TEXT          Data dolfer contains SDRFs and ibaq CSVs. [required]
-  -o, --output TEXT               Output file after batch effect removal. [required]
-  -c, --covariate TEXT            Indicators included in covariate consideration when datasets are merged.
-  --organism TEXT                 Organism to keep in input data.
-  --covariate_to_keep TEXT        Keep tissue parts from metadata, e.g. 'LV,RV,LA,RA'.
-  --non_missing_percent_to_keep FLOAT
-                                  non-missing values in each group.
-  --n_components TEXT             Number of principal components to be computed.
-  --min_cluster TEXT              The minimum size of clusters.
-  --min_sample_num TEXT           The minimum number of samples in a neighborhood for a point to be considered as a core point.
-  --n_iter TEXT                   Number of iterations to be performed.
-  -v, --verbose / -q, --quiet     Output debug information.
-  --help                          Show this message and exit.
-```
-
-#### 1. Impute Missing Values
-Remove proteins missing more than 30% of all samples. Users can keep tissue parts interested, and data will be converted to a expression matrix (samples as columns, proteins as rows), then missing values will be imputed with `KNNImputer`. 
-
-#### 2. Remove Outliers
-When outliers are removed, multiple hierarchical clustering is performed using `hdbscan.HDBSCAN`, where outliers are labeled -1 in the PCA plot. When clustering is performed, the default cluster minimum value and the minimum number of neighbors of the core point are the minimum number of samples in all datasets.
-
-*Users can skip this step and do outliers removal manually.*
 
-#### 3. Batch Effect Correction
-Using `pycombat` for batch effect correction, and batch is set to `datasets` (refers specifically to PXD ids) and the covariate should be `tissue_part`.
+For cellular protein copy number calculation, the uniprot accession of histones was obtained from species first, and the molecular mass of DNA was calculated. Then the dataframe was grouped according to different conditions, and the copy number, molar number and mass of proteins were calculated. In the calculation of protein concentration, the volume is calculated according to the cell protein concentration first, and then the protein mass is divided by the volume to calculate the intracellular protein concentration.
 
-### Citation
+### How to cite ibaqpy
 
 Wang H, Dai C, Pfeuffer J, Sachsenberg T, Sanchez A, Bai M, Perez-Riverol Y. Tissue-based absolute quantification using large-scale TMT and LFQ experiments. Proteomics. 2023 Oct;23(20):e2300188. doi: [10.1002/pmic.202300188](https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/10.1002/pmic.202300188). Epub 2023 Jul 24. PMID: 37488995.
 
diff --git a/commands/__init__.py b/commands/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/ibaq/__init__.py b/ibaqpy/__init__.py
similarity index 100%
rename from ibaq/__init__.py
rename to ibaqpy/__init__.py
diff --git a/__init__.py b/ibaqpy/commands/__init__.py
similarity index 100%
rename from __init__.py
rename to ibaqpy/commands/__init__.py
diff --git a/commands/compute_tpa.py b/ibaqpy/commands/compute_tpa.py
similarity index 92%
rename from commands/compute_tpa.py
rename to ibaqpy/commands/compute_tpa.py
index e6a34c1..00652ec 100644
--- a/commands/compute_tpa.py
+++ b/ibaqpy/commands/compute_tpa.py
@@ -1,22 +1,17 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
 import os
-
 import click
 import matplotlib.pyplot as plt
 import pandas as pd
 from matplotlib.backends.backend_pdf import PdfPages
 from pyopenms import *
 
-from ibaq.ibaqpy_commons import (
+from ibaqpy.ibaq.ibaqpy_commons import (
     CONDITION,
     NORM_INTENSITY,
     PROTEIN_NAME,
     SAMPLE_ID,
     plot_box_plot,
     plot_distributions,
-    print_help_msg,
     get_accession,
 )
 
@@ -33,12 +28,20 @@ def handle_nonstandard_aa(aa_seq: str) -> (list, str):
     return nonstandard_aa_lst, considered_seq
 
 
-@click.command()
-@click.option("-f", "--fasta", help="Protein database")
+@click.command("tpa", short_help="Compute TPA values.")
+@click.option(
+    "-f",
+    "--fasta",
+    help="Protein database",
+    required=True,
+    type=click.Path(exists=True),
+)
 @click.option(
     "-p",
     "--peptides",
     help="Peptide identifications with intensities following the peptide intensity output",
+    required=True,
+    type=click.Path(exists=True),
 )
 @click.option("-r", "--ruler", help="Whether to use ProteomicRuler", is_flag=True)
 @click.option("-n", "--ploidy", help="Ploidy number", default=2)
@@ -48,7 +51,7 @@ def handle_nonstandard_aa(aa_seq: str) -> (list, str):
 @click.option(
     "--verbose",
     help="Print addition information about the distributions of the intensities, number of peptides remove "
-    "after normalization, etc.",
+         "after normalization, etc.",
     is_flag=True,
 )
 @click.option(
@@ -57,15 +60,15 @@ def handle_nonstandard_aa(aa_seq: str) -> (list, str):
     default="TPA-QCprofile.pdf",
 )
 def tpa_compute(
-    fasta: str,
-    peptides: str,
-    ruler: bool,
-    organism: str,
-    ploidy: int,
-    cpc: float,
-    output: str,
-    verbose: bool,
-    qc_report: str,
+        fasta: str,
+        peptides: str,
+        ruler: bool,
+        organism: str,
+        ploidy: int,
+        cpc: float,
+        output: str,
+        verbose: bool,
+        qc_report: str,
 ) -> None:
     """
     This command computes the protein copies and concentrations according to a file output of peptides with the
@@ -81,9 +84,6 @@ def tpa_compute(
     :param qc_report: PDF file to store multiple QC images.
     :return:
     """
-    if peptides is None or fasta is None:
-        print_help_msg(tpa_compute)
-        exit(1)
 
     data = pd.read_csv(
         peptides, sep=",", usecols=[PROTEIN_NAME, NORM_INTENSITY, SAMPLE_ID, CONDITION]
@@ -100,7 +100,7 @@ def tpa_compute(
     proteins = res[PROTEIN_NAME].unique().tolist()
     proteins = sum([i.split(";") for i in proteins], [])
 
-    # calculate molecular weight of quantified proteins
+    # calculate the molecular weight of quantified proteins
     mw_dict = dict()
     fasta_proteins = list()  # type: list[FASTAEntry]
     FASTAFile().load(fasta, fasta_proteins)
@@ -115,7 +115,7 @@ def tpa_compute(
                 mw = AASequence().fromString(seq).getMonoWeight()
                 mw_dict[accession] = mw
                 print(
-                    f"Nonstandard amimo acids found in {accession}: {error_aa}, ignored!"
+                    f"Nonstandard amino acids found in {accession}: {error_aa}, ignored!"
                 )
 
     res = res[res[PROTEIN_NAME].isin(mw_dict.keys())]
@@ -237,7 +237,3 @@ def proteomic_ruler(df):
             pdf.savefig(box)
             pdf.close()
         res.to_csv(output, index=False)
-
-
-if __name__ == "__main__":
-    tpa_compute()
diff --git a/commands/datasets_merger.py b/ibaqpy/commands/datasets_merger.py
similarity index 72%
rename from commands/datasets_merger.py
rename to ibaqpy/commands/datasets_merger.py
index 6954dd6..613725f 100644
--- a/commands/datasets_merger.py
+++ b/ibaqpy/commands/datasets_merger.py
@@ -1,24 +1,9 @@
-#!/usr/bin/env python
 import re
 
 import click
+from ibaqpy.ibaq.combiner import Combiner
 
-from ibaq import __version__
-from ibaq.combiner import Combiner
 
-CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"])
-
-
-@click.group(context_settings=CONTEXT_SETTINGS)
-def cli():
-    """
-    This is the main tool that gives access to all commands.
-    """
-
-
-@click.version_option(
-    version=__version__, package_name="ibaqpy", message="%(package)s %(version)s"
-)
 @click.command("datasets_merge", short_help="Merge ibaq results from compute_ibaq")
 @click.option(
     "--data_folder",
@@ -75,19 +60,19 @@ def cli():
 )
 @click.pass_context
 def datasets_merge(
-    ctx,
-    data_folder: str,
-    output: str,
-    covariate: str,
-    organism: str,
-    covariate_to_keep: list,
-    non_missing_percent_to_keep: float,
-    skip_outliers_removal: bool,
-    n_components: int,
-    min_cluster: int,
-    min_sample_num: int,
-    n_iter: int,
-    verbose: bool,
+        ctx,
+        data_folder: str,
+        output: str,
+        covariate: str,
+        organism: str,
+        covariate_to_keep: list,
+        non_missing_percent_to_keep: float,
+        skip_outliers_removal: bool,
+        n_components: int,
+        min_cluster: int,
+        min_sample_num: int,
+        n_iter: int,
+        verbose: bool,
 ):
     if covariate_to_keep:
         covariate_to_keep = re.split(",\s*", covariate_to_keep)
@@ -98,14 +83,3 @@ def datasets_merge(
     combiner.batch_correction(n_components, covariate_to_keep)
     result = combiner.df_corrected
     result.to_csv(output, sep=",", index=True)
-
-
-cli.add_command(datasets_merge)
-
-
-def main():
-    cli()
-
-
-if __name__ == "__main__":
-    main()
diff --git a/commands/features2peptides.py b/ibaqpy/commands/features2peptides.py
similarity index 51%
rename from commands/features2peptides.py
rename to ibaqpy/commands/features2peptides.py
index de2ba2f..f20cc31 100644
--- a/commands/features2peptides.py
+++ b/ibaqpy/commands/features2peptides.py
@@ -1,13 +1,15 @@
-#!/usr/bin/env python
+
 import click
-from ibaq.peptide_normalization import peptide_normalization
+from ibaqpy.ibaq.peptide_normalization import peptide_normalization
 
 
-@click.command()
+@click.command("features2peptides", short_help="Convert features to parquet file.")
 @click.option(
     "-p",
     "--parquet",
     help="Parquet file import generated by quantms.io",
+    required=True,
+    type=click.Path(exists=True),
 )
 @click.option(
     "-s", "--sdrf", help="SDRF file import generated by quantms", default=None
@@ -63,7 +65,9 @@
     help="Save normalized peptides to parquet",
     is_flag=True,
 )
+@click.pass_context
 def features2parquet(
+    ctx,
     parquet: str,
     sdrf: str,
     min_aa: int,
@@ -78,8 +82,36 @@ def features2parquet(
     log2: bool,
     save_parquet: bool,
 ) -> None:
-    peptide_normalization(**click.get_current_context().params)
-
+    """
+    Convert features to parquet file.
+    :param ctx: Click context
+    :param parquet: Parquet file import generated by quantms.io
+    :param sdrf: SDRF file import generated by quantms
+    :param min_aa: Minimum number of amino acids to filter peptides
+    :param min_unique: Minimum number of unique peptides to filter proteins
+    :param remove_ids: Remove specific protein ids from the analysis using a file with one id per line
+    :param remove_decoy_contaminants: Remove decoy and contaminants proteins from the analysis
+    :param remove_low_frequency_peptides: Remove peptides that are present in less than 20% of the samples
+    :param output: Output file including all properties for normalization
+    :param skip_normalization: Skip normalization step
+    :param nmethod: Normalization method used to normalize feature intensities for tec (options: mean, median, iqr, none)
+    :param pnmethod: Peptide normalization method used to normalize peptides intensities for all samples (options:globalMedian,conditionMedian)
+    :param log2: Log2 transformation of peptide intensity values before normalization
+    :param save_parquet: Save normalized peptides to parquet
+    """
 
-if __name__ == "__main__":
-    features2parquet()
+    peptide_normalization(
+        parquet=parquet,
+        sdrf=sdrf,
+        min_aa=min_aa,
+        min_unique=min_unique,
+        remove_ids=remove_ids,
+        remove_decoy_contaminants=remove_decoy_contaminants,
+        remove_low_frequency_peptides=remove_low_frequency_peptides,
+        output=output,
+        skip_normalization=skip_normalization,
+        nmethod=nmethod,
+        pnmethod=pnmethod,
+        log2=log2,
+        save_parquet=save_parquet,
+    )
diff --git a/commands/merge_condition_files.py b/ibaqpy/commands/merge_condition_files.py
similarity index 75%
rename from commands/merge_condition_files.py
rename to ibaqpy/commands/merge_condition_files.py
index c0898c6..0454fb3 100644
--- a/commands/merge_condition_files.py
+++ b/ibaqpy/commands/merge_condition_files.py
@@ -1,20 +1,9 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
 
 import os
+import click
+import pandas as pd
 
-from ibaq.ibaqpy_commons import *
-
-
-def print_help_msg(command) -> None:
-    """
-    Print help information
-    :param command: command to print helps
-    :return: print help
-    """
-    with click.Context(command) as ctx:
-        click.echo(command.get_help(ctx))
+from ibaqpy.ibaq.ibaqpy_commons import CONDITION
 
 
 @click.command()
@@ -44,7 +33,3 @@ def merge_condition_generation(input: str, output: str, pattern: str) -> None:
         g.to_csv(
             f"{output}/{k}-grouped-Intensities.csv", index=False
         )  # '{}.csv'.format(k)
-
-
-if __name__ == "__main__":
-    merge_condition_generation()
diff --git a/commands/peptides2proteins.py b/ibaqpy/commands/peptides2proteins.py
similarity index 74%
rename from commands/peptides2proteins.py
rename to ibaqpy/commands/peptides2proteins.py
index 23737bc..58fa2af 100644
--- a/commands/peptides2proteins.py
+++ b/ibaqpy/commands/peptides2proteins.py
@@ -1,15 +1,22 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
+
 import click
-from ibaq.compute_ibaq import ibaq_compute
+from ibaqpy.ibaq.compute_ibaq import ibaq_compute
 
 
-@click.command()
-@click.option("-f", "--fasta", help="Protein database to compute IBAQ values")
+@click.command("peptides2proteins", short_help="Compute IBAQ values for proteins")
+@click.option(
+    "-f",
+    "--fasta",
+    help="Protein database to compute IBAQ values",
+    required=True,
+    type=click.Path(exists=True),
+)
 @click.option(
     "-p",
     "--peptides",
     help="Peptide identifications with intensities following the peptide intensity output",
+    required=True,
+    type=click.Path(exists=True),
 )
 @click.option(
     "-e",
@@ -41,7 +48,9 @@
     help="PDF file to store multiple QC images",
     default="IBAQ-QCprofile.pdf",
 )
+@click.pass_context
 def peptides2proteins(
+    click_context,
     fasta: str,
     peptides: str,
     enzyme: str,
@@ -55,6 +64,7 @@ def peptides2proteins(
     """
     This command computes the IBAQ values for a file output of peptides with the format described in
     peptide_contaminants_file_generation.py.
+    :param click_context: Click context
     :param min_aa: Minimum number of amino acids to consider a peptide.
     :param max_aa: Maximum number of amino acids to consider a peptide.
     :param fasta: Fasta file used to perform the peptide identification.
@@ -66,8 +76,14 @@ def peptides2proteins(
     :param qc_report: PDF file to store multiple QC images.
     :return:
     """
-    ibaq_compute(**click.get_current_context().params)
-
-
-if __name__ == "__main__":
-    peptides2proteins()
+    ibaq_compute(
+        fasta=fasta,
+        peptides=peptides,
+        enzyme=enzyme,
+        normalize=normalize,
+        min_aa=min_aa,
+        max_aa=max_aa,
+        output=output,
+        verbose=verbose,
+        qc_report=qc_report,
+    )
diff --git a/commands/tsne_visualization.py b/ibaqpy/commands/tsne_visualization.py
similarity index 96%
rename from commands/tsne_visualization.py
rename to ibaqpy/commands/tsne_visualization.py
index e787965..2fdf5fc 100644
--- a/commands/tsne_visualization.py
+++ b/ibaqpy/commands/tsne_visualization.py
@@ -1,6 +1,4 @@
-# import libraries
 import glob
-import math
 
 import click
 import matplotlib.pyplot as plt
@@ -10,7 +8,7 @@
 from sklearn.decomposition import PCA
 from sklearn.manifold import TSNE
 
-from ibaq.ibaqpy_commons import IBAQ_LOG, IBAQ_NORMALIZED, PROTEIN_NAME, SAMPLE_ID
+from ibaqpy.ibaq.ibaqpy_commons import PROTEIN_NAME, SAMPLE_ID, IBAQ_LOG
 
 
 # function to compute principal components
@@ -134,6 +132,7 @@ def plot_tsne(df, x_col, y_col, hue_col, file_name):
     "-o",
     "--pattern",
     help="Protein file pattern",
+    # TODO: I think we should use instead of pattern the structure of quantms.io for absolute quantification
     required=False,
     default="proteins.tsv",
 )
diff --git a/ibaqpy/ibaq/__init__.py b/ibaqpy/ibaq/__init__.py
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/ibaqpy/ibaq/__init__.py
@@ -0,0 +1,1 @@
+
diff --git a/ibaq/combiner.py b/ibaqpy/ibaq/combiner.py
similarity index 97%
rename from ibaq/combiner.py
rename to ibaqpy/ibaq/combiner.py
index e996dfd..0e4fcd4 100644
--- a/ibaq/combiner.py
+++ b/ibaqpy/ibaq/combiner.py
@@ -5,8 +5,8 @@
 import numpy as np
 import pandas as pd
 
-from ibaq.ibaqpy_commons import load_feature, load_sdrf
-from ibaq.utils import (
+from ibaqpy.ibaq.ibaqpy_commons import load_feature, load_sdrf
+from ibaqpy.ibaq.utils import (
     apply_batch_correction,
     compute_pca,
     fill_samples,
@@ -32,6 +32,16 @@ def __init__(
         self, data_folder: os.PathLike, covariate: str = None, organism: str = "HUMAN"
     ):
         """Generate concated IbaqNorm and metadata."""
+        self.df_pca = compute_pca(self.df_corrected.T, n_components=5)
+        self.df_corrected = None
+        self.batch_index = get_batch_info_from_sample_names(self.df.columns.tolist())
+        self.df_pca = None
+        self.df_filtered_outliers = None
+        self.batch_index = get_batch_info_from_sample_names(self.df.columns)
+        self.samples_number = None
+        self.datasets = None
+        self.samples = self.df.columns.tolist()
+        self.proteins = self.df["ProteinName"].unique().tolist()
         logger.info("Combining SDRFs and ibaq results ...")
         self.data_folder = Path(data_folder)
         if not self.data_folder.exists() or not self.data_folder.is_dir():
@@ -83,7 +93,6 @@ def imputer(self, covariate_to_keep: list = None):
         self.df = filter_missing_value_by_group(
             self.df, col="ProteinName", non_missing_percent_to_keep=0.3
         )
-        self.proteins = self.df["ProteinName"].unique().tolist()
 
         # TODO: Data for imputation should take samples as columns, proteins as rows. [Expression Matrix]
         # Also need to fill the proteins didn't show in original results for each sample.
@@ -101,7 +110,6 @@ def imputer(self, covariate_to_keep: list = None):
             self.df = fill_samples(self.df, self.proteins)
             self.df = impute_missing_values(self.df)
 
-        self.samples = self.df.columns.tolist()
         self.datasets = list(set([sample.split("-")[0] for sample in self.samples]))
         print(self.df.head)
 
@@ -122,7 +130,6 @@ def outlier_removal(
         min_samples = round(np.median(list(self.samples_number.values())))
         if min_samples == 1:
             min_samples = 2
-        self.batch_index = get_batch_info_from_sample_names(self.df.columns)
         # apply iterative outlier removal
         self.df_filtered_outliers = iterative_outlier_removal(
             self.df,
@@ -196,7 +203,6 @@ def batch_correction(
         # Before apply batch correction, filter out batches with just one sample (otherwise the batch correction will fail).
         batch_index = get_batch_info_from_sample_names(self.df.columns.tolist())
         self.df = remove_single_sample_batches(self.df, batch_index)
-        self.batch_index = get_batch_info_from_sample_names(self.df.columns.tolist())
 
         # get covariate information from metadata.
         columns = self.df.columns.tolist()
@@ -220,7 +226,6 @@ def batch_correction(
 
         # plot PCA of corrected data
         # transpose the dataframe to get samples as rows and features as columns
-        self.df_pca = compute_pca(self.df_corrected.T, n_components=5)
         # add batch information to the dataframe
         self.df_pca["batch"] = self.df_pca.index.str.split("-").str[0]
 
diff --git a/ibaq/compute_ibaq.py b/ibaqpy/ibaq/compute_ibaq.py
similarity index 97%
rename from ibaq/compute_ibaq.py
rename to ibaqpy/ibaq/compute_ibaq.py
index 6b3eb15..048b268 100644
--- a/ibaq/compute_ibaq.py
+++ b/ibaqpy/ibaq/compute_ibaq.py
@@ -1,15 +1,11 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
 import math
-
 import matplotlib.pyplot as plt
 import pandas as pd
 from matplotlib.backends.backend_pdf import PdfPages
 from pandas import DataFrame, Series
 from pyopenms import *
 
-from ibaq.ibaqpy_commons import (
+from ibaqpy.ibaq.ibaqpy_commons import (
     CONDITION,
     IBAQ,
     IBAQ_LOG,
@@ -20,11 +16,9 @@
     SAMPLE_ID,
     plot_box_plot,
     plot_distributions,
-    print_help_msg,
     get_accession,
 )
 
-
 def normalize(group):
     group[IBAQ_NORMALIZED] = group[IBAQ] / group[IBAQ].sum()
     return group
@@ -78,8 +72,7 @@ def ibaq_compute(
     :return:
     """
     if peptides is None or fasta is None:
-        print_help_msg(ibaq_compute)
-        exit(1)
+        raise ValueError("Fasta and peptides files are required")
 
     fasta_proteins = list()  # type: list[FASTAEntry]
     protein_accessions = list()
diff --git a/ibaq/ibaqpy_commons.py b/ibaqpy/ibaq/ibaqpy_commons.py
similarity index 96%
rename from ibaq/ibaqpy_commons.py
rename to ibaqpy/ibaq/ibaqpy_commons.py
index f75b9d5..2dd2337 100644
--- a/ibaq/ibaqpy_commons.py
+++ b/ibaqpy/ibaq/ibaqpy_commons.py
@@ -1,6 +1,4 @@
-#!/usr/bin/env python
 import os
-import click
 import matplotlib
 import numpy as np
 import pandas as pd
@@ -139,16 +137,6 @@
 }
 
 
-def print_help_msg(command: click.Command):
-    """
-    Print the help of the command
-    :param command: click command object
-    :return: None
-    """
-    with click.Context(command) as ctx:
-        click.echo(command.get_help(ctx))
-
-
 def get_accession(identifier: str) -> str:
     """
     Get protein accession from the identifier  (e.g. sp|P12345|PROT_NAME)
diff --git a/ibaq/normalization_methods.py b/ibaqpy/ibaq/normalization_methods.py
similarity index 97%
rename from ibaq/normalization_methods.py
rename to ibaqpy/ibaq/normalization_methods.py
index fb316e6..4c9aaef 100644
--- a/ibaq/normalization_methods.py
+++ b/ibaqpy/ibaq/normalization_methods.py
@@ -1,6 +1,6 @@
 import pandas as pd
 from sklearn.preprocessing import quantile_transform
-from ibaq.ibaqpy_commons import SAMPLE_ID, NORM_INTENSITY, TECHREPLICATE
+from ibaqpy.ibaq.ibaqpy_commons import SAMPLE_ID, NORM_INTENSITY, TECHREPLICATE
 
 
 def normalize_run(df, reps, method):
diff --git a/ibaq/peptide_normalization.py b/ibaqpy/ibaq/peptide_normalization.py
similarity index 90%
rename from ibaq/peptide_normalization.py
rename to ibaqpy/ibaq/peptide_normalization.py
index d2d446d..bceb268 100644
--- a/ibaq/peptide_normalization.py
+++ b/ibaqpy/ibaq/peptide_normalization.py
@@ -1,11 +1,11 @@
-#!/usr/bin/env python
 import pandas as pd
 import os
 import re
 import numpy as np
 import duckdb
-from ibaq.normalization_methods import normalize_run, normalize
-from ibaq.ibaqpy_commons import (
+from ibaqpy.ibaq.normalization_methods import normalize_run
+
+from ibaqpy.ibaq.ibaqpy_commons import (
     BIOREPLICATE,
     TECHREPLICATE,
     CHANNEL,
@@ -27,7 +27,6 @@
     ITRAQ4plex,
     ITRAQ8plex,
     parquet_map,
-    print_help_msg,
 )
 
 
@@ -67,7 +66,6 @@ def analyse_sdrf(sdrf_path: str) -> tuple:
     4. choice: A dictionary caontains key-values between channel
         names and numbers.
     :param sdrf_path: File path of SDRF.
-    :param compression: Whether compressed.
     :return:
     """
     sdrf_df = pd.read_csv(sdrf_path, sep="\t")
@@ -107,12 +105,12 @@ def get_label(labels: list) -> (str, dict):
         label = "LFQ"
     elif "TMT" in ",".join(labels) or "tmt" in ",".join(labels):
         if (
-            len(labels) > 11
-            or "TMT134N" in labels
-            or "TMT133C" in labels
-            or "TMT133N" in labels
-            or "TMT132C" in labels
-            or "TMT132N" in labels
+                len(labels) > 11
+                or "TMT134N" in labels
+                or "TMT133C" in labels
+                or "TMT133N" in labels
+                or "TMT132C" in labels
+                or "TMT132N" in labels
         ):
             choice = TMT16plex
         elif len(labels) == 11 or "TMT131C" in labels:
@@ -134,13 +132,12 @@ def get_label(labels: list) -> (str, dict):
 
 
 def remove_contaminants_entrapments_decoys(
-    dataset: pd.DataFrame, protein_field=PROTEIN_NAME
+        dataset: pd.DataFrame, protein_field=PROTEIN_NAME
 ) -> pd.DataFrame:
     """
     This method reads a file with a list of contaminants and high abudant proteins and
     remove them from the dataset.
     :param dataset: Peptide intensity DataFrame
-    :param contaminants_file: contaminants file
     :param protein_field: protein field
     :return: dataset with the filtered proteins
     """
@@ -153,7 +150,7 @@ def remove_contaminants_entrapments_decoys(
 
 
 def remove_protein_by_ids(
-    dataset: pd.DataFrame, protein_file: str, protein_field=PROTEIN_NAME
+        dataset: pd.DataFrame, protein_file: str, protein_field=PROTEIN_NAME
 ) -> pd.DataFrame:
     """
     This method reads a file with a list of contaminants and high abudant proteins and
@@ -171,13 +168,16 @@ def remove_protein_by_ids(
 
 
 def parquet_common_process(
-    data_df: pd.DataFrame, label: str, choice: dict
+        data_df: pd.DataFrame, label: str, choice: dict
 ) -> pd.DataFrame:
-    """Apply common process on data.
-
+    """
+    Apply a common process on data.
     :param data_df: Feature data in dataframe.
+    :param label: Label type of the experiment.
+    :param choice: Choice dict for a label type.
     :return: Processed data.
     """
+
     data_df = data_df.rename(columns=parquet_map)
     data_df[PROTEIN_NAME] = data_df.apply(lambda x: ";".join(x[PROTEIN_NAME]), axis=1)
     if label == "LFQ":
@@ -247,7 +247,7 @@ def merge_fractions(dataset: pd.DataFrame) -> pd.DataFrame:
 
 
 def get_peptidoform_normalize_intensities(
-    dataset: pd.DataFrame, higher_intensity: bool = True
+        dataset: pd.DataFrame, higher_intensity: bool = True
 ) -> pd.DataFrame:
     """
     Select the best peptidoform for the same sample and the same replicates. A peptidoform is the combination of
@@ -366,17 +366,17 @@ def get_label(labels: list) -> (str, dict):
         """
         choice = None
         if len(labels) == 1 and (
-            "LABEL FREE" in ",".join(labels) or "label free" in ",".join(labels)
+                "LABEL FREE" in ",".join(labels) or "label free" in ",".join(labels)
         ):
             label = "LFQ"
         elif "TMT" in ",".join(labels) or "tmt" in ",".join(labels):
             if (
-                len(labels) > 11
-                or "TMT134N" in labels
-                or "TMT133C" in labels
-                or "TMT133N" in labels
-                or "TMT132C" in labels
-                or "TMT132N" in labels
+                    len(labels) > 11
+                    or "TMT134N" in labels
+                    or "TMT133C" in labels
+                    or "TMT133N" in labels
+                    or "TMT132C" in labels
+                    or "TMT132N" in labels
             ):
                 choice = TMT16plex
             elif len(labels) == 11 or "TMT131C" in labels:
@@ -399,7 +399,8 @@ def get_label(labels: list) -> (str, dict):
     def get_report_from_database(self, samples: list, columns: list = None):
         """
         This function loads the report from the duckdb database for a group of ms_runs.
-        :param runs: A list of ms_runs
+        :param columns: A list of columns
+        :param samples: A list of samples
         :return: The report
         """
         cols = ','.join(columns) if columns is not None else '*'
@@ -413,11 +414,11 @@ def get_report_from_database(self, samples: list, columns: list = None):
 
     def iter_samples(self, file_num: int = 20, columns: list = None):
         """
-        :params file_num: The number of files being processed at the same time(default 20)
+        :params file_num: The number of files being processed at the same time (default 20)
         :yield: _description_
         """
         ref_list = [
-            self.samples[i : i + file_num]
+            self.samples[i: i + file_num]
             for i in range(0, len(self.samples), file_num)
         ]
         for refs in ref_list:
@@ -471,7 +472,8 @@ def get_median_map(self):
     def get_report_condition_from_database(self, cons: list, columns: list = None):
         """
         This function loads the report from the duckdb database for a group of ms_runs.
-        :param runs: A list of ms_runs
+        :param columns: A list of columns
+        :param cons: A list of conditions in
         :return: The report
         """
         cols = ','.join(columns) if columns is not None else '*'
@@ -486,7 +488,7 @@ def get_report_condition_from_database(self, cons: list, columns: list = None):
     def iter_conditions(self, conditions: int = 10, columns: list = None):
         condition_list = self.get_unique_conditions()
         ref_list = [
-            condition_list[i : i + conditions]
+            condition_list[i: i + conditions]
             for i in range(0, len(condition_list), conditions)
         ]
         for refs in ref_list:
@@ -503,7 +505,7 @@ def get_unique_conditions(self):
     def get_median_map_to_condition(self):
         med_map = {}
         for cons, batch_df in self.iter_conditions(
-            1000, ["condition", "sample_accession", "intensity"]
+                1000, ["condition", "sample_accession", "intensity"]
         ):
             for con in cons:
                 meds = (
@@ -517,26 +519,25 @@ def get_median_map_to_condition(self):
 
 
 def peptide_normalization(
-    parquet: str,
-    sdrf: str,
-    min_aa: int,
-    min_unique: int,
-    remove_ids: str,
-    remove_decoy_contaminants: bool,
-    remove_low_frequency_peptides: bool,
-    output: str,
-    skip_normalization: bool,
-    nmethod: str,
-    pnmethod: str,
-    log2: bool,
-    save_parquet: bool,
+        parquet: str,
+        sdrf: str,
+        min_aa: int,
+        min_unique: int,
+        remove_ids: str,
+        remove_decoy_contaminants: bool,
+        remove_low_frequency_peptides: bool,
+        output: str,
+        skip_normalization: bool,
+        nmethod: str,
+        pnmethod: str,
+        log2: bool,
+        save_parquet: bool,
 ) -> None:
     if os.path.exists(output):
-        exit(f"{output} already exist!")
+        raise FileExistsError("The output file already exists.")
 
     if parquet is None:
-        print_help_msg(peptide_normalization)
-        exit(1)
+        raise FileNotFoundError("The file does not exist.")
 
     print("Loading data..")
     F = Feature(parquet)
@@ -559,12 +560,12 @@ def peptide_normalization(
             if not skip_normalization:
                 if pnmethod == "globalMedian":
                     dataset_df.loc[:, "intensity"] = (
-                        dataset_df["intensity"] / med_map[sample]
+                            dataset_df["intensity"] / med_map[sample]
                     )
                 elif pnmethod == "conditionMedian":
                     con = dataset_df["condition"].unique()[0]
                     dataset_df.loc[:, "intensity"] = (
-                        dataset_df["intensity"] / med_map[con][sample]
+                            dataset_df["intensity"] / med_map[con][sample]
                     )
             dataset_df = dataset_df[PARQUET_COLUMNS]
             dataset_df = parquet_common_process(dataset_df, label, choice)
@@ -574,8 +575,10 @@ def peptide_normalization(
                 lambda x: len(set(x[PEPTIDE_CANONICAL])) >= min_unique
             )
             dataset_df.rename(columns={INTENSITY: NORM_INTENSITY}, inplace=True)
+            ## TODO: @PingZheng I think the removal of these features must happen before the normalization step.
+            ## TODO: The entire point to catch any outliers and remove them before normalization.
 
-            # Remove high abundant, entrapments, contaminants proteins and the outliers
+            # Remove high abundant entrapment's, contaminants, proteins and the outliers
             if remove_ids is not None:
                 dataset_df = remove_protein_by_ids(dataset_df, remove_ids)
             if remove_decoy_contaminants:
@@ -601,9 +604,9 @@ def peptide_normalization(
 
             # TODO: Normalize at feature level between ms runs (technical repetitions)
             if (
-                not skip_normalization
-                and nmethod != "none"
-                and technical_repetitions > 1
+                    not skip_normalization
+                    and nmethod != "none"
+                    and technical_repetitions > 1
             ):
                 print(f"{str(sample).upper()}: Normalize intensities of features.. ")
                 dataset_df = normalize_run(dataset_df, technical_repetitions, nmethod)
diff --git a/ibaq/utils.py b/ibaqpy/ibaq/utils.py
similarity index 99%
rename from ibaq/utils.py
rename to ibaqpy/ibaq/utils.py
index 3f4416b..310c1ab 100644
--- a/ibaq/utils.py
+++ b/ibaqpy/ibaq/utils.py
@@ -6,10 +6,10 @@
 import matplotlib.pyplot as plt
 import pandas as pd
 import seaborn as sns
+from sklearn.cluster._hdbscan import hdbscan
 from sklearn.decomposition import PCA
 from sklearn.impute import KNNImputer
 
-import hdbscan
 from combat.pycombat import pycombat
 
 logging.basicConfig(
diff --git a/ibaqpy/ibaqpyc.py b/ibaqpy/ibaqpyc.py
new file mode 100644
index 0000000..79c3f12
--- /dev/null
+++ b/ibaqpy/ibaqpyc.py
@@ -0,0 +1,44 @@
+import click
+
+from ibaqpy.commands.features2peptides import features2parquet
+from ibaqpy.commands.compute_tpa import tpa_compute
+from ibaqpy.commands.merge_condition_files import merge_condition_generation
+from ibaqpy.commands.peptides2proteins import peptides2proteins
+from ibaqpy.commands.tsne_visualization import tsne_visualization
+import ibaqpy.__init__ as __init__
+
+CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"])
+
+
+@click.group(context_settings=CONTEXT_SETTINGS)
+@click.version_option(
+    version=__init__.__version__,
+    package_name="ibaqpy",
+    message="%(package)s %(version)s",
+)
+def cli():
+    """
+    Main entry point for the CLI
+    """
+
+
+cli.add_command(tpa_compute)
+cli.add_command(features2parquet)
+cli.add_command(merge_condition_generation)
+cli.add_command(peptides2proteins)
+cli.add_command(tsne_visualization)
+
+
+def main():
+    """
+    Main function to run the CLI
+    """
+    try:
+        cli()
+    except SystemExit as e:
+        if e.code != 0:
+            raise
+
+
+if __name__ == "__main__":
+    main()
diff --git a/requirements.txt b/requirements.txt
index f8a921d..1ffbba3 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,4 +10,5 @@ qnorm
 scipy>=1.10
 seaborn
 typing_extensions==4.6.3
-setuptools==67.8.0
\ No newline at end of file
+setuptools==67.8.0
+combat
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 4012a9a..8007177 100644
--- a/setup.py
+++ b/setup.py
@@ -38,15 +38,9 @@ def readme():
         "qnorm",
         "seaborn",
         "typing_extensions",
+        "combat",
     ],
-    scripts=[
-        "commands/features2peptides.py",
-        "commands/peptides2proteins.py",
-        "commands/compute_tpa.py",
-        "commands/datasets_merger.py",
-        "commands/tsne_visualization.py",
-        "commands/merge_condition_files.py",
-    ],
+    entry_points={"console_scripts": ["ibaqpy=ibaqpy.ibaqpyc:main"]},
     classifiers=[
         "Development Status :: 4 - Beta",
         "Environment :: Console",


Test patch:
diff --git a/tests/PXD003947/PXD003947-featrue.parquet b/tests/PXD003947/PXD003947-feature.parquet
similarity index 100%
rename from tests/PXD003947/PXD003947-featrue.parquet
rename to tests/PXD003947/PXD003947-feature.parquet
diff --git a/tests/test_ibaqpy.py b/tests/test_ibaqpy.py
index 8788055..903101f 100644
--- a/tests/test_ibaqpy.py
+++ b/tests/test_ibaqpy.py
@@ -1,13 +1,13 @@
 from unittest import TestCase
 
-from ibaq.peptide_normalization import peptide_normalization
-from ibaq.compute_ibaq import ibaq_compute
+from ibaqpy.ibaq.peptide_normalization import peptide_normalization
+from ibaqpy.ibaq.compute_ibaq import ibaq_compute
 
 
 class TestIbaqpy(TestCase):
     def test_feature_assembly(self):
         args = {
-            "parquet": __package__ + "PXD003947/PXD003947-featrue.parquet",
+            "parquet": __package__ + "PXD003947/PXD003947-feature.parquet",
             "sdrf": __package__ + "PXD003947/PXD003947.sdrf.tsv",
             "min_aa": 7,
             "min_unique": 2,