diploshictest 0.2
DiploS/HIC
License: MIT

# diploS/HIC
This repo contains the implementation for `diploS/HIC` as described in Kern and Schrider (2018; https://doi.org/10.1534/g3.118.200262), along 
with its associated support scripts. `diploS/HIC` uses a deep convolutional neural network to identify
hard and soft selective sweep in population genomic data. 

The workflow for analysis using `diploS/HIC` consists of four basic parts. 1) Generation of a training set for `diploS/HIC` 
using simulation. 2) `diploS/HIC` training and performance evaluation. 3) Calculation of `dipoS/HIC` feature vectors from genomic data.
4) prediction on empirical data using the trained network. The software provided here can handle the last three parts; population
genetic simulations must be performed using separate software such as discoal (https://github.com/kern-lab/discoal) 

## Installation
`diploS/HIC` has a number of dependencies that should be straightforward to install using python package managers
such as `conda` or `pip`. The complete list of dependencies looks like this:

- numpy
- scipy
- pandas
- scikit-allel
- scikit-learn
- tensorflow
- keras

## Install on linux
I'm going to focus on the steps involved to install on a linux machine using Anaconda as our python source / main
package manager. First download and install Anaconda

```
$ wget https://repo.continuum.io/archive/Anaconda3-5.0.1-Linux-x86_64.sh
$ bash Anaconda3-5.0.1-Linux-x86_64.sh
```
That will give us the basics (numpy, scipy, scikit-learn, etc). Next lets install scikit-allel using `conda`
```
$ conda install -c conda-forge scikit-allel
```
That's easy. Installing tensorflow and keras can be slightly more touchy. You will need to determine if 
you want to use a CPU-only implementation (probably) or a GPU implementation of tensorflow. See
https://www.tensorflow.org/install/install_linux for install instructions. I'm going to install the
CPU version for simplicity. tensorflow and keras are the libraries which handle the deep learning
portion of things so it is important to make sure the versions of these libraries play well together 
```
$ pip install tensorflow 
$ pip install keras
```
Note that because I'm using the Anaconda version of python, pip will only install this in the anaconda directory
which is a good thing. Okay that should be the basics of dependencies. Now we are ready to install `diploS/HIC` itself
```
$ git clone https://github.com/kern-lab/diploSHIC.git
$ cd diploSHIC 
$ python setup.py install
```
Assuming all the dependencies were installed this should be all set

## Usage
The main program that you will interface with is `diploSHIC.py`. This script has four run modes that allow the user to 
perform each of the main steps in the supervised machine learning process. We will briefly lay out the modes of use
and then will provide a complete example of how to use the program for fun and profit.

`diploSHIC.py` uses the `argparse` module in python to try to give the user a complete, command line based help menu. 
We can see the top level of this help by typing
```
$ python diploSHIC.py -h
usage: diploSHIC.py [-h] {train,predict,fvecSim,fvecVcf} ...

calculate feature vectors, train, or predict with diploSHIC

possible modes (enter 'python diploSHIC.py modeName -h' for modeName's help message:
  {fvecSim,makeTrainingSets,train,fvecVcf,predict}
                        sub-command help
    fvecSim             Generate feature vectors from simulated data
    makeTrainingSets    Combine feature vectors from muliple fvecSim runs into
                        5 balanced training sets
    train               train and test a shic CNN
    fvecVcf             Generate feature vectors from data in a VCF file
    predict             perform prediction using an already-trained SHIC CNN

optional arguments:
  -h, --help            show this help message and exit
```
### before running diploSHIC: simulating training/testing data
All flavors of S/HIC require simulated data for training (and ideally, testing). Users can select whatever simulator 
they prefer and parameterize them however they wish. We have included an example script in this respository 
(generateSimLaunchScript.py) which demonstrates how a training set can be simulated with discoal (available at 
https://github.com/kern-lab/discoal).

### feature vector generation modes
The first task in our pipeline is generating feature vectors from simulation data (or empirical data) to
use with the CNN that we will train and then use for prediction. The `diploSHIC.py` script eases this 
process with two run modes

#### fvecSim mode
The fvecSim run mode is used for turning ms-style output into feature vectors compatible with `diploSHIC.py`. The
help message from this mode looks like this
```
$ python diploSHIC.py fvecSim -h
usage: diploSHIC.py fvecSim [-h] [--totalPhysLen TOTALPHYSLEN]
                            [--numSubWins NUMSUBWINS]
                            [--maskFileName MASKFILENAME]
                            [--chrArmsForMasking CHRARMSFORMASKING]
                            [--unmaskedFracCutoff UNMASKEDFRACCUTOFF]
                            [--outStatsDir OUTSTATSDIR]
                            [--ancFileName ANCFILENAME] [--pMisPol PMISPOL]
                            shicMode msOutFile fvecFileName

required arguments:
  shicMode              specifies whether to use original haploid SHIC (use
                        'haploid') or diploSHIC ('diploid')
  msOutFile             path to simulation output file (must be same format
                        used by Hudson's ms)
  fvecFileName          path to file where feature vectors will be written

optional arguments:
  -h, --help            show this help message and exit
   --totalPhysLen TOTALPHYSLEN
                        Length of simulated chromosome for converting infinite
                        sites ms output to finite sites (default=1100000)
  --numSubWins NUMSUBWINS
                        The number of subwindows that our chromosome will be
                        divided into (default=11)
  --maskFileName MASKFILENAME
                        Path to a fasta-formatted file that contains masking
                        information (marked by 'N'). If specified, simulations
                        will be masked in a manner mirroring windows drawn
                        from this file.
  --chrArmsForMasking CHRARMSFORMASKING
                        A comma-separated list (no spaces) of chromosome arms
                        from which we want to draw masking information (or
                        'all' if we want to use all arms. Ignored if
                        maskFileName is not specified.
  --unmaskedFracCutoff UNMASKEDFRACCUTOFF
                        Minimum fraction of unmasked sites, if masking
                        simulated data
  --outStatsDir OUTSTATSDIR
                        Path to a directory where values of each statistic in
                        each subwindow are recorded for each rep
  --ancFileName ANCFILENAME
                        Path to a fasta-formatted file that contains inferred
                        ancestral states ('N' if unknown). This is used for
                        masking, as sites that cannot be polarized are masked,
                        and we mimic this in the simulted data. Ignored in
                        diploid mode which currently does not use ancestral
                        state information
  --pMisPol PMISPOL     The fraction of sites that will be intentionally
                        polarized to better approximate real data
```
This mode takes three arguments and then offers many options. The arguments are the "shicMode", i.e. whether
to calculate the haploid or diploid summary statistics, the name of the input file, and the name of the output file. 
The various options allow one to account for missing data (via masking), unfolding the site frequency spectrum via the ancestral
states file (haploid only), and a mis-polarization rate of that unfolded site frequency spectrum. Please see the example usage below
for a fleshed out example of how to use these features.

#### fvecVcf mode
The fvecVcf mode is used for calculating feature vectors from data that is stored as a VCF file. 
The help message from this mode is as follows
```
$ python diploSHIC.py fvecVcf -h
usage: diploSHIC.py fvecVcf [-h] [--targetPop TARGETPOP]
                            [--sampleToPopFileName SAMPLETOPOPFILENAME]
                            [--winSize WINSIZE] [--numSubWins NUMSUBWINS]
                            [--maskFileName MASKFILENAME]
                            [--unmaskedFracCutoff UNMASKEDFRACCUTOFF]
                            [--ancFileName ANCFILENAME]
                            [--statFileName STATFILENAME]
                            [--segmentStart SEGMENTSTART]
                            [--segmentEnd SEGMENTEND]
                            shicMode chrArmVcfFile chrArm chrLen

required arguments:
  shicMode              specifies whether to use original haploid SHIC (use
                        'haploid') or diploSHIC ('diploid')
  chrArmVcfFile         VCF format file containing data for our chromosome arm
                        (other arms will be ignored)
  chrArm                Exact name of the chromosome arm for which feature
                        vectors will be calculated
  chrLen                Length of the chromosome arm
  fvecFileName          path to file where feature vectors will be written

optional arguments:
  -h, --help            show this help message and exit
  --targetPop TARGETPOP
                        Population ID of samples we wish to include
  --sampleToPopFileName SAMPLETOPOPFILENAME
                        Path to tab delimited file with population
                        assignments; format: SampleID popID
  --winSize WINSIZE     Length of the large window (default=1100000)
  --numSubWins NUMSUBWINS
                        Number of sub-windows within each large window
                        (default=11)
  --maskFileName MASKFILENAME
                        Path to a fasta-formatted file that contains masking
                        information (marked by 'N'); must have an entry with
                        title matching chrArm
  --unmaskedFracCutoff UNMASKEDFRACCUTOFF
                        Fraction of unmasked sites required to retain a
                        subwindow
  --ancFileName ANCFILENAME
                        Path to a fasta-formatted file that contains inferred
                        ancestral states ('N' if unknown); must have an entry
                        with title matching chrArm. Ignored for diploid mode
                        which currently does not use ancestral state
                        information.
  --statFileName STATFILENAME
                        Path to a file where statistics will be written for
                        each subwindow that is not filtered out
  --segmentStart SEGMENTSTART
                        Left boundary of region in which feature vectors are
                        calculated (whole arm if omitted)
  --segmentEnd SEGMENTEND
                        Right boundary of region in which feature vectors are
                        calculated (whole arm if omitted)
```
This mode takes five arguments and again has many options. The required arguments are the "shicMode", i.e. whether
to calculate the haploid or diploid summary statistics, the name of the input file, which chromosome to arm to calculate
statistics for, the length of that chromosome, and the name of the output file.

### training the CNN and prediction
Once we have feature vector files ready to go we can train and test our CNN and then finally do prediction on empirical data.

### formatting our training set
Before entering train mode we need to consolidate our training set into 5 files, one for each class. This is done using the
makeTrainingSets mode whose help message is as follows:
```
$ python diploSHIC.py makeTrainingSets -h
usage: diploSHIC.py makeTrainingSets [-h]
                                     neutTrainingFileName
                                     softTrainingFilePrefix
                                     hardTrainingFilePrefix
                                     sweepTrainingWindows
                                     linkedTrainingWindows outDir

required arguments:
  neutTrainingFileName  Path to our neutral feature vectors
  softTrainingFilePrefix
                        Prefix (including higher-level path) of files
                        containing soft training examples; files must end with
                        '_$i.$ext' where $i is the subwindow index of the
                        sweep and $ext is any extension.
  hardTrainingFilePrefix
                        Prefix (including higher-level path) of files
                        containing hard training examples; files must end with
                        '_$i.$ext' where $i is the subwindow index of the
                        sweep and $ext is any extension.
  sweepTrainingWindows  comma-separated list of windows to classify as sweeps
                        (usually just '5' but without the quotes)
  linkedTrainingWindows
                        list of windows to treat as linked to sweeps (usually
                        '0,1,2,3,4,6,7,8,9,10' but without the quotes)
  outDir                path to directory where the training sets will be
                        written

optional arguments:
  -h, --help            show this help message and exit
```
#### train mode
Here is the help message for the train mode of `diploSHIC.py`
```
$ python diploSHIC.py train -h
usage: diploSHIC.py train [-h] [--epochs EPOCHS] [--numSubWins NUMSUBWINS]
                          trainDir testDir outputModel

required arguments:
  trainDir              path to training set files
  testDir               path to test set files, can be same as trainDir
  outputModel           file name for output model, will create two files one
                        with structure one with weights

optional arguments:
  -h, --help            show this help message and exit
  --epochs EPOCHS       max epochs for training CNN (default = 100)
  --numSubWins NUMSUBWINS
                        number of subwindows that our chromosome is divided
                        into (default = 11)
```
As you will see in a moment train mode is used for training the deep learning classifier. Its required
arguments are trainDir (the directory where the training feature vectors
are kept), testDir (the directory where the testing feature vectors are kept), and outputModel the file name for the trained
network. One note -- `diploSHIC.py` expects five files named `hard.fvec`, `soft.fvec`, `neut.fvec`, `linkedSoft.fvec`, and 
`linkedHard.fvec` in the training and testing directories. The training and testing directory can be the same directory in 
which case 20% of the training examples are held out for use in testing and validation.

train mode has two options, the number of subwindows used for the feature vectors and the number of training epochs for the
network.

### predict mode
Once a classifier has been trained, one uses the predict mode of `diploSHIC.py` to classify empirical data. Here is the help
statement
```
$ python diploSHIC.py predict -h
usage: diploSHIC.py predict [-h] [--numSubWins NUMSUBWINS]
                            modelStructure modelWeights predictFile
                            predictFileOutput

required arguments:
  modelStructure        path to CNN structure .json file
  modelWeights          path to CNN weights .h5 file
  predictFile           input file to predict
  predictFileOutput     output file name

optional arguments:
  -h, --help            show this help message and exit
  --numSubWins NUMSUBWINS
                        number of subwindows that our chromosome is divided
                        into (default = 11)
```
The predict mode takes as input the two model files output by the train mode, an input file of empirical feature 
vectors, and a file name for the prediction output. 

#### a quick example of the train/predict cycle
We have supplied in the repo some example data that can give you a quick run through the train/predict cycle (we will also
shortly provide a soup-to-nuts example that starts by calculating feature vectors from simulations and ends with prediction of 
genomic data). Let's quickly give that code a spin. The directories `testing/` and `training/` each contain appropriately
formatted diploid feature vectors that are ready to be fed into diploSHIC. First we will train the diploSHIC CNN, but we will
restrict the number of training epochs to 10 to keep things relatively brief (this runs in less than 5 minutes on our server). 
```
$ python diploSHIC.py train training/ testing/ fooModel --epochs 10
```
as it runs a bunch of information monitoring the training of the network will apear. We are tracking the loss and accuracy in the
validation set. When optimization is complete our trained network will be contained in two files, `fooModel.json` and 
`fooModel.weights.hdf5`. The last bit of output from `diploSHIC.py` gives us information about the loss and accuracy on
the held out test data. From the above run my looks like this:
```
evaluation on test set:
diploSHIC loss: 0.404791
diploSHIC accuracy: 0.846800
```
Not bad. In practice I would set the `--epochs` value much higher than 10- the default setting of 100 should suffice in most cases.
Now that we have a trained model we can make predictions on some empirical data. In the repo there is a file called `testEmpirical.fvec`
that we will use as input
```
$ python diploSHIC.py predict fooModel.json fooModel.weights.hdf5 testEmpirical.fvec testEmpirical.preds
```
the output predictions will be saved in `testEmpirical.preds` and should be straightforward to interpret.

### A complete test case
In the interest of showing the user the whole enchilada when it comes to the workflow, I've provided the user
with a more detailed example on the wiki of this repo. That example can be found here: https://github.com/kern-lab/diploSHIC/wiki/A-soup-to-nuts-example