dgmfinder 0.0.0
Package to detect diversity generating mechanisms from sequencing data alone.
License: Apache License 2.0
Keywords: python,crispr,dgr

# dgmfinder

## Installation

To install dgmfinder simply run the following commands:

```python
git clone https://github.com/jordiabante/dgmfinder.git
cd dgmfinder
conda create -n dgmfinder python=3.6.8
conda activate dgmfinder
pip install -e .
```

## Usage

To run a single-sample analysis

```python
# import bio module
from dgmfinder import bio

# define input file and output dir
outdir = "/path/to/outdir/"
fqfile = "/path/to/fastq.gz"

# configure run
config = bio.Config(
    outdir=outdir,              # directory where output files will be stored
    kmer_size=27,               # k-mer size used in the analysis
    min_smp_sz=5,               # minimum sample size to compute p-value
    max_smp_sz=50,              # maximum number of sequences sampled per
    lmer_size=7,                # l-mer size used to compute jaccard similarity between k-mers
    jsthrsh=0.25                # jaccard similarity threshold used to collapse the observed sequences
    max_fastq_reads=5000000,    # maximum number of FASTQ reads to process
    annot_fasta=[]              # array containing fasta files to use with blast
)

# run analysis
bio.dgmfinder_single_sample_analysis(fqfile,config)
```

If `len(annot_fasta)>0`, then `blast` must be available on the path.

## Output files

### Anchor sequences table

A table with suffix `_anchors.txt.gz` is produced containing information about the interesting anchors detected (keys in old convention). The file contains the following columns:

```bash
    sequence id | assembly of {anchor1,anchor2,...} | max_c_up | max_n_up | max_efct_sz_up | max_efct_sz_qval_up | max_kmer_up | max_c_dn | max_n_dn | max_efct_sz_dn | max_efct_sz_qval_dn | max_kmer_dn | A% | C% | G% | T% | {anchor1,anchor2,...}
```

where `up/dn` indicate the position of the HVR with respect to the anchor and:

* `max_c_*`: number of clusters formed for the maximizing anchor in the set in `*` direction.
* `max_n_*`: corresponding number of target sequences observed.
* `max_efct_sz_*`: corresponding effect size.
* `max_efct_sz_qval_*`: corresponding adjusted p-value.
* `max_kmer_*`: corresponding k-mer sequence.

If the `len(annot_fasta)>0`, then two extra columns will be added to the previous table, for each direction (upstream, downstream) and for each FASTA in `annot_fasta`, containing the lowest e-value and the corresponding hit in the FASTA (sequence in FASTA resulting in lowest e-value), and the output will be stored in a new table with suffix `_anchors_annot.txt.gz` (NA will be assigned when e-value>1). For example, if we pass `annot_fasta=[fasta1]` we will see four extra columns:

```bash
    sequence id | ... | {anchor1,anchor2,...} | best_eval_up_fasta1 | best_hit_up_fasta1 | best_eval_dn_fasta1 | best_hit_dn_fasta1 
```

The intermediate XML files produced by blast are also stored for further analysis.

### Re-running annotation

In some cases we might want to update the set of FASTA files we want to blast the results against. Say, for example, that we want to re-run the annotation with FASTA files `f1.fasta`, `f2.fasta`, and `f3.fasta`, with our output `SRRXYZ_anchors.txt.gz` (note the suffix is `_anchors.txt.gz`). In that case, we can use the following python code:

```python
from dgmfinder import bio

anchorfile = "/path/to/SRRXYZ_anchors.txt.gz"
annot_fasta = ["/path/to/annotations/f1.fa", "/path/to/annotations/f2.fa", "/path/to/annotations/f3.fa"]
config = bio.Config(annot_fasta=annot_fasta)

bio.dgmfinder_single_sample_analysis_annotation(anchorfile,config)
```

### Anchor sequences FASTA

Three FASTA files are produced:

1. FASTA file with suffix `_assemb_anchors.fasta`: assembled anchor sequences.
2. FASTA file with suffix `_max_anchor_up.fasta`: maximizing anchor sequence upstream.
3. FASTA file with suffix `_max_anchor_dn.fasta`: maximizing anchor sequence downstream.

Note that not all anchor sequences in 2 and 3 are necessarily significant.

### Target sequences table

For each anchor in the set `{anchor1,anchor2,...}`, the target sequences are stored in a file with suffix `_targets.txt.gz` containing the following columns:

```bash
    anchor | upstream/downstream | distance | target | number of instances observed
```

## Command line

A set of useful command line tools can be found in the `bash/bin/` folder. To use the commands, add the path to this folder to your environment `PATH`, and to check usage simply type their name on the terminal. Currently, we have the following commands:

* `dgmfinder_header`: adds header to `_anchors_annot.txt.gz` file. E.g., if `annot_fasta=[a.fasta,b.fasta,c.fasta]`, then

```bash
    dgmfinder_header test_data_anchors_annot.txt.gz a,b,c
```

* `dgmfinder_get_context`: prints out the sequences sorrounding a given kmer. E.g., to look for the sorrounding sequences around k-mer ATGATGA in FASTQ file `sample_prefix.fastq.gz` we can run

```bash
    dgmfinder_get_context sample_prefix.fastq.gz ATGATGA
```

* `fastq_trim`: wrapper around Trimgalore to trim FASTQ files. E.g., if we have a sample `sample_prefix.fastq.gz` we can run

```bash
    fastq_trim sample_prefix
```