---
# System prepended metadata

title: Craig and the origin of the marine picocyanos*
tags: [Craig cyanos]

---

---
tags: Craig cyanos
---

# Craig and the origin of the marine picocyanos*
> \* Available wherever books are sold.


[toc]

## Programs/versions and conda environment setup

The primary conda environment to do the following work can be setup and activated as follows:

```bash
conda create -n cyanos -c conda-forge -c bioconda -c defaults -c astrobiomike gtotree=1.5.36 bit=1.8.11 ppanggolin=1.1.96 scoary=1.6.16 fastani=1.31

conda activate cyanos
```

Anvi'o environment installed separately, noted to be activated when used:

```bash
conda create -n anvio-6.2 -c conda-forge -c bioconda -c defaults anvio=6.2
```

# Working with all Cyanos

## Getting NCBI RefSeq Cyanobacteria accessions

Searching for all NCBI RefSeq "[complete](https://www.ncbi.nlm.nih.gov/assembly/help/#glossary)" and "[representative](https://www.ncbi.nlm.nih.gov/refseq/about/prokaryotes/#representative_genomes)" genomes retrieved 174 accessions as searched on 10-Sept-2020:

```bash
esearch -query 'Cyanobacteria[ORGN] AND "latest refseq"[filter] AND "complete genome"[filter] AND (latest[filter] AND all[filter] NOT anomalous[filter])' \
        -db assembly | esummary | xtract -pattern DocumentSummary -def "NA" -element AssemblyAccession > Cyano-refseq-complete-accs-10-Sept-2020.txt
        
wc -l Cyano-refseq-complete-accs-10-Sept-2020.txt
    # 174
```

These are what we will be working with to start.

Getting their taxids in a table too to be able to get full lineage info for all:

```bash
esearch -query 'Cyanobacteria[ORGN] AND "latest refseq"[filter] AND "complete genome"[filter] AND (latest[filter] AND all[filter] NOT anomalous[filter])' \
        -db assembly | esummary | xtract -pattern DocumentSummary -def "NA" -element AssemblyAccession,Taxid > Cyano-refseq-complete-accs-and-taxids-10-Sept-2020.tsv
```

Getting lineage info with my toolkit:

```bash
cut -f 2 Cyano-refseq-complete-accs-and-taxids-10-Sept-2020.tsv > target-taxids.tmp

bit-get-lineage-from-taxids -i target-taxids.tmp -o lineages.tsv

  # making sure still in same order before sticking together
diff <( cut -f 2 Cyano-refseq-complete-accs-and-taxids-10-Sept-2020.tsv ) <( cut -f 1 lineages.tsv | tail -n +2 )

paste <( cat <( printf "accession\n" ) <( cut -f 1 Cyano-refseq-complete-accs-and-taxids-10-Sept-2020.tsv ) ) <( cut -f 2- lineages.tsv ) | sed 's/specific_name/species/'> Cyano-refseq-complete-accs-and-lineages-10-Sept-2020.tsv

rm lineages.tsv target-taxids.tmp
```

And here's our table with accessions and full lineages:

```bash
column -ts $'\t' Cyano-refseq-complete-accs-and-lineages-10-Sept-2020.tsv | head | sed 's/^/# /'
# accession        domain    phylum         class          order                  family                  genus                           species
# GCF_014489415.1  Bacteria  Cyanobacteria  NA             Nostocales             Aphanizomenonaceae      Cylindrospermopsis              Cylindrospermopsis curvispora
# GCF_014304815.1  Bacteria  Cyanobacteria  NA             Synechococcales        Synechococcaceae        Synechococcus                   Synechococcus sp. MIT S9220
# GCF_014304795.1  Bacteria  Cyanobacteria  NA             Synechococcales        Synechococcaceae        Synechococcus                   Synechococcus sp. WH 8101
# GCF_014280235.1  Bacteria  Cyanobacteria  NA             Synechococcales        Synechococcaceae        Cyanobium                       Cyanobium sp. NS01
# GCF_014280215.1  Bacteria  Cyanobacteria  NA             Synechococcales        Synechococcaceae        Synechococcus                   Synechococcus sp. A15-127
# GCF_014280195.1  Bacteria  Cyanobacteria  NA             Synechococcales        Synechococcaceae        Synechococcus                   Synechococcus sp. A15-24
# GCF_014280175.1  Bacteria  Cyanobacteria  NA             Synechococcales        Synechococcaceae        Synechococcus                   Synechococcus sp. A15-28
# GCF_014280115.1  Bacteria  Cyanobacteria  NA             Synechococcales        Synechococcaceae        Synechococcus                   Synechococcus sp. A15-44
# GCF_014280095.1  Bacteria  Cyanobacteria  NA             Synechococcales        Synechococcaceae        Synechococcus                   Synechococcus sp. A15-60
```

## Phylogenomic trees
These are using the Cyanobacteria-specific HMM packaged with GToTree that targets 251 single-copy genes.

### Standard [FastTree](http://www.microbesonline.org/fasttree/)
```bash
GToTree -a Cyano-refseq-complete-accs-10-Sept-2020.txt -H Cyanobacteria \
        -t -L Species,Strain -j 15 -o Cyano-refseq-gtotree-out
```

**Circular interactive tree can be explored [here](https://itol.embl.de/tree/7184244121413481600300947).**

**Normal interactive tree can be explored [here](https://itol.embl.de/tree/7184244121479331600301304).**


<a href="https://i.imgur.com/nIjc89S.png"><img src="https://i.imgur.com/nIjc89S.png"></a>


Used interactive tree to get the IDs of those that are in long-branch marine clade and put them in a file called `long-branch-marine-buggers.tmp`.

Formatting to get them to be regular accessions again:

```bash
cut -f 1,2 -d " " long-branch-marine-buggers.tmp | tr " " "_" > long-branch-marine-accs.txt && rm long-branch-marine-buggers.tmp
```

This is 57 of the total 174:

```bash
wc -l long-branch-marine-accs.txt
  # 57
```

### IQTREE2 with mixed-model
This will take a while to finish, but when we get to talking about individual gene trees as compared to the phylogenomic tree, it'll be good to have a more robustly built tree.

```bash
iqtree -s Cyano-refseq-gtotree-out/Aligned_SCGs_mod_names.faa \
       -spp Cyano-refseq-gtotree-out/run_files/Partitions.txt \
       -m MFP -bb 1000 -nt 20 -pre Cyano-gtotree-refseq-iqtree-mixed-model-out
```

*RUNNING* on kuat server in `gtotree-cyano` screen

## Pangenomics

Doing within [anvio](http://merenlab.org/software/anvio/) which uses [mcl](https://micans.org/mcl/) for clustering based on amino-acid alignment bitscores.


```bash
conda activate anvio-6.2
```

Main setup and processing is in the following script:

```bash
cat run-anvio-pan.sh
```
```bash
printf "name\tcontigs_db_path\n" > external-genomes.tsv

for acc in $(cat Cyano-refseq-complete-accs-10-Sept-2020.txt)
do
    printf "\n  On: ${acc}\n"
    gunzip -c ${acc}.gb.gz > ${acc}.gb
    
    anvi-script-process-genbank -i ${acc}.gb --output-fasta ${acc}.fa --output-gene-calls ${acc}-gene-calls.tsv \
                                --output-functions ${acc}-functions.tsv

    rm ${acc}.gb
    
    anvi-gen-contigs-database -f ${acc}.fa -o ${acc}-contigs.db -L 0 --external-gene-calls ${acc}-gene-calls.tsv 

    anvi-import-functions -c ${acc}-contigs.db -i ${acc}-functions.tsv

    # genome names can't have periods in anvio, so changing them to underscores in here too
    new_acc=$(echo $acc | sed 's/\./_/')

    printf "$new_acc\t${acc}-contigs.db\n" >> external-genomes.tsv

done

# generating genomes storage file
anvi-gen-genomes-storage -e external-genomes.tsv --gene-caller NCBI_PGAP \
                         -o Cyano-RefSeq-174-GENOMES.db 

# running pangenome
anvi-pan-genome -g Cyano-RefSeq-174-GENOMES.db --exclude-partial-gene-calls -o Cyano-RefSeq-174-anvio-pan -T 20 -n cyanos
```

```bash
bash run-anvio-pan.sh
```

Generating summary table:

```bash
anvi-script-add-default-collection -p Cyano-RefSeq-174-anvio-pan/cyanos-PAN.db -C all

anvi-summarize -p Cyano-RefSeq-174-anvio-pan/cyanos-PAN.db -g Cyano-RefSeq-174-GENOMES.db -C all -o Cyano-RefSeq-174-pan-summary
```

Getting those GCs that show up in at least 2 genomes so can filter down the pangenome and possibly visualize (it's too big right now):

```bash
gunzip Cyano-RefSeq-174-pan-summary/cyanos_gene_clusters_summary.txt.gz

awk -F $'\t' ' { if ( $6 > 1 ) print $2 } ' Cyano-RefSeq-174-pan-summary/cyanos_gene_clusters_summary.txt | tail -n +2 | sort -u > min-2-GCs.txt

wc -l min-2-GCs.txt
# 43579

## using the "X" to be portable as it will work on darwin shell's too
sed 's/$/Xno_singletons/' min-2-GCs.txt | tr "X" "\t" > no-singletons-collection.tsv

anvi-import-collection -p Cyano-RefSeq-174-anvio-pan/cyanos-PAN.db -C no_singletons no-singletons-collection.tsv
```

Splitting pangenome to make on with no singletons:

```bash
anvi-split -p Cyano-RefSeq-174-anvio-pan/cyanos-PAN.db -g Cyano-RefSeq-174-GENOMES.db -C no_singletons -o Cyano-RefSeq-174-anvio-pan-no-singletons-split
```


No good there either, there are a lot of gene clusters with this many genomes. Since we mostly care about the clear differences between the two groups (i.e., no biggie if we see 5 pico marines with a gene that no others have – given our current scope), going to limit it to just those gene clusters with at least 10 genomes contributing:

```bash
awk -F $'\t' ' { if ( $6 > 9 ) print $2 } ' Cyano-RefSeq-174-pan-summary/cyanos_gene_clusters_summary.txt | tail -n +2 | sort -u > min-10-GCs.txt

wc -l min-10-GCs.txt
# 9695

## using the "X" to be portable as it will work on darwin shell's too
sed 's/$/Xmin_ten/' min-10-GCs.txt | tr "X" "\t" > min-ten-collection.tsv

anvi-import-collection -p Cyano-RefSeq-174-anvio-pan/cyanos-PAN.db -C min_ten min-ten-collection.tsv
```

Splitting pangenome to make one with minimum 10:

```bash
anvi-split -p Cyano-RefSeq-174-anvio-pan/cyanos-PAN.db -g Cyano-RefSeq-174-GENOMES.db -C min_ten -o Cyano-RefSeq-174-anvio-pan-min-ten-split
```

Visualizing:

```bash
  # ssh'd in like: ssh -L 8080:localhost:8080 <user>@<host>
anvi-display-pan -p Cyano-RefSeq-174-anvio-pan-min-ten-split/min_ten/PAN.db -g Cyano-RefSeq-174-GENOMES.db --server-only -P 8080
```

Can't share the interactive, but here's a snapshot of the overview:

<a href="https://i.imgur.com/APU3ZpV.png"><img src="https://i.imgur.com/APU3ZpV.png"></a>


Also running a new pangenome with min 10:

```bash
anvi-pan-genome -g Cyano-RefSeq-174-GENOMES.db --exclude-partial-gene-calls -o Cyano-RefSeq-174-anvio-pan-min-10 -T 30 -n cyanos --min-occurrence 10
```


### Looking for gene-cluster enrichment/depletion in marine pico clade genomes

Anvi'o has a [built-in program](http://merenlab.org/2016/11/08/pangenomics-v2/#making-sense-of-functions-in-your-pangenome) to search for enriched functions. 

```bash
conda activate anvio-6.2
```

Generating and importing required table as [needed](http://merenlab.org/2017/12/11/additional-data-tables/#layers-additional-data-table):

```bash
# doing with the "X" first so portable/works on darwin shell's too, and changing . to _ to match how they names needed to be in anvio
sed 's/$/Xmarine/' long-branch-marine-accs.txt | tr "X" "\t" | tr "." "_" > p1.tmp

# getting rest
comm -23 <( sort -V Cyano-refseq-complete-accs-10-Sept-2020.txt ) <( sort -V long-branch-marine-accs.txt ) | sed 's/$/Xother/' | tr "X" "\t" | tr "." "_" > p2.tmp

# combining and adding header
cat <( printf "samples\tclade\n" ) p1.tmp p2.tmp > trait-tab-for-anvio.tsv

rm *.tmp
```

Importing:

```bash
anvi-import-misc-data -p Cyano-RefSeq-174-anvio-pan/cyanos-PAN.db --target-data-table layers trait-tab-for-anvio.tsv
```

Running functional enrichment test based on clade:

```bash
anvi-get-enriched-functions-per-pan-group -p Cyano-RefSeq-174-anvio-pan/cyanos-PAN.db -g Cyano-RefSeq-174-GENOMES.db --category clade --annotation-source NCBI_PGAP -o Clade-enriched-functions.txt --functional-occurrence-table-output Function-occurrence-frequencies.tsv
```

#### Enriched/depleted functions in marine pico clade genomes

This was pitting the blue-colored clade highlighted in this snapshot from above against the rest (green and grey):

<a href="https://i.imgur.com/nIjc89S.png"><img src="https://i.imgur.com/nIjc89S.png"></a>

There are 57 within the long branch (blue), and 117 outside of it. 

The full table of enriched functions can be downloaded by clicking [here](https://ndownloader.figshare.com/files/24711836). 

> "Prop. in marine clade" is the proportion of the 57 pico marine cyano genomes that had that function annotated within them, and "Prop. other" is the proportion of the other 117 genomes that had that function annotated within them. 

**Some enriched in marine clade**
| Function | Prop. marine clade | Prop. other |
|:---------|------:|------:|
| carboxysome peptide B | 	1 | 	0 |
| potassium transporter TrkG | 	1 | 	0 |
| ATP adenylyltransferase | 	1 | 	0 |
| non-canonical purine NTP pyrophosphatase | 	1 | 	0 |
| glutamate racemase (murI) | 	1 | 	0 |
| NADPH-dependent assimilatory sulfite reductase hemoprotein subunit | 	1 | 	0 |
| carboxysome shell carbonic anhydrase | 	1 | 	0 |
| anthranilate synthase component I family protein | 	1 | 	0 |
| DUF3136 domain-containing protein | 	1 | 	0 |
| glutathione-disulfide reductase (gorA) | 	1 | 	0 |
| DUF1957 domain-containing protein | 	1 | 	0 |
| NAD\(P\)H-quinone oxidoreductase subunit O | 	1 | 	0 |
| NAD\(P\)H-quinone oxidoreductase subunit L | 	1 | 	0 |
| histidine phosphotransferase | 	1 | 	0 |
| DUF3188 domain-containing protein | 	1 | 	0 |
| carboxysome peptide A | 	1 | 	0 |
| carboxysome shell protein | 	1 | 	0 |
| AIR synthase | 	1 | 	0 |
| precorrin-2 C(20)-methyltransferase (cobI) | 	1 | 	0 |
| DUF296 domain-containing protein | 	0.9825 | 	0 |
| arylesterase | 	0.9825 | 	0 |
| DUF3721 domain-containing protein | 	0.9825 | 	0 |
| DUF1651 domain-containing protein | 	0.9825 | 	0 |
| SprT family zinc-dependent metalloprotease | 	0.9825 | 	0 |
| dolichol kinase | 	0.9825 | 	0 |
| uracil phosphoribosyltransferase | 	0.9474 | 	0 |
| malate:quinone oxidoreductase | 	0.9474 | 	0 |
| photosystem II reaction center PsbP family protein | 	1 | 	0.0342 |
| glucose-6-phosphate dehydrogenase assembly protein OpcA | 	1 | 	0.0342 |
| DUF814 domain-containing protein | 	0.9298 | 	0 |
| phosphonate ABC transporter | 	0.9298 | 	0 |
| thiamine-phosphate kinase (thiL) | 	0.9649 | 	0.0171 |
| LptF/LptG family permease | 	1 | 	0.0427 |
| thylakoid membrane photosystem I accumulation factor | 	1 | 	0.0427 |
| ClC family H(+)/Cl(-) exchange transporter | 	0.9123 | 	0 |
| EF-1 guanine nucleotide exchange domain-containing protein | 	0.9123 | 	0 |
| homoserine O-succinyltransferase | 	0.9123 | 	0 |
| DUF3764 family protein | 	0.8772 | 	0 |
| Fe-S cluster containing protein | 	0.8772 | 	0 |
| COX15/CtaA family protein | 	0.8596 | 	0 |
| Nif11 domain/cupin domain-containing protein | 	0.8596 | 	0 |
| PLP-dependent transferase | 	0.8596 | 	0 |
| N-acetylglucosamine-6-phosphate deacetylase | 	0.8596 | 	0 |
| S-isoprenylcysteine methyltransferase | 	0.8596 | 	0 |
| Tic20 | 	0.8596 | 	0 |


**Some depleted in marine clade**
| Function | Prop. marine clade | Prop. other |
|:---------|------:|------:|
| Tab2/Atab2 family RNA-binding protein | 	0 | 	1 |
| ArsA family ATPase | 	0 | 	1 |
| TIGR04376 family protein | 	0 | 	1 |
| VacB/RNase II family 3'-5' exoribonuclease | 	0 | 	1 |
| penicillin-binding protein 1A | 	0 | 	1 |
| lipoprotein signal peptidase | 	0 | 	1 |
| heat-inducible transcriptional repressor HrcA (hrcA) | 	0 | 	1 |
| glutamate racemase | 	0 | 	1 |
| cytochrome C biogenesis protein CcdA | 	0 | 	1 |
| KH domain-containing protein | 	0 | 	1 |
| DNA double-strand break repair nuclease NurA | 	0 | 	1 |
| NAD\(P\)H-quinone oxidoreductase subunit N (ndhN) | 	0 | 	1 |
| DUF4330 domain-containing protein | 	0 | 	1 |
| NAD\(P\)H-quinone oxidoreductase subunit O (ndhO) | 	0 | 	1 |
| MlaE family lipid ABC transporter permease subunit | 	0 | 	1 |
| cytochrome c biogenesis protein | 	0 | 	1 |
| rRNA pseudouridine synthase | 	0 | 	0.9915 |
| S-layer homology domain-containing protein | 	0 | 	0.9915 |
| iron-containing alcohol dehydrogenase family protein | 	0 | 	0.9915 |
| ATP-dependent helicase | 	0 | 	0.9915 |
| M48 family metallopeptidase | 	0 | 	0.9915 |
| nitrogenase | 	0 | 	0.9915 |
| class 1 fructose-bisphosphatase (fbp) | 	0 | 	0.9915 |
| DNA mismatch repair endonuclease MutL (mutL) | 	0 | 	0.9915 |
| homocysteine biosynthesis protein | 	0 | 	0.9915 |
| UDP-glucose/GDP-mannose dehydrogenase family protein | 	0 | 	0.9829 |
| PP2C family protein-serine/threonine phosphatase | 	0 | 	0.9829 |
| threonylcarbamoyl-AMP synthase | 	0 | 	0.9829 |
| TRC40/GET3/ArsA family transport-energizing ATPase | 	0 | 	0.9829 |
| energy-coupling factor ABC transporter ATP-binding protein | 	0 | 	0.9829 |


The full table of enriched functions can be downloaded by clicking [here](https://ndownloader.figshare.com/files/24711836). 

---

## Looking for light/salinity-relevant genes
**COMING SOON**

* Can get targets from those not significant from enrichment analyis
    * If in all, align and do ancestral reconstruction, and do codon alignments and dN/dS

---

## Things to think about what what to try to get at them
* What makes the branch so long?
    * Separate from gene presence/absence differences, as they aren't relevant to the long branch business, there must be a good amount of selection in shared genes in addition to neutral divergence over time. 
    * Get dN/dS of every gene possible (could start with just the 251 in the phylogenomic tree), look at the distribution across all genes – which genes have a dN/dS above the Nth percentile? This *should* lead us to relevant targets i think
    * Related, does the distribution of dN/dS for all genes of out-group genomes look different than the distribution of marine pico-cyanos'?

---

```bash
conda create -n paml -c conda-forge -c bioconda -c defaults paml=4.9 pal2nal=14.1 muscle=3.8.1551 clustalo=1.2.4

conda activate paml
```


"Relative rate ratio test"...


---

# Focusing on marine pico clade and their closest relatives only

Making list of accessions from things done above and looking at the tree for the 4 outgroup ones we want:

```bash
cat long-branch-marine-accs.txt <( printf "GCF_003957805.1\nGCF_000012525.1\nGCF_000817325.1\nGCF_000010065.1\n" ) > target-picos-subset.txt
```

## Phylogenomic trees

These 4 tree variations all look virtually identical with minor differences in branch length that seem to be basically equivalent in all places, and little change in support values. Here's a snapshot of one, they all have links to their interactive trees below.

<a href="https://i.imgur.com/Jjic8aD.png"><img src="https://i.imgur.com/Jjic8aD.png"></a>

### Standard FastTree with GToTree cyanobacterial 251 SCG-set HMM
```bash
GToTree -a target-picos-subset.txt -H Cyanobacteria -t -L Species,Strain \
        -j 15 -o target-picos-gtotree-standard
```

**[Interactive tree](https://itol.embl.de/tree/7184244121410631600459049)**


### IQTREE2 with mixed-model with GToTree cyanobacterial 251 SCG-set HMM
```bash
iqtree -s target-picos-gtotree-standard/Aligned_SCGs_mod_names.faa \
       -spp target-picos-gtotree-standard/run_files/Partitions.txt \
       -m MFP -bb 1000 -nt 20 -pre Picocyano-gtotree-refseq-iqtree-mixed-model-out
```

**[Interactive tree](https://itol.embl.de/tree/7184244121411321600459068)**

### Making new SCG-set specific for these (probably overkill)
Probably overkill for the resolution we need to delineate these, but i wrote the program recently so want to use it:

```bash
gtt-gen-SCG-HMMs -a target-picos-subset.txt -n 20 -o picocyanos
```

```
  New SCG-HMMs holding 370 target genes written to:
    picocyanos/picocyanos.hmm
```

> Not terribly too many more considering the breadth of diversity we cut out.

Running standard FastTree one with those:

```bash
GToTree -a target-picos-subset.txt -H picocyanos.hmm -t -L Species,Strain \
        -j 20 -o target-picos-gtotree-picocyano-HMM
```

**[Interactive tree](https://itol.embl.de/tree/7184244121410911600459057)**

Running and IQTREE2 mixed-model one too, because why the hell not?

```bash
iqtree -s target-picos-gtotree-picocyano-HMM/Aligned_SCGs_mod_names.faa \
       -spp target-picos-gtotree-picocyano-HMM/run_files/Partitions.txt \
       -m MFP -bb 1000 -nt 20 -pre Picocyano-gtotree-refseq-picocyano-specific-HMM-iqtree-mixed-model-out
```

**[Interactive tree](https://itol.embl.de/tree/7184244121411531600459072)**


### Making new one with latest GToTree 
```
mamba create -n gtotree -c conda-forge -c bioconda -c defaults -c astrobiomike gtotree=1.5.42

conda activate gtotree
```

And keeping individual alignment files for Sarah/Amy's phylo work with new `-k` flag:

```
GToTree -a target-picos-subset.txt -H picocyanos.hmm -t -L Species,Strain -k\
        -j 20 -o target-picos-gtotree-picocyano-HMM-1.5.42
```        

> **NOTE** 
> I needed to change GCF_000063505.1 and GCF_000063525.1 to their GCA counterparts, as they were suppressed from RefSeq due to not having strain identifiers. But I wanted to keep them in here because they are WH7803 and RCC307.

## Pangenomics
Doing within [anvio](http://merenlab.org/software/anvio/) which uses [mcl](https://micans.org/mcl/) for clustering based on amino-acid alignment bitscores.


```bash
conda activate anvio-6.2
```

Main setup and processing is in the following script:

```bash
cat run-anvio-pan.sh
```
```bash
printf "name\tcontigs_db_path\n" > external-genomes.tsv

for acc in $(cat target-picos-subset.txt)
do
    printf "\n  On: ${acc}\n"
    gunzip -c ../genbank-files/${acc}.gb.gz > ${acc}.gb

    anvi-script-process-genbank -i ../genbank-files/${acc}.gb --output-fasta ${acc}.fa --output-gene-calls ${acc}-gene-calls.tsv \
                                --output-functions ${acc}-functions.tsv

    rm ${acc}.gb

    anvi-gen-contigs-database -f ${acc}.fa -o ${acc}-contigs.db -L 0 --external-gene-calls ${acc}-gene-calls.tsv

    anvi-import-functions -c ${acc}-contigs.db -i ${acc}-functions.tsv

    # genome names can't have periods in anvio, so changing them to underscores in here too
    new_acc=$(echo $acc | sed 's/\./_/')

    printf "$new_acc\t${acc}-contigs.db\n" >> external-genomes.tsv

done

# generating genomes storage file
anvi-gen-genomes-storage -e external-genomes.tsv --gene-caller NCBI_PGAP \
                         -o Picocyano-RefSeq-61-GENOMES.db

# running pangenome
anvi-pan-genome -g Picocyano-RefSeq-61-GENOMES.db --exclude-partial-gene-calls -o Picocyano-RefSeq-61-pan -T 5 -n picocyanos

# generating summary table
anvi-script-add-default-collection -p Picocyano-RefSeq-61-pan/picocyanos-PAN.db -C all
anvi-summarize -p Picocyano-RefSeq-61-pan/picocyanos-PAN.db -g Picocyano-RefSeq-61-GENOMES.db -C all -o Picocyano-RefSeq-61-pan-summary

# running pangenome with no singletons
anvi-pan-genome -g Picocyano-RefSeq-61-GENOMES.db --exclude-partial-gene-calls -o Picocyano-RefSeq-61-pan-no-singletons -T 5 -n picocyanos --min-occurrence 2

# generating summary table
anvi-script-add-default-collection -p Picocyano-RefSeq-61-pan-no-singletons/picocyanos-PAN.db -C all
anvi-summarize -p Picocyano-RefSeq-61-pan-no-singletons/picocyanos-PAN.db -g Picocyano-RefSeq-61-GENOMES.db -C all -o Picocyano-RefSeq-61-no-singletons-pan-summary
```

```bash
bash run-anvio-pan.sh
```

Making a split of only those GCs with greater than 3 genomes contributing (to thin it out more visually but still keep those unique to the 4 non-marine buggers):

```bash
gunzip Picocyano-RefSeq-61-no-singletons-pan-summary/picocyanos_gene_clusters_summary.txt.gz

awk -F $'\t' ' { if ( $6 > 3 ) print $2 } ' Picocyano-RefSeq-61-no-singletons-pan-summary/picocyanos_gene_clusters_summary.txt | tail -n +2 | sort -u > min-4-GCs.txt

wc -l min-4-GCs.txt
# 6,925

## using the "X" to be portable as it will work on darwin shell's too
sed 's/$/Xmin_4/' min-4-GCs.txt | tr "X" "\t" > min-4-collection.tsv

anvi-import-collection -p Picocyano-RefSeq-61-pan-no-singletons/picocyanos-PAN.db -C min_4_GCs min-4-collection.tsv

anvi-split -p Picocyano-RefSeq-61-pan-no-singletons/picocyanos-PAN.db -g Picocyano-RefSeq-61-GENOMES.db -C min_4_GCs -o Picocyano-RefSeq-61-pan-min-4-split
```

Peeking at this one:

```bash
  # ssh'd in like: ssh -L 8080:localhost:8080 <user>@<host>
anvi-display-pan -p Picocyano-RefSeq-61-pan-min-4-split/min_4/PAN.db -g Picocyano-RefSeq-61-GENOMES.db --server-only -P 8080
```

Can't share the interactive, but here's a screenshot:

<a href="https://i.imgur.com/scgVx2d.png"><img src="https://i.imgur.com/scgVx2d.png"></a>


### Functional enrichment testing

Anvi'o has a [built-in program](http://merenlab.org/2016/11/08/pangenomics-v2/#making-sense-of-functions-in-your-pangenome) to search for enriched functions. 

```bash
conda activate anvio-6.2
```

Generating and importing required table as [needed](http://merenlab.org/2017/12/11/additional-data-tables/#layers-additional-data-table):

```bash
# doing with the "X" first so portable/works on darwin shells too, and changing . to _ to match how they names needed to be in anvio
sed 's/$/Xmarine/' long-branch-marine-accs.txt | tr "X" "\t" | tr "." "_" > p1.tmp

# getting rest
comm -23 <( sort -V target-picos-subset.txt ) <( sort -V long-branch-marine-accs.txt ) | sed 's/$/Xother/' | tr "X" "\t" | tr "." "_" > p2.tmp

# combining and adding header
cat <( printf "samples\tclade\n" ) p1.tmp p2.tmp > trait-tab-for-anvio.tsv

rm *.tmp
```

Importing:

```bash
anvi-import-misc-data -p Picocyano-RefSeq-61-pan/picocyanos-PAN.db --target-data-table layers trait-tab-for-anvio.tsv
```

Running functional enrichment test based on clade:

```bash
anvi-get-enriched-functions-per-pan-group -p Picocyano-RefSeq-61-pan/picocyanos-PAN.db \
                                          -g Picocyano-RefSeq-61-GENOMES.db --category clade \
                                          --annotation-source NCBI_PGAP \
                                          -o Marine-clade-enriched-functions.tsv \
                                          --functional-occurrence-table-output Function-occurrence-frequencies.tsv
```

---
> **Full table available for download by clicking [here](https://ndownloader.figshare.com/files/24723917).**
---

Perfectly split ones in this sub-table (there are lots more significant still in the full table):

**Functions only in all marine clade genomes**
| Function | Prop. marine clade | Prop. other |
|:---------|------:|------:|
| carboxysome shell carbonic anhydrase | 1 | 0 |
| 4-hydroxybenzoate polyprenyltransferase | 1 | 0 |
| riboflavin synthase | 1 | 0 |
| chorismate lyase | 1 | 0 |
| carotenoid oxygenase family protein | 1 | 0 |
| energy-coupling factor transporter transmembrane protein EcfT | 1 | 0 |
| carboxysome peptide B | 1 | 0 |
| SDR family NAD\(P\)-dependent oxidoreductase | 1 | 0 |
| photosystem II reaction center PsbP family protein | 1 | 0 |
| histidine phosphotransferase | 1 | 0 |
| DUF2834 domain-containing protein | 1 | 0 |
| DUF3136 domain-containing protein | 1 | 0 |
| carboxylating nicotinate-nucleotide diphosphorylase (nadC) | 1 | 0 |
| NAD\(P\)H-quinone oxidoreductase subunit L | 1 | 0 |
| phosphate ABC transporter substrate-binding protein PstS (pstS) | 1 | 0 |
| dTMP kinase | 1 | 0 |
| carboxysome peptide A | 1 | 0 |
| signal peptidase II | 1 | 0 |
| threonine ammonia-lyase, biosynthetic (ilvA) | 1 | 0 |
| glutamate racemase (murI) | 1 | 0 |
| aspartoacylase | 1 | 0 |
| S1 RNA-binding domain-containing protein | 1 | 0 |
| DUF4090 family protein | 1 | 0 |
| DUF1957 domain-containing protein | 1 | 0 |
| CIA30 family protein | 1 | 0 |
| Re/Si-specific NAD\(P\)(+) transhydrogenase subunit alpha | 1 | 0 |
| DUF1092 family protein | 1 | 0 |
| carboxysome shell protein | 1 | 0 |
| dihydroorotase (pyrC) | 1 | 0 |
| apolipoprotein N-acyltransferase | 1 | 0 |
| AIR synthase | 1 | 0 |
| sirohydrochlorin chelatase | 1 | 0 |
| potassium transporter TrkG | 1 | 0 |
| histidinol-phosphate aminotransferase family protein | 1 | 0 |
| LptF/LptG family permease | 1 | 0 |
| agmatinase (speB) | 1 | 0 |
| class I SAM-dependent RNA methyltransferase | 1 | 0 |
| Y-family DNA polymerase | 1 | 0 |
| GTP cyclohydrolase I | 1 | 0 |
| SufE family protein | 1 | 0 |
| nucleotide sugar dehydrogenase | 1 | 0 |
| glucose-6-phosphate dehydrogenase assembly protein OpcA | 1 | 0 |
| putative selenate ABC transporter substrate-binding protein | 1 | 0 |
| anthranilate synthase component I family protein | 1 | 0 |
| thylakoid membrane photosystem I accumulation factor | 1 | 0 |
| aminotransferase class IV | 1 | 0 |
| galactose mutarotase | 1 | 0 |
| repressor LexA (lexA) | 1 | 0 |
| YihY/virulence factor BrkB family protein | 1 | 0 |
| DUF3188 domain-containing protein | 1 | 0 |
| DoxX family protein | 1 | 0 |
| adenosine kinase | 1 | 0 |
| BMC domain-containing protein | 1 | 0 |
| ATP adenylyltransferase | 1 | 0 |
| DUF21 domain-containing protein | 1 | 0 |
| non-canonical purine NTP pyrophosphatase | 1 | 0 |
| Nif11-like leader peptide family natural product precursor | 1 | 0 |


**Functions only in non-marine clade 4 genomes**
| Function | Prop. marine clade | Prop. other |
|:---------|------:|------:|
| 4Fe-4S ferredoxin | 0 | 1 |
| energy-coupling factor ABC transporter ATP-binding protein | 0 | 1 |
| MlaE family lipid ABC transporter permease subunit | 0 | 1 |
| DUF2605 domain-containing protein | 0 | 1 |
| tetracycline resistance MFS efflux pump | 0 | 1 |
| gamma-glutamyltranspeptidase | 0 | 1 |
| alkene reductase (nemA) | 0 | 1 |
| HIT family protein | 0 | 1 |
| GTP 3',8-cyclase MoaA (moaA) | 0 | 1 |
| bidirectional hydrogenase complex protein HoxE (hoxE) | 0 | 1 |
| toxin-antitoxin system HicB family antitoxin | 0 | 1 |
| 3-methyladenine DNA glycosylase | 0 | 1 |
| lipoprotein signal peptidase | 0 | 1 |
| 5-(carboxyamino)imidazole ribonucleotide synthase (purK) | 0 | 1 |
| zinc-dependent alcohol dehydrogenase family protein | 0 | 1 |
| TrbI/VirB10 family protein | 0 | 1 |
| DUF896 domain-containing protein | 0 | 1 |
| AMIN domain-containing protein | 0 | 1 |
| coproporphyrinogen III oxidase | 0 | 1 |
| valine--pyruvate transaminase | 0 | 1 |
| histidinol-phosphate transaminase | 0 | 1 |
| polyphosphate:AMP phosphotransferase (pap) | 0 | 1 |
| flagellar biosynthesis protein FlgM | 0 | 1 |
| DUF1802 family protein | 0 | 1 |
| PilI type IV pilus biogenesis protein | 0 | 1 |
| Rpn family recombination-promoting nuclease/putative transposase | 0 | 1 |
| RNA-binding region RNP-1 | 0 | 1 |
| DUF4281 domain-containing protein | 0 | 1 |
| precorrin-3B synthase (cobG) | 0 | 1 |
| lipase | 0 | 1 |
| DUF2156 domain-containing protein | 0 | 1 |
| glycoside hydrolase family 57 protein | 0 | 1 |
| DUF1350 domain-containing protein | 0 | 1 |
| hydroxyneurosporene-O-methyltransferase | 0 | 1 |
| cytochrome C biogenesis protein CcdA | 0 | 1 |
| argonaute PAZ domain-containing protein | 0 | 1 |
| nitrogenase | 0 | 1 |
| flavin oxidoreductase | 0 | 1 |
| serine/threonine phosphatase | 0 | 1 |
| beta-glucosidase | 0 | 1 |
| iron-containing alcohol dehydrogenase family protein | 0 | 1 |
| DUF1853 family protein | 0 | 1 |
| TIGR02652 family protein | 0 | 1 |
| thiol:disulfide interchange protein TxlA | 0 | 1 |
| G-D-S-L family lipolytic protein | 0 | 1 |
| glyoxalase-like domain protein | 0 | 1 |
| acetate/propionate family kinase | 0 | 1 |
| glycosyl hydrolase family 15 | 0 | 1 |
| DUF4330 domain-containing protein | 0 | 1 |
| DUF87 domain-containing protein | 0 | 1 |
| M48 family metallopeptidase | 0 | 1 |
| cb-type cytochrome C oxidase subunit I | 0 | 1 |
| dihydrofolate reductase | 0 | 1 |
| cobalamin biosynthesis protein | 0 | 1 |
| heterocyst frequency control protein PatD (patD) | 0 | 1 |
| DUF98 domain-containing protein | 0 | 1 |
| DUF445 family protein | 0 | 1 |
| Vat family streptogramin A O-acetyltransferase (vat) | 0 | 1 |
| H(+)/Cl(-) exchange transporter ClcA (clcA) | 0 | 1 |
| DUF72 domain-containing protein | 0 | 1 |
| carboxylating nicotinate-nucleotide diphosphorylase | 0 | 1 |
| YbjN domain-containing protein | 0 | 1 |
| asparagine--tRNA ligase (asnS) | 0 | 1 |
| phosphate signaling complex protein PhoU (phoU) | 0 | 1 |
| serine/threonine protein kinase | 0 | 1 |
| glycogen debranching protein | 0 | 1 |
| YtxH domain-containing protein | 0 | 1 |
| DUF29 family protein | 0 | 1 |
| Zn(II)-sensing metalloregulatory transcriptional repressor SmtB (smtB) | 0 | 1 |
| DUF370 domain-containing protein | 0 | 1 |
| DUF3291 domain-containing protein | 0 | 1 |
| competence protein ComE | 0 | 1 |
| pantothenate metabolism flavoprotein | 0 | 1 |
| SLC26A/SulP transporter family protein | 0 | 1 |
| Fe(3+) ABC transporter substrate-binding protein | 0 | 1 |
| ArsA family ATPase | 0 | 1 |
| YjgP/YjgQ family permease | 0 | 1 |
| Y-family DNA polymerase (umuC) | 0 | 1 |
| peptide ABC transporter substrate-binding protein | 0 | 1 |
| multidrug ABC transporter permease | 0 | 1 |
| 8-amino-7-oxononanoate synthase | 0 | 1 |
| threonine-phosphate decarboxylase | 0 | 1 |
| ATP-dependent Zn protease | 0 | 1 |
| acyl-CoA dehydrogenase | 0 | 1 |
| threonylcarbamoyl-AMP synthase | 0 | 1 |
| carbon dioxide concentrating mechanism protein CcmL | 0 | 1 |
| curli assembly protein CsgF | 0 | 1 |
| HTH-type transcriptional activator CmpR | 0 | 1 |
| fibrillin | 0 | 1 |
| TolB family protein | 0 | 1 |
| lignostilbene alpha-beta-dioxygenase | 0 | 1 |
| caspase family protein | 0 | 1 |
| pre-16S rRNA-processing nuclease YqgF (ruvX) | 0 | 1 |
| TRC40/GET3/ArsA family transport-energizing ATPase | 0 | 1 |
| P-loop NTPase family protein | 0 | 1 |
| 1,4-dihydroxy-2-naphthoyl-CoA hydrolase | 0 | 1 |
| carboxymethylenebutenolidase | 0 | 1 |
| apolipoprotein N-acyltransferase (lnt) | 0 | 1 |
| K+ transporter Trk | 0 | 1 |
| proton-translocating NADH-quinone oxidoreductase subunit NdhD3 (ndhD3) | 0 | 1 |
| DUF4168 domain-containing protein | 0 | 1 |
| dTMP kinase (tmk) | 0 | 1 |
| 6-phosphofructokinase | 0 | 1 |
| VacB/RNase II family 3'-5' exoribonuclease | 0 | 1 |
| agmatine deiminase family protein | 0 | 1 |
| molybdopterin synthase sulfur carrier subunit MoaD (moaD) | 0 | 1 |
| S-layer homology domain-containing protein | 0 | 1 |
| DNA damage-inducible protein DinB | 0 | 1 |
| E3 ubiquitin ligase | 0 | 1 |
| nitrate/nitrite transporter NrtS (nrtS) | 0 | 1 |
| KH domain-containing protein | 0 | 1 |
| Gldg family protein | 0 | 1 |
| DUF2358 domain-containing protein | 0 | 1 |
| riboflavin synthase (ribE) | 0 | 1 |
| FAD-dependent hydroxylase | 0 | 1 |
| lipoprotein-like | 0 | 1 |
| OstA family protein | 0 | 1 |
| PstS family phosphate ABC transporter substrate-binding protein (pstS) | 0 | 1 |
| RNA polymerase sigma factor SigF | 0 | 1 |
| cytosine deaminase | 0 | 1 |
| PD40 domain-containing protein | 0 | 1 |
| bifunctional molybdenum cofactor biosynthesis protein MoaC/MoaB | 0 | 1 |
| phycobilisome maturation protein | 0 | 1 |
| DNA mismatch repair endonuclease MutL (mutL) | 0 | 1 |
| CbiQ family ECF transporter T component | 0 | 1 |
| RNA metabolism | 0 | 1 |
| DNA repair protein | 0 | 1 |
| ATP-dependent DNA helicase RecQ | 0 | 1 |
| SAM-dependent chlorinase/fluorinase | 0 | 1 |
| heme ABC exporter ATP-binding protein CcmA (ccmA) | 0 | 1 |
| GTP cyclohydrolase I FolE (folE) | 0 | 1 |
| type IV pilus assembly protein PilM (pilM) | 0 | 1 |
| Rubisco chaperonin | 0 | 1 |
| carbon dioxide-concentrating mechanism protein CcmK | 0 | 1 |
| chemotaxis protein CheA | 0 | 1 |
| purine-binding chemotaxis protein CheW | 0 | 1 |
| iron-regulated protein A | 0 | 1 |
| nicotinate phosphoribosyltransferase | 0 | 1 |
| LysE family translocator | 0 | 1 |
| class 1 fructose-bisphosphatase | 0 | 1 |
| 2-succinylbenzoate--CoA ligase | 0 | 1 |
| elongation factor G | 0 | 1 |
| glucose-6-phosphate dehydrogenase assembly protein OpcA (opcA) | 0 | 1 |
| SUF system NifU family Fe-S cluster assembly protein | 0 | 1 |
| type 2 isopentenyl-diphosphate Delta-isomerase | 0 | 1 |
| JAB domain-containing protein (radC) | 0 | 1 |
| molybdenum-pterin-binding protein | 0 | 1 |
| D-glycero-beta-D-manno-heptose-7-phosphate kinase | 0 | 1 |
| DUF4240 domain-containing protein | 0 | 1 |
| glutamate racemase | 0 | 1 |
| 3'(2'),5'-bisphosphate nucleotidase | 0 | 1 |
| allophycocyanin alpha chain-like protein | 0 | 1 |
| D-hexose-6-phosphate mutarotase | 0 | 1 |
| DNA double-strand break repair nuclease NurA | 0 | 1 |
| WYL domain-containing protein | 0 | 1 |
| ADP-ribosylglycohydrolase | 0 | 1 |
| branched chain amino acid aminotransferase | 0 | 1 |
| heavy-metal-associated domain-containing protein | 0 | 1 |
| riboflavin deaminase | 0 | 1 |
| homoserine O-succinyltransferase (metA) | 0 | 1 |
| carbon-nitrogen hydrolase | 0 | 1 |
| N-acetylglucosamine-6-phosphate deacetylase (nagA) | 0 | 1 |
| tRNA glutamyl-Q(34) synthetase GluQRS | 0 | 1 |
| NAD(P)H-quinone oxidoreductase subunit L (ndhL) | 0 | 1 |
| DUF3488 domain-containing transglutaminase family protein | 0 | 1 |
| FGGY-family carbohydrate kinase | 0 | 1 |
| carboxypeptidase | 0 | 1 |
| homocysteine biosynthesis protein | 0 | 1 |
| peptidase M42 | 0 | 1 |
| mercuric reductase | 0 | 1 |
| 5',5'-P-1P-4-tetraphosphate phosphorylase | 0 | 1 |
| tRNA-binding protein | 0 | 1 |
| cytochrome c biogenesis protein | 0 | 1 |
| heat-inducible transcriptional repressor HrcA (hrcA) | 0 | 1 |
| purple acid phosphatase | 0 | 1 |
| lipid-A-disaccharide synthase | 0 | 1 |
| Tab2/Atab2 family RNA-binding protein | 0 | 1 |
| DUF975 family protein | 0 | 1 |
| SMC family ATPase | 0 | 1 |
| bifunctional sterol desaturase/short chain dehydrogenase | 0 | 1 |
| MATE family efflux transporter | 0 | 1 |
| carbon dioxide concentrating mechanism protein CcmM | 0 | 1 |
| rod shape-determining protein MreD (mreD) | 0 | 1 |
| chemotaxis protein CheW | 0 | 1 |
| murein transglycosylase | 0 | 1 |
| PLP-dependent aminotransferase family protein | 0 | 1 |
| M1 family metallopeptidase | 0 | 1 |
| phycobilisome degradation protein NblA | 0 | 1 |
| DUF2281 domain-containing protein | 0 | 1 |
| 2-diaminobutyrate decarboxylase | 0 | 1 |
| endonuclease | 0 | 1 |
| lipid kinase | 0 | 1 |
| late competence development ComFB family protein | 0 | 1 |
| citramalate synthase (cimA) | 0 | 1 |
| cobalt-precorrin-6A reductase | 0 | 1 |
| folate/biopterin family MFS transporter (fbt) | 0 | 1 |
| TIGR04376 family protein | 0 | 1 |
| membrane protease subunit | 0 | 1 |
| 5-oxoprolinase | 0 | 1 |
| HypC/HybG/HupF family hydrogenase formation chaperone (hypC) | 0 | 1 |
| DUF3916 domain-containing protein | 0 | 1 |
| AAA-associated domain-containing protein | 0 | 1 |
| DUF4349 domain-containing protein | 0 | 1 |
| NADP-dependent malic enzyme | 0 | 1 |
| aquaporin family protein | 0 | 1 |
| 50S ribosome-binding GTPase | 0 | 1 |
| periplasmic heavy metal sensor | 0 | 1 |
| aminodeoxychorismate synthase component I (pabB) | 0 | 1 |

### Function counts
It is also worth looking at this table of annotation counts per genome. Lots of interesting things beyond presence/absence (which is what the above considers), such as large gene-copy number differences between the bugs (e.g. "response regulator", ABC transporter ATP-binding protein", "helix-turn-helix transcriptional regulator"). 

**Table can be downloaded by clicking [here](https://ndownloader.figshare.com/files/24725045).**


## Some gene trees
First looking for things that are in all genomes, but in two gene clusters (as that means they are annotated the same, but sequence-wise formed two clusters)

```bash
awk -F $'\t' ' { OFS=FS } { if ( $8 == 1 && $9 == 1 && length($7) == 24 ) print $1,$7 } ' Marine-clade-enriched-functions.tsv > potentially-interesting-functions.tmp
```

Want to work with things only in single-copy in each genome, so total number of found genes for the two gene clusters for each needs to by 61:

```bash
cut -f 2 potentially-interesting-functions.tmp | sed 's/, /\\|/' > GC-greps.tmp

for i in $(cat GC-greps.tmp)
do
    zgrep -w "$i" Picocyano-RefSeq-61-pan-summary/picocyanos_gene_clusters_summary.txt.gz | wc -l
    
done > counts.tmp

paste potentially-interesting-functions.tmp counts.tmp | awk -F $'\t' ' $3 == 61 ' | cut -f 1,2 > potentially-interesting-functions.tsv

rm *.tmp
```

Here's making a version that I could paste in here as a table:

```bash
  # note, sed tab replace as below won't work on darwin shell
sed 's/\t/ | /' potentially-interesting-functions.tsv | sed 's/^/| /' | sed 's/$/ |/'
```

> **NOTE**  
> Those in the table below that are linked are the ones I made trees for (just picking some). There are snapshots of them below if wanting to just scroll through them.

| Function | GC IDs |
|:---- |----|
| [divergent PAP2 family protein](https://itol.embl.de/tree/7184244121264731600483373) | GC_00000690, GC_00005009 |
| DUF3134 domain-containing protein | GC_00000626, GC_00005097 |
| UDP-N-acetylmuramoyl-L-alanyl-D-glutamate--2, 6-diaminopimelate ligase | GC_00000839, GC_00003650 |
| TIGR01777 family protein | GC_00000666, GC_00005347 |
| orotate phosphoribosyltransferase | GC_00000752, GC_00006492 |
| bifunctional diaminohydroxyphosphoribosylaminopyrimidine deaminase/5-amino-6-(5-phosphoribosylamino)uracil reductase RibD (ribD) | GC_00000734, GC_00006499 |
| glycine oxidase ThiO (thiO) | GC_00000751, GC_00005720 |
| quinone-dependent dihydroorotate dehydrogenase | GC_00000743, GC_00006501 |
| DUF3086 domain-containing protein | GC_00000753, GC_00005876 |
| F0F1 ATP synthase subunit delta | GC_00000662, GC_00005162 |
| tRNA (adenosine(37)-N6)-dimethylallyltransferase MiaA (miaA) | GC_00000740, GC_00006303 |
| competence/damage-inducible protein A | GC_00000847, GC_00003441 |
| [ferredoxin:protochlorophyllide reductase (ATP-dependent) subunit B](https://itol.embl.de/tree/7184244121426141600479628) | GC_00000733, GC_00005550 |
| homoserine kinase | GC_00000715, GC_00005313 |
| [ribosome maturation factor RimM (rimM)](https://itol.embl.de/tree/718424412194671600485122) | GC_00000842, GC_00003687 |
| low molecular weight phosphotyrosine protein phosphatase | GC_00000654, GC_00005105 |
| glutamate--cysteine ligase (gshA) | GC_00000774, GC_00004594 |
| adenosylhomocysteinase | GC_00000744, GC_00006430 |
| tRNA (cytidine(34)-2'-O)-methyltransferase | GC_00001120, GC_00002430 |
| [photosystem II reaction center protein T](https://itol.embl.de/tree/7184244121293351600483557) | GC_00000623, GC_00005194 |
| TIGR04168 family protein | GC_00000721, GC_00005887 |
| prephenate dehydratase (pheA) | GC_00000709, GC_00006240 |
| 7-carboxy-7-deazaguanine synthase QueE | GC_00000678, GC_00005819 |
| DUF3727 domain-containing protein | GC_00000695, GC_00005216 |
| ATP-dependent DNA helicase RecG (recG) | GC_00000831, GC_00003483 |
| preprotein translocase subunit SecG (secG) | GC_00000758, GC_00004397 |
| DUF2470 domain-containing protein | GC_00000820, GC_00003438 |
| protoheme IX farnesyltransferase | GC_00000739, GC_00005507 |
| [photosystem II reaction center protein J](https://itol.embl.de/tree/7184244121308301600483721) | GC_00000631, GC_00005635 |
| Holliday junction branch migration DNA helicase RuvB (ruvB) | GC_00000850, GC_00003534 |
| cytochrome b6-f complex subunit PetM (petM) | GC_00000625, GC_00005672 |
| DUF2499 domain-containing protein | GC_00000828, GC_00003247 |
| [DNA polymerase III subunit gamma/tau](https://itol.embl.de/tree/7184244121323361600483932) | GC_00000838, GC_00003383 |
| DUF2232 domain-containing protein | GC_00000712, GC_00005006 |
| [photosystem II reaction center X protein](https://itol.embl.de/tree/7184244121124951600482043) | GC_00000918, GC_00002588 |
| dihydroxy-acid dehydratase (ilvD) | GC_00000750, GC_00005388 |
| [iron ABC transporter permease](https://itol.embl.de/tree/7184244121347311600484072) | GC_00000706, GC_00004930 |
| ferredoxin:protochlorophyllide reductase (ATP-dependent) subunit N | GC_00000685, GC_00006200 |
| penicillin-binding protein 2 (mrdA) | GC_00000747, GC_00006111 |
| chorismate mutase (aroH) | GC_00000833, GC_00003295 |
| ribulose bisphosphate carboxylase small subunit | GC_00000621, GC_00005543 |
| uroporphyrinogen-III synthase | GC_00000688, GC_00005055 |
| DUF4346 domain-containing protein | GC_00000826, GC_00003481 |
| NAD\(P\)(+) transhydrogenase (Re/Si-specific) subunit beta | GC_00000644, GC_00005697 |
| TIGR01548 family HAD-type hydrolase | GC_00000735, GC_00005082 |
| [photosynthesis system II assembly factor Ycf48](https://itol.embl.de/tree/7184244121129421600485361) | GC_00000702, GC_00006034 |
| YbaB/EbfC family nucleoid-associated protein | GC_00000629, GC_00005086 |
| [photosystem II manganese-stabilizing polypeptide](https://itol.embl.de/tree/7184244121382861600484222) | GC_00001048, GC_00002498 |
| penicillin-binding protein 2 | GC_00000749, GC_00005532 |
| UDP-N-acetylmuramate dehydrogenase (murB) | GC_00000841, GC_00003436 |
| [50S ribosomal protein L29](https://itol.embl.de/tree/7184244121416951600484371) | GC_00000825, GC_00003386 |
| diacylglycerol kinase family protein | GC_00000835, GC_00003539 |
| tRNA pseudouridine(38-40) synthase TruA (truA) | GC_00000718, GC_00006514 |
| DUF3611 family protein | GC_00000670, GC_00005385 |
| endonuclease III (nth) | GC_00000855, GC_00003014 |
| [sulfate adenylyltransferase (sat)](https://itol.embl.de/tree/718424412129211600484782) | GC_00000728, GC_00005215 |
| [30S ribosome-binding factor RbfA (rbfA)](https://itol.embl.de/tree/7184244121444431600484491) | GC_00000823, GC_00003790 |
| DUF2854 domain-containing protein | GC_00000686, GC_00005768 |
| segregation/condensation protein A | GC_00000714, GC_00005802 |
| C-3',4' desaturase CrtD (crtD) | GC_00000720, GC_00005852 |
| dihydropteroate synthase (folP) | GC_00000829, GC_00003240 |
| [nickel pincer cofactor biosynthesis protein LarB (larB)](https://itol.embl.de/tree/7184244121481701600484648) | GC_00000646, GC_00006438 |
| bifunctional riboflavin kinase/FAD synthetase | GC_00000701, GC_00005487 |
| glutathione-disulfide reductase (gorA) | GC_00000844, GC_00003522 |
| glycine cleavage system aminomethyltransferase GcvT (gcvT) | GC_00000836, GC_00003379 |
| tRNA uridine-5-carboxymethylaminomethyl(34) synthesis GTPase MnmE (mnmE) | GC_00000849, GC_00003361 |
| UDP-3-O-acyl-N-acetylglucosamine deacetylase | GC_00000682, GC_00005964 |
| dethiobiotin synthase (bioD) | GC_00000642, GC_00005447 |
| DUF1816 domain-containing protein | GC_00000705, GC_00004942 |

Making an iToL mapping file to color all easily with one of my `bit` programs:

```bash
bit-gen-iToL-map -g non-marine-clade-labels.txt -w labels -c green -o iToL-cols-p1.txt
bit-gen-iToL-map -g marine-clade-labels.txt -w labels -c blue -o iToL-cols-p2.txt

cat iToL-cols-p1.txt <( tail -n +5 iToL-cols-p2.txt ) > iToL-cols.txt

rm iToL-cols-p*
```

Ended up modifying the colors anyway, here's the file if i want it elsewhere or lose it:

```
TREE_COLORS
SEPARATOR TAB
DATA

GCF_000010065.1_Synechococcus_elongatus_PCC_6301	label	#2a7926	bold
GCF_000012525.1_Synechococcus_elongatus_PCC_7942	label	#2a7926	bold
GCF_000817325.1_Synechococcus_sp_UTEX_2973_UTEX_2973	label	#2a7926	bold
GCF_003957805.1_Synechococcus_elongatus_UTEX_3055	label	#2a7926	bold
GCF_000018065.1_Prochlorococcus_marinus_MIT_9215	label	#b10c34	bold
GCF_000757845.1_Prochlorococcus_sp_MIT_0604_MIT_0604	label	#b10c34	bold
GCF_000015645.1_Prochlorococcus_marinus_AS9601	label	#b10c34	bold
GCF_000015965.1_Prochlorococcus_marinus_MIT_9301	label	#b10c34	bold
GCF_000015665.1_Prochlorococcus_marinus_MIT_9515	label	#b10c34	bold
GCF_000011465.1_Prochlorococcus_marinus_MED4	label	#b10c34	bold
GCF_000012645.1_Prochlorococcus_marinus_MIT_9312	label	#b10c34	bold
GCF_000063525.1_Synechococcus_sp_RCC307	label	#434da7	bold
GCF_014279855.1_Synechococcus_sp_Minos11_MINOS11	label	#434da7	bold
GCF_001885215.1_Synechococcus_sp_SynAce01_SynAce01	label	#434da7	bold
GCF_008807075.1_Synechococcus_sp_RSCCF101_RSCCF101	label	#434da7	bold
GCF_014280235.1_Cyanobium_sp_NS01_NS01	label	#8c712d	bold
GCF_900088535.1_Cyanobium_sp_NIES-981_NIES-981	label	#8c712d	bold
GCF_000316515.1_Cyanobium_gracile_PCC_6307	label	#8c712d	bold
GCF_000179235.2_Synechococcus_sp_CB0101_CB0101	label	#434da7	bold
GCF_014217875.1_Synechococcus_sp_LTW-R_LTW-R	label	#434da7	bold
GCF_000007925.1_Prochlorococcus_marinus_CCMP1375_SS120	label	#b10c34	bold
GCF_000757865.1_Prochlorococcus_sp_MIT_0801_MIT_0801	label	#b10c34	bold
GCF_000012465.1_Prochlorococcus_marinus_NATL2A	label	#b10c34	bold
GCF_000015685.1_Prochlorococcus_marinus_NATL1A	label	#b10c34	bold
GCF_000011485.1_Prochlorococcus_marinus_MIT9313	label	#b10c34	bold
GCF_004209775.1_Synechococcus_sp_WH_8101_WH_8101	label	#434da7	bold
GCF_014304795.1_Synechococcus_sp_WH_8101_WH_8101	label	#434da7	bold
GCF_014279595.1_Synechococcus_sp_RS9909_RS9909	label	#434da7	bold
GCF_000012505.1_Synechococcus_sp_CC9902_CC9902	label	#434da7	bold
GCF_014279555.1_Synechococcus_sp_SYN20_SYN20	label	#434da7	bold
GCF_014279775.1_Synechococcus_sp_PROS-9-1_PROS-9-1	label	#434da7	bold
GCF_014279835.1_Synechococcus_sp_MVIR-18-1_MVIR-18-1	label	#434da7	bold
GCF_000014585.1_Synechococcus_sp_CC9311_CC9311	label	#434da7	bold
GCF_014279655.1_Synechococcus_sp_ROS8604_ROS8604	label	#434da7	bold
GCF_000063505.1_Synechococcus_sp_WH_7803	label	#434da7	bold
GCF_014279955.1_Synechococcus_sp_BMK-MC-1_BMK-MC-1	label	#434da7	bold
GCF_014279795.1_Synechococcus_sp_PROS-7-1_PROS-7-1	label	#434da7	bold
GCF_014279875.1_Synechococcus_sp_MEDNS5_MEDNS5	label	#434da7	bold
GCF_000737535.1_Synechococcus_sp_KORDI-100_KORDI_100	label	#434da7	bold
GCF_000737575.1_Synechococcus_sp_KORDI-49_KORDI_49	label	#434da7	bold
GCF_014280215.1_Synechococcus_sp_A15-127_A15-127	label	#434da7	bold
GCF_000012625.1_Synechococcus_sp_CC9605_CC9605	label	#434da7	bold
GCF_000737595.1_Synechococcus_sp_KORDI-52_KORDI_52	label	#434da7	bold
GCF_014280075.1_Synechococcus_sp_A15-62_A15-62	label	#434da7	bold
GCF_014279895.1_Synechococcus_sp_M16_1_M16.1	label	#434da7	bold
GCF_014279635.1_Synechococcus_sp_RS9902_RS9902	label	#434da7	bold
GCF_014279535.1_Synechococcus_sp_TAK9802_TAK9802	label	#434da7	bold
GCF_014280115.1_Synechococcus_sp_A15-44_A15-44	label	#434da7	bold
GCF_000161795.2_Synechococcus_sp_WH_8109_WH_8109	label	#434da7	bold
GCF_014279755.1_Synechococcus_sp_PROS-U-1_PROS-U-1	label	#434da7	bold
GCF_014280175.1_Synechococcus_sp_A15-28_A15-28	label	#434da7	bold
GCF_000195975.1_Synechococcus_sp_WH_8102_WH_8102	label	#434da7	bold
GCF_001182765.1_Synechococcus_sp_WH_8103	label	#434da7	bold
GCF_014280015.1_Synechococcus_sp_A18-46_1_A18-46.1	label	#434da7	bold
GCF_014280055.1_Synechococcus_sp_A18-40_A18-40	label	#434da7	bold
GCF_014280195.1_Synechococcus_sp_A15-24_A15-24	label	#434da7	bold
GCF_014279975.1_Synechococcus_sp_BIOS-U3-1_BIOS-U3-1	label	#434da7	bold
GCF_014304815.1_Synechococcus_sp_MIT_S9220_MIT_S9220	label	#434da7	bold
GCF_014279815.1_Synechococcus_sp_NOUM97013_NOUM97013	label	#434da7	bold
GCF_014280035.1_Synechococcus_sp_A18-25c_A18-25c	label	#434da7	bold
GCF_014280095.1_Synechococcus_sp_A15-60_A15-60	label	#434da7	bold
```


**Script to automate generating these specific single-gene trees**

```bash
cat ad-hoc-make-gene-tree.sh
```

```bash
  # getting wanted gene AA seqs from each genome:
zgrep -w "${1}\|${2}" Picocyano-RefSeq-61-pan-summary/picocyanos_gene_clusters_summary.txt.gz | \
      awk -F $'\t' ' BEGIN { OFS="\n" } { print ">"$4,$15 } ' | tr -d "-" > ${3}.faa

  # swapping labels so they have tax info instead of just accession
grep ">" ${3}.faa | tr -d ">" | cut -f 1,2 -d "_" > ids.tmp

grep -v ">" ${3}.faa > seqs.tmp

for i in $(cat ids.tmp)
do
    grep -m 1 "^$i" target-picos-gtotree-standard/Genomes_summary_info.tsv | cut -f 2 | sed 's/^/>/'
    
done > new-ids.tmp

paste -d "\n" new-ids.tmp seqs.tmp > ${3}.faa

rm ids.tmp seqs.tmp new-ids.tmp

  # aligning
muscle -in ${3}.faa -out ${3}_aln.faa

  # treeing
fasttree ${3}_aln.faa > ${3}.tre
```

Just picking some from that table to peek at...

### divergent PAP2 family protein

```bash
bash ad-hoc-make-gene-tree.sh GC_00000690 GC_00005009 divergent_PAP2_family_protein
```

<a href="https://i.imgur.com/P33iQ1H.png"><img src="https://i.imgur.com/P33iQ1H.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121264731600483373).**


### ferredoxin:protochlorophyllide reductase (ATP-dependent) subunit B

```bash
bash ad-hoc-make-gene-tree.sh GC_00000733 GC_00005550 ferredoxin_protochlorophyllide_reductace_subunit_B
```

<a href="https://i.imgur.com/jU8zDF2.png"><img src="https://i.imgur.com/jU8zDF2.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121426141600479628).**


### ribosome maturation factor RimM (rimM)

```bash
bash ad-hoc-make-gene-tree.sh GC_00000842 GC_00003687 ribosome_maturation_factor_RimM
```

<a href="https://i.imgur.com/2srtd6b.png"><img src="https://i.imgur.com/2srtd6b.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/718424412194671600485122).**


### photosystem II reaction center protein T

```bash
bash ad-hoc-make-gene-tree.sh GC_00000623 GC_00005194 photosystem_II_reaction_center_protein_T
```

<a href="https://i.imgur.com/MdoSxs1.png"><img src="https://i.imgur.com/MdoSxs1.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121293351600483557).**


### photosystem II reaction center protein J

```bash
bash ad-hoc-make-gene-tree.sh GC_00000631 GC_00005635 photosystem_II_reaction_center_protein_J
```

<a href="https://i.imgur.com/lAIMa6m.png"><img src="https://i.imgur.com/lAIMa6m.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121308301600483721).**


### photosystem II reaction center X protein

```bash
bash ad-hoc-make-gene-tree.sh GC_00000918 GC_00002588 photosystem_II_reaction_center_X
```

<a href="https://i.imgur.com/Puwh8B2.png"><img src="https://i.imgur.com/Puwh8B2.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121124951600482043).**


### photosynthesis system II assembly factor Ycf48

```bash
bash ad-hoc-make-gene-tree.sh GC_00000702 GC_00006034 photosynthesis_system_II_assembly_factor_Ycf48
```

<a href="https://i.imgur.com/6eX2ZOv.png"><img src="https://i.imgur.com/6eX2ZOv.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121129421600485361).**


### photosystem II manganese-stabilizing polypeptide

```bash
bash ad-hoc-make-gene-tree.sh GC_00001048 GC_00002498 photosystem_II_manganese-stabilizing_polypeptide
```

<a href="https://i.imgur.com/yesbIar.png"><img src="https://i.imgur.com/yesbIar.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121382861600484222).**


### iron ABC transporter permease

```bash
bash ad-hoc-make-gene-tree.sh GC_00000706 GC_00004930 iron_ABC_transporter_permease
```

<a href="https://i.imgur.com/r7bVX0a.png"><img src="https://i.imgur.com/r7bVX0a.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121347311600484072).**


### DNA polymerase III subunit gamma/tau

```bash
bash ad-hoc-make-gene-tree.sh GC_00000838 GC_00003383 DNA_polymerase_III_subunit_gamma_tau
```

<a href="https://i.imgur.com/nkFAD7O.png"><img src="https://i.imgur.com/nkFAD7O.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121323361600483932).**


### 50S ribosomal protein L29

```bash
bash ad-hoc-make-gene-tree.sh GC_00000825 GC_00003386 50S_ribosomal_protein_L29
```

<a href="https://i.imgur.com/WpmvnOp.png"><img src="https://i.imgur.com/WpmvnOp.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121416951600484371).**


### 30S ribosome-binding factor RbfA (rbfA)

```bash
bash ad-hoc-make-gene-tree.sh GC_00000823 GC_00003790 30S_ribosome-binding_factor_RbfA
```

<a href="https://i.imgur.com/kxKI2Rf.png"><img src="https://i.imgur.com/kxKI2Rf.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121444431600484491).**


### sulfate adenylyltransferase (sat)

```bash
bash ad-hoc-make-gene-tree.sh GC_00000728 GC_00005215 sulfate_adenylyltransferase
```

<a href="https://i.imgur.com/OIUdBwa.png"><img src="https://i.imgur.com/OIUdBwa.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/718424412129211600484782).**


### nickel pincer cofactor biosynthesis protein LarB (larB)

```bash
bash ad-hoc-make-gene-tree.sh GC_00000646 GC_00006438 nickel_pincer_cofactor_biosynthesis_protein_LarB
```

<a href="https://i.imgur.com/63XJlLT.png"><img src="https://i.imgur.com/63XJlLT.png"></a>

**Interactive tree [here](https://itol.embl.de/tree/7184244121481701600484648).**


<!-- ### 

```bash
bash ad-hoc-make-gene-tree.sh 
```

<a href=""><img src=""></a>

**Interactive tree [here]().**
 -->



