Sequence mining and genomic neighborhood analysis of mCpol

This repository contains code relevant to the publication A miniature CRISPR-Cas10 enzyme confers immunity by inhibitory signaling, by the authors Erin E. Doherty, Benjamin A. Adler, Peter H. Yoon, Kendall Hsieh, Kenneth Loi, Emily G. Armbruster, Arushi Lahiri, Cydni S. Bolling, Xander E. Wilcox, Amogha Akkati, Anthony T. Iavarone, Joe Pogliano, Jennifer A. Doudna. 

READ ME: 

This work flow uses python3, with following modules from biopython: Bio: Entrez, SeqIO and  Bio.SeqFeature: SeqFeature, FeatureLocation

The workflow also relies on mmseqs2, and was tested with version (6f45232ac8daca14e354ae320a4359056ec524c2) 

An MSA (MCPOL_CUSTOM_MSA_HMMER_INPUT.fasta) generated using the DALI/T-coffee method published previously (Yoon et al., 2024) and used as query for a HMMER search on MPI web server against NR_AUG_2019 (0 iterations hhblits, e-value e-10). 

This resulted in 512 total hits of full-length sequences. 

A custom python script was used to download the genomic loci for these hits, with a 10kb window (20kb total)

python /home/petery/scripts/ncbi_protein_id_to_genbank_window.py 11262024_hmmer_e10_filtered.txt 11262024_hmmer_e10_filtered.genbank 11262024_hmmer_e10_filtered.notfound --window 10000 && 

A custom script was also used to reorient the hits to be facing the same direction and to color the gene of interest for easy visualization in Geneious.

python /home/petery/scripts/color_protein_hits_genbank_flip_orientation.py 11262024_hmmer_e10_filtered.txt 11262024_hmmer_e10_filtered.genbank 11262024_hmmer_e10_filtered.genbank.colored
/groups/doudna/projects/yoon/cas10/mcpol/neighbor2/11262024_hmmer_e10_filtered.genbank.colored 

Resulting genbank was imported to geneious, then exported as FASTA in Geneious to be used in MMseqs.

Ran mmseqs:
mmseqs easy-cluster 11262024_hmmer_e10_filtered.genbank_len_20.fasta mcpol_20_0505 tmpDir -c 0.5 --min-seq-id 0.5 --dbtype 2 --cluster-mode 2

Took representative ID of mmseqs

cut -f 1 mcpol_20_0505_cluster.tsv > mcpol_20_0505_cluster.rep.txt && sort mcpol_20_0505_cluster.rep.txt | uniq > mcpol_20_0505_cluster.rep.unique.txt

Filtered to keep only mmseqs representatives using a custom python script: 

python /home/petery/scripts/filter_genbank_accession_v2.py 11262024_hmmer_e10_filtered.genbank.colored mcpol_20_0505_cluster.rep.unique.txt 11262024_hmmer_e10_filtered.genbank.colored.mmseqs_0505.gb

Imported to geneious. Then exported as FASTA to be used in DefenseFinder (webserver)

Defense finder outputs were then fed into this python script to annotate the existing genbank:

python /home/petery/scripts/add_annotations_defense_finder_genes.py -f defense-finder-proteins-3ce3e2281e4506fb83dab3f561a6edec-44c769zv.faa -g 11262024_hmmer_e10_filtered.genbank.colored.mmseqs_0505.gb -o 11262024_hmmer_e10_filtered.genbank.colored.mmseqs_0505.defense.annotated.gb -hmmer defense-finder-genes-3ce3e2281e4506fb69a9511a1b24d8ef-28jndi6v.tsv