If you use the a2z model in your work, please cite the following article:
Wrightsman, T., Marand, A. P., Crisp, P. A., Springer, N. M., & Buckler, E. S. (2022). Modeling chromatin state from sequence across angiosperms using recurrent convolutional neural networks. The Plant Genome, 15, e20249. https://doi.org/10.1002/tpg2.20249
The authors recommend mambaforge as a faster alternative to conda.
If using mambaforge, use the mamba create command instead of conda create.
git clone https://github.com/twrightsman/a2z-regulatory
cd a2z-regulatory
conda env create -f envs/a2z.yml
conda activate a2z
pip install src/python/a2z
After setting up the environment, you can then use the scripts in data/ to download all the necessary data for the experiments in the results/ folder.
This is only useful for better git-diff output on Jupyter Notebooks.
git config --local --add diff.jupyternotebook.command "git-nbdiffdriver diff"
git config --local --add merge.jupyternotebook.command "git-nbmergedriver merge %O %A %B %L %P"
echo "*.ipynb diff=jupyternotebook merge=jupyternotebook" >> .git/info/attributes
Use predict_genome.py -h to see all available options.
(base)$ curl 'https://zenodo.org/record/5724562/files/model-accessibility-full.h5?download=1' > model-accessibility-full.h5
(base)$ conda activate a2z
(a2z)$ src/python/scripts/predict_genome.py --stride 50 model-accessibility-full.h5 genome.fa > preds.bed
Now you will have a BED3+1 file with predictions in 600bp windows sliding 50bp at a time.
You can do simple "peak-calling" with bedtools.
(a2z)$ samtools faidx genome.fa
(a2z)$ bedtools slop -g genome.fa.fai -b -150 < preds.bed > preds.trimmed.bed
(a2z)$ bedtools map -a <(bedtools makewindows -g genome.fa.fai -w 50) -b preds.trimmed.bed -c 4 -o max | awk '$4 >= 0.9' | bedtools merge > preds.gte90.merged.bed
The above will trim 150bp off each end of the predictions and take the max prediction of all overlapping windows every 50bp and then merge touching 50bp chunks with predictions greater than or equal to 0.90 into whole regions. Trimming is optional but recommended because the size of true accessible regions called from the ATAC-seq data is close to 300bp.
You can also mask out repeats with bedtools and an annotation.
(a2z)$ grep -v '^#' repeats.gff3 | awk -v 'OFS=\t' '{print $1, $4 - 1, $5}' > repeats.bed
(a2z)$ bedtools subtract -A -f 0.5 -a preds.gte90.merged.bed -b repeats.bed > preds.gte90.merged.masked.bed
The a2z leaf chromatin and DNA methylation models are also available in the Kipoi model zoo.
The following will make leaf chromatin accessibility predictions on 600bp sliding windows with a stride of 50bp.
You can adjust the window stride by modifying the -s parameter to bedtools.
(base)$ ls
a2z-regulatory
(base)$ mkdir work
(base)$ cd work
(base)$ curl https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1/GCF_000001735.4_TAIR10.1_genomic.fna.gz | gzip -cd > tair.fa
(base)$ conda activate a2z
(a2z)$ samtools faidx tair.fa
(a2z)$ bedtools makewindows -g tair.fa.fai -w 600 -s 50 | awk '$3-$2 == 600' > windows.bed
(a2z)$ conda env create -f ../a2z-regulatory/envs/kipoi.yml
(a2z)$ conda activate a2z-kipoi
(a2z-kipoi)$ kipoi env export a2z-chromatin/a2z-accessibility -o env.yml
(a2z-kipoi)$ conda env create -f env.yml
(a2z-kipoi)$ conda activate kipoi-a2z-chromatin__a2z-accessibility
(kipoi-a2z-chromatin__a2z-accessibility)$ kipoi predict a2z-chromatin/a2z-accessibility -n 8 --batch_size=128 --dataloader_args='{"intervals_file": "windows.bed", "fasta_file": "tair.fa", "num_chr_fasta": true}' -o a2z.acc.preds.tsv