# Source codes, installation guides, and step by step installing commands
Manuscript: Predicting Protein Synergistic Effect in Arabidopsis using Epigenome Profiling
Author: Chih-Hung Hsieh, Ya-Ting Sabrina Chang, Ming-Ren Yen, Jo-Wei Allison Hsieh, Pao-Yang Chen


# System Requirements


## 1_ChIPpeakcall

### Introduction
All Chip-Seq libraries in the database were processed by a unified python scripting pipeline. To get started, download the script below and input your CSV file with at least two columns: Sample Name, Sample Fastq, Control Name and Control Fastq (if any). Output includes two folders namely "peakcall" (narrowPeak files) and "BigWig" (bw files).

This script is used for the peak calling pipeline of histone ChIP-seq data.
We recommend using the provided [Docker images](https://hub.docker.com/r/dppss90008/qhistone-pipeline) to run our peak calling pipeline.

### Installation and Step by step running command

**STEP1**: Pull docker images from docker hub
```
docker pull dppss90008/qhistone-pipeline
```
**STEP2**: Download bowtie2 Index File
```
wget https://genome-idx.s3.amazonaws.com/bt/TAIR10.zip
unzip TAIR10.zip
```
**STEP3**: Create a CSV file
The file scheduled the samples to be processed
1.  The file name of the CSV file should be Peak-ProcessTable.csv  
      Example of [Peak-ProcessTable.csv](https://qhistone.paoyang.ipmb.sinica.edu.tw/static/download/Peak-ProcessTable.csv)
2.  If the sample have multiple SRR files, use ";" to delimite them.
3.  If the sample have no input control, left blank.

**STEP4**: Prepare the ChIP-seq data
-   The file scheduled the SRR file to be downloaded
-   The file name of the CSV file should be SRAfile.txt  
      Example of [SRAfile.txt](https://qhistone.paoyang.ipmb.sinica.edu.tw/static/download/SRAfile.txt)
```
docker run -it --rm \
  -v /path/of/SRAfiles/:/SRAFiles/ \
  dppss90008/qhistone-pipeline \
  prefetch --option-file /SRAFiles/SRAfile.txt --output-directory /SRAFiles
```
**STEP5**: Run the code
```
docker run -it --rm \
  -v /path/of/Peak-ProcessTable.csv:/source \
  -v /path/of/TAIR10-Bowtie2Index:/Bowtie2Index \
  -v /path/of/SRAfiles/:/SRAFiles/ \
  dppss90008/qhistone-pipeline \
  conda run -n chipseq python /PeakCallingPipeline/Peak-Calling-Pipeline-Version4.py \
    --cores 70 \ #  CPU cores for speeding up the pipeline
    --bowtie2index TAIR10 \ # bowtie2 index file downloaded from step 2
    --wkdir /source \
```

### Source data for running peak-calling pipeline

1. The metadata of histone related ChIP-seqs and 1,566 protein ChIP-seqs for preparing "Peak-ProcessTable.csv" (STEP3) and "SRAfile.txt" (STEP4) are ava
2. After running the peak-calling pipeline, we obtained results from 1,534 histone-related ChIP-seq experiments and 1,566 protein ChIP-seq experiments. The data is available at
https://zenodo.org/doi/10.5281/zenodo.11228366.

## 2_ChIPdatabase

### Installation

We suggest useing the lastest version of python.
Install the packages: 
```pip install pandas numpy scikit-learn joblib scipy billiard pybedtools matplotlib mplcursors seaborn```


### Introduction

"Predicting the epigenome profile of a protein query (phase I)"

In phase 1, we download and preprocess 1,534 ChIP-seq datasets of histone modifications/variants from NCBI, resulting in a total of 596,343 bins for the TAIR10 genome. We represent histone modification presence in the genome (i.e., peak locations) with a Boolean digit system, where "1" indicates a "marking" and "0" indicates "no marking." Consequently, a matrix is created to store all Boolean digits from the 1,534 histone 488 modification datasets (i.e., 1,534 ChIP-seq libraries x 598,343 bins). The final matrix of QHistone covered 1,534 libraries, significantly more than the sum of the other two major 490 databases of Arabidopsis histone modifications.

Scripts in 2_ChIPdatabase are used for preparing matrix for machine learning.

### Step by step running command

**STEP1**: Create 200 bp Binning TAIR10 Genome

Script: 0_make_genomic_bins.sh
Input: TAIR10.chromInfo
Output: TAIR10_200bp.bed
```
sh 0_make_genomic_bins.sh
```
**STEP2**: Flatten the narrowPeak data on 200bp genomic bins

Script: 1_Peak_200BP.py
Input: narrowPeaks in peakcall folder
(Peak calling results are available at
https://zenodo.org/doi/10.5281/zenodo.11228366.)
Output: 200BP/{sample}_200.id 
```
python 1_Peak_200BP.py
```
**STEP3**: Create matrix

Create matrixes of each peak data “1” represented "marking” whereas “0” represented “no marking”
Script: 2_Peak_Flat_200BP.py
Input: 200BP/{sample}_200.id ; TAIR10_200bp.bed
Output: 200BP_Genome/{sample}_200.id_Genome
```
python 2_Peak_Flat_200BP.py
```

**STEP4** Create a matrix npy and metadata metrix for machine learning

Script: 3_Prepare_MetaData.py
Input: QHistone_Histone_ChIP_metadata.csv、200BP_Genome/{sample}_200.id_Genome
Output: Genome.npy、Genome_metadata.csv
The metadata is available at https://zenodo.org/doi/10.5281/zenodo.11228366.
```
python 3_Prepare_MetaData.py
```
**ALTERNATIVE-STEP1**: Flatten the narrowPeak data on Gene/TE

Code for binning peaks (Gene and TE)

Script: 4_Peak_Flat_GeneTE.py
Input: TAIR10_TE_enhancer.bed、Genome_metadata.csv、peakcall/{sample}_peaks.narrowPeak
Output: Gene_Genome/{sample}_Gene、TE_Genome/{sample}_TE
```
python 4_Peak_Flat_GeneTE.py
```
**ALTERNATIVE-STEP2**: Create a matrix npy for machine learning

Code for binning peaks (Gene and TE)

Script: 5_Prepare_MetaData-Gene-TE.py
Input: Gene_Genome/{sample}_Gene、TE_Genome/{sample}_TE、Genome_metadata.csv
Output: Gene.npy、TE.npy
```
python 5_Prepare_MetaData-Gene-TE.py
```

## 3_ChIP_SVM

### Introduction
In phase 2, we implemented SVM and searched for the best SVM model with optimal parameters to cluster 1,534 ChIP-seq libraries. We used the machine learning Python package scikit-learn (version 0.23.2) to implement SVM, taking the matrix from phase 1 as the input. To determine the optimal parameters, including regularization parameter, kernel type, and kernel coefficient, we explored various parameter combinations using grid search coupled with 5-fold cross-validation. We assessed model performance by calculating overall accuracy, F1 score, recall, and precision. In the 5-fold cross-validation, the dataset of 1,534 preprocessed ChIP-seqs is randomly split into 5 groups. The model is trained on 4 groups in each cycle, while the fifth group, containing the actual answers, is used for performance evaluation. This cycle is repeated five times, with different groups used for testing each time, thereby providing a comprehensive evaluation of the model's accuracy, F1 score, recall, and precision across multiple trials. As a result, the linear SVM model with a regularization parameter C=0.1 demonstrated the best performance in the test datasets.

### Step by step running command

**STEP1**: GridSearch for finding best parameter with 5 fold CV
Script: 1_SVM_GridSearch.py
Input: Genome.npy (in 2_ChIPdatabase folder)、Genome_metadata.csv (in 2_ChIPdatabase folder)
Output: GridSummary.txt (result)、{kernel}_Grid_Accuracy.png (plot)
```
python 1_SVM_GridSearch.py
```
**STEP2**: ROC Analysis for 4 SVM kernals
Script: 2_SVM_ROC_4kernels.py
Input: Genome.npy (in 2_ChIPdatabase folder)、Genome_metadata.csv (in 2_ChIPdatabase folder)
Output: ROC_kernel_score.txt (accuracy for 4 kernals、{kernel}_{fpr/tpr/roc}.npy (result)、ROC4kernel.png (plot)
```
python 2_SVM_ROC_4kernels.py
```

**STEP3**: Final model for prediction
Script: 3_SVM_model.py
Input: Genome.npy (in 2_ChIPdatabase folder)、Genome_metadata.csv (in 2_ChIPdatabase folder)、QHistone_Histone_ChIP_metadata.csv (in 2_ChIPdatabase folder)
Output: svm_allregions_report_matrix_remove.txt (classification report)、ChIP_linear_C01_clf_remove.pkl (SVM model)、ChIP_linear_C01_clf_noH3_remove.pkl (SVM model, no H3 version)
```
# Note: Remove comment in code line 50 and 87 for no H3 version
python 3_SVM_model.py
```

**STEP4**: Recursive Feature Elimination (Find the import features for predictions)
Script: 4_ChIP_RFE.py
Input: 
  Genome.npy (in 2_ChIPdatabase folder)
  Genome_metadata.csv (in 2_ChIPdatabase folder)
  QHistone_Histone_ChIP_metadata.csv (in 2_ChIPdatabase folder)
output: 
  RFE_CV5_c01_1.txt (5foldCV)
  RFE_CV5_c01_1000.txt
  RFE_CV5_c01_10000.txt
  RFE_LOO_c01_1.txt (LOOCV)
  rfe_linear_c01_step1_RandomState42.npy (feature importance)
  rfe_linear_c01_step1000_RandomState42.npy (feature importance)
  rfe_linear_c01_step10000_RandomState42.npy (feature importance)
```
python 4_ChIP_RFE.py
```
**STEP5**: Recursive Feature Elimination (Make critical modol for prediction)
Script: 4_ChIP_RFE_model.py
Input: 
  Genome.npy (in 2_ChIPdatabase folder)
  Genome_metadata.csv (in 2_ChIPdatabase folder)
  QHistone_Histone_ChIP_metadata.csv (in 2_ChIPdatabase folder)
  rfe_linear_c01_step1_RandomState42.npy (feature importance)
  rfe_linear_c01_step1000_RandomState42.npy (feature importance)
  rfe_linear_c01_step10000_RandomState42.npy (feature importance)
Output: 
  ChIP_linear_C1_critical303_remove_clf.pkl
  ChIP_linear_C1_critical303_noH3_remove_clf.pkl

```
# Note: Remove comment in code line 53 and 143 for no H3 version
python 4_ChIP_RFE_model.py
```
**STEP6**: Final model for prediction (Gene TE version)
Script: 5_ChIP_GeneTE.py
Input: 
  TE.npy (in 2_ChIPdatabase folder)
  Gene.npy (in 2_ChIPdatabase folder)
  Genome_metadata.csv (in 2_ChIPdatabase folder)
  QHistone_Histone_ChIP_metadata.csv (in 2_ChIPdatabase folder)
Output: 
  ChIP_linear_C01_clf_TE_noH3_remove.pkl (SVM model, no H3 TE version)
  ChIP_linear_C01_clf_Gene_noH3_remove.pkl (SVM model, no H3 Gene version)
```
python 5_ChIP_GeneTE
```

**STEP7**: Enrichment analysis of critical regions
Sciprt: 6_EnrichmentAnalysis303.ipynb
Input: 
  TAIR10_TE_enhancer.bed
  TAIR10_200bp.bed
  rfe_linear_c01_step1_RandomState42.npy (feature importance)
  rfe_linear_c01_step1000_RandomState42.npy (feature importance)
  rfe_linear_c01_step10000_RandomState42.npy (feature importance)
Output: 
  critical_regions303.bed
  CriticalRegion303_enrichment.png 

```
Use jupyterNotebook to open 6_EnrichmentAnalysis303.ipynb and run the commands
```

**STEP8**: Plot RFE
Script: 6_RFE_plot.ipynb
Input: 
  RFE_CV5_c01_1.txt (5foldCV)
  RFE_CV5_c01_1000.txt
  RFE_CV5_c01_10000.txt
  RFE_LOO_c01_1.txt (LOOCV)
Output: 
  SVM-RFE (5-fold cv)_C01.png
  SVM-RFE_LOO_5FOLD.png
```
Use jupyterNotebook to open 6_RFE_plot.ipynb and run the commands
```

## 4_ChIP_tSNE_PCA_MDS
**STEP1**: ChIP tSNE_PCA_MDS Analysis and Plotting / TSNE_PCA_MDS_Analysis.py
Script: a. ChIP tSNE_PCA_MDS Analysis and Plotting.ipynb
```
Use jupyterNotebook to open a. ChIP tSNE_PCA_MDS Analysis and Plotting.ipynb and run the commands
```
**STEP2**: tSNE with new Query
Script: b. tSNE with new Query.ipynb
```
Use jupyterNotebook to open b. tSNE with new Query.ipynb and run the commands
```

## 5_ChIP_hi-cluster
Script: HiCluster.ipynb
```
Use jupyterNotebook to open HiCluster.ipynb and run the commands
```

## 6_ProteinChIPs

### Introduction
QHistone has developed epigenome profiles for 616 Arabidopsis proteins, with their ChIP-seq data accessible in the GEO database. To discover potential protein partners or functionally associated pairs within these libraries, we conducted 2.4M pairwise comparisons of their epigenome profiles. The outcomes are depicted based on their KLD scores, MI values, and the overlapping percentage of binding locations between two proteins to aid in evaluating the functional consistency of these pairs.

### Step by step running command

**STEP1**: Flatten the narrowPeak data on 200bp genomic bins

Script: 1_Peak_200BP.py
Input: narrowPeaks in peakcall folder
(Peak calling results are available at
https://zenodo.org/doi/10.5281/zenodo.11228366.)
Output: 200BP/{sample}_200.id 
```
python 1_Peak_200BP.py
```

**STEP2**: Create matrix

Create matrixes of each peak data “1” represented "binding” whereas “0” represented “no binding”
Script: 2_Peak_Flat_200BP.py
Input: 200BP/{sample}_200.id ; TAIR10_200bp.bed
Output: 200BP_Genome/{sample}_200.id_Genome
```
python 2_Peak_Flat_200BP.py
```

**STEP3**: Predict epigenome profile of all protein ChIPs
Script: 3_SVM_Result.py
Input: 200BP_Genome/{sample}_Genome
Ouput: SVM_result/allregion/predict_res_{sample}_200_allregion.csv
	   SVM_result/critical/predict_res_{sample}_200_critical.csv
```
python 3_SVM_Result.py
```
**STEP4**: Pairwise comparision of epigenome profile

Scirpt: 4_pairwise.py
Input: all_genomic_noH3 
Output: protein_ChIP_pairwise.csv
```
python 4_pairwise.py
```

**STEP5**: Plot Pairwise comparision of all protein ChIPs

Scirpt: 5_Plot_pairwise.py
Input: 
  protein_ChIP_pairwise.csv
  TF_TableManagement_hsieh_frip_version7.csv
Output: 
  MI-R_Ovp_1566x1565_hex.png
  KLD_Ovp_1566x1565_hex.png
```
python 5_Plot_pairwise.py
```
## 7_QhistoneWeb
Source codes of QHistone website. For better use our tools,
Please use "https://qhistone.paoyang.ipmb.sinica.edu.tw/".

## 8_Comprehensive_histone_Analysis

Script: 1. Histone enrichment analysis.ipynb 
```
Use jupyterNotebook to open 1. Histone enrichment analysis.ipynb and run the commands
```
Script: 2. Histone overlapping analysis.ipynb
```
Use jupyterNotebook to open 2. Histone overlapping analysis.ipynb and run the commands
```
Script: 3. MetaAnalysis of Histone mark.ipynb
```
Use jupyterNotebook to open 3. MetaAnalysis of Histone mark.ipynb and run the commands
```