Diversity and dynamics within oral bacterial populations
Genome diversity & dynamcis of the oral microbiome
Background
We set out to investigate how genomic diversity changed within oral bacterial populations in healthy mouths (i.e., absent from relatively drastic selective events like antibiotics, hospitalization, etc.). We used the tongue dorsum as a model oral microbial system as it represents a rich oral community on a shedding surface intimately connected at the interface of the human body and exogenous perturbations. This document describes the full computaitonal methods necessary to reproduce the results described in the biorxiv preprint.
Briefly, the sequence data produced here was obtained by recruiting 17 healthy donors who self-sampled their tongue dorsa with a plastic tongue over a two-week window on days 1, 2, 5, 7, and 14. A few volunters missed a day or sampled on an alternate day (e.g., day 3 instead of day 2); all told, we collected 81 samples. For each sample we prepared a metagenomic library with Kapa HyperPrep kits, and the pooled library of all 81 samples was sequenced on 10 NextSeq HO PE150 runs.
Scripts included here that were run on the cluster maintain the requested cluster parameters (e.g., memory and time).
Code background
Raw reads were sequenced in 2017 at the Harvard Bauer Sequencing Core.
Each lane was in a separate directory, e.g., run1 and run2 each
contained 81 fastq pairs. Note that these included substantial human
DNA; for participant privacy we excluded human DNA from the raw data
uploaded to NCBI. For each run, the following script was run:
De-human
#!/bin/bash
mkdir deHumanized
# # cycle through all
fastas=$(ls A* | grep 'A*.fastq')
#fastas="AT15-TD_S58.R1.fastq.gz" # # re-do specific files
for file in $(echo "$fastas" | sed 's/.R[12].fastq.gz//' | uniq); do
samp=$(echo "$file" | sed 's/-.*$//g')
R1=$(echo "$file.R1.fastq.gz")
R2=$(echo "$file.R2.fastq.gz")
echo "#!/bin/bash
#SBATCH -N 1 #node
#SBATCH -n 1 #one cpu
#SBATCH --mem=8000
#SBATCH -t 0-6:00:00
#SBATCH -p shared
#SBATCH --job-name=deHum-$samp
#SBATCH -o odyssey_deHuman-$samp.out
#SBATCH -e odyssey_deHuman-$samp.err
#SBATCH --mail-type=END
#SBATCH --mail-user=dutter@g.harvard.edu
module load bib
# map vs human database, sam to bam
bowtie2 -x /n/holyscratch01/cavanaugh_lab/dutter/catalyst_metagenome_raw/hg38_bw2/hg38_noalt -1 $R1 -2 $R2 -S $file-mapped_and_unmapped.sam
samtools view -bS $file-mapped_and_unmapped.sam > $file-mapped_and_unmapped.bam
# get unmapped reads (non-human) to new bam
samtools view -b -f 12 -F 256 $file-mapped_and_unmapped.bam > $file-bothEndsUnmapped.bam
# bam back to fasta
bamtools sort -in $file-bothEndsUnmapped.bam -byname | bamtools convert -format fastq -out deHumanized/$file-dehuman.fastq
rm $file*bam $file*sam
" > $file-dehuman.sh
sbatch $file-dehuman.sh && rm $file-dehuman.sh
sleep 20
done
To keep the working directory clear, the original reads were then put
away with mkdir raw_reads && mv *fastq.gz raw_reads
This produced fastq files of all reads not aligning to the human genome. However, the resultant fastq files were interleaved and so we manually split them up with the following script from inside the directory containing each run’s reads:
#!/bin/bash
#SBATCH -N 1 #node
#SBATCH -n 1 #one cpu
#SBATCH --mem=12000
#SBATCH -t 0-2:00:00
#SBATCH -p general
#SBATCH --job-name="splitPE"
#SBATCH -o odyssey_split.out
#SBATCH -e odyssey_split.err
#SBATCH --mail-type=END
#SBATCH --mail-user=dutter@g.harvard.edu
for fasta in $(ls deHumanized/*-dehuman.fastq); do
prefix=$(echo "$fasta" | sed 's/-dehuman.fastq//')
grep -A 3 '^@N.*/1$' $fasta | grep -v '^--$' > $prefix.R1.fastq
grep -A 3 '^@N.*/2$' $fasta | grep -v '^--$' > $prefix.R2.fastq
done
These raw non-human reads were deposited in NCBI at SRA Accessions TO ADD upon manuscript publication.
The files have prefixes in the format IdDay that specify the donor’s
unique id, e.g. AE, and the sample day, e.g. 1. So sample
AB15_QUALITY-PASSED_R1.fastq.gz is the sample obtained from donor AB
on day 15.
QC
Reads were then filtered according to the Minoche et al. 2011 recommendations with illumina-utils:
#!/bin/bash
#SBATCH -N 1 #node
#SBATCH -n 1 #one cpu
#SBATCH --mem=12000
#SBATCH -t 5-06:00:00
#SBATCH -p general
#SBATCH --job-name="iu-filter"
#SBATCH -o odyssey_filter.out
#SBATCH -e odyssey_filter.err
#SBATCH --mail-type=END
#SBATCH --mail-user=dutter@g.harvard.edu
source ~/virtual-envs/anvio-dev-github/bin/activate
mkdir filtered-minoche
# get read names
R1=$(ls A*gz | grep 'R1.fastq' | sed 's/^/, /' | tr -d '\n' | sed 's/^, //')
R2=$(ls A*gz | grep 'R2.fastq' | sed 's/^/, /' | tr -d '\n' | sed 's/^, //')
# make config file
echo "[general]
project_name = catalyst_TD
researcher_email = dutter@g.harvard.edu
input_directory = $PWD
output_directory = $PWD/filtered-minoche
[files]
pair_1 = $R1
pair_2 = $R2
" > filter_config.txt
# filter reads
iu-filter-quality-minoche filter_config.txt
Then we moved back into the top-level directory, prepared a directory
with mkdir filtered-dehuman-reads-all-bzip, and finally concatenated
reads from different runs for the same sample with:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 12 # 1 cores from each
#SBATCH --contiguous
#SBATCH --mem=18G #per node
#SBATCH -t 0-12:00:00
#SBATCH -p shared
#SBATCH --job-name='concat'
#SBATCH -o odyssey_concat.out
#SBATCH -e odyssey_concat.err
#SBATCH --mail-type=NONE
for fastq in $(ls run1/deHumanized/*fastq | sed 's,^.*/,,g'); do
cat run*/deHumanized/$fastq > filtered-dehuman-reads-all-bzip/$fastq
done
MAG generation
To construct reference genomes for the bacterial populations dwelling in these populations, we used IDBA-UD to construct Metagenome Assembled Genomes (MAGs).
We co-assembled MAGs by mouth, that is, all metagenomes from one donor (4-7 timepoints) were combined and used to produce a single assembly. This approach was taken to strike a balance between assembling each sample independently, which may miss some low-abundance oral bacteria, and overloading the assembler with the maximum complexity of multiple closely related popualtions from different people’s mouths.
Co-assembly (IDBA-UD)
We used IDBA-UD as we were intially interested in being able to use the same assembler to co-assemble individual donors and also co-assemble the whole dataset. The whole dataset co-assembly did not end up being retained as part of the analysis, but at the time factored into our choice for assembly. At the time of performing these steps in 2017, the best two assemblers that could run on the entire dataset given the 500GB memory accessible were megahit and IDBA-UD (metaSPAdes was attempted but required too much memory). As IDBA-UD outperformed megahit empirically based on size and subsequent MAG quality (data not included), we chose the IDBA-UD assemblies to analyze in detail.
This scripts, and most of the future scripts in the MAG curation
workflow, start with a variable indiv that specifies the individual
being studied. For the rest of the MAG workflow we will use AE as an
example, however, this entire assembly and MAG refinement process was
performed separately for the other 16 donors as well.
So for donor AE, we made a separate directory like AE_coassembly
inside the top-level directory. Here is the script for generating the
co-assembly for donor AE:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 12 # 1 cores from each
#SBATCH --contiguous
#SBATCH --mem=80G #per node
#SBATCH -t 5-00:00:00
#SBATCH -p shared
#SBATCH --job-name='concat'
#SBATCH -o odyssey_concat.out
#SBATCH -e odyssey_concat.err
#SBATCH --mail-type=NONE
source ~/virtual-envs/anvio-dev-github/bin/activate
indiv="AE"
cat ../filtered-dehuman-reads-all-bzip/$indiv*QUALITY_PASSED_R1.fastq > raw_$indiv-R1.fastq
cat ../filtered-dehuman-reads-all-bzip/$indiv*QUALITY_PASSED_R2.fastq > raw_$indiv-R2.fastq
fq2fa --merge --filter raw_$indiv-R1.fastq raw_$indiv-R2.fastq interleaved_$indiv.fa
idba_ud -r interleaved_$indiv.fa -o idba_ud-$indiv
Generating contigs database
Contigs were then read into anvi’o
for management during subsequent binning, curation, and annotation of
MAGs. We kept only contigs with a minimum length (--min-len flag) of
2500 nucleotides (nt) as shorter contig lengths may introduce problems
for binning metrics like tetranucleotide frequency and differential
coverage. We generated contig databses with this script:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 1 # 1 cores from each
#SBATCH --contiguous
#SBATCH --mem=16G #per node
#SBATCH -t 0-24:00:00
#SBATCH -p shared
#SBATCH --job-name="contigsDB"
#SBATCH -o odyssey_cont.out
#SBATCH -e odyssey_cont.err
#SBATCH --mail-type=END
# need to activate venv directly because node receiving job order doesn't play with bashrc aliases
source ~/virtual-envs/anvio-dev-github/bin/activate
indiv="AE"
# clean up fasta deflines in contig file at the start for smooth downstream
# min-len 2500 = only keep contigs >= 2.5kb since 1kb gives poor tetranuc freq results
anvi-script-reformat-fasta $indiv-idba_ud/contig.fa --simplify-names -o $indiv-contigs.fa
anvi-script-reformat-fasta $indiv-contigs.fa --min-len 2500 -o $indiv-contigs-2500nt.fa
# do the contig gen
anvi-gen-contigs-database -f $indiv-contigs-2500nt.fa -o $indiv-CONTIGS.db
During this step, genes were also called with Prodigal as part of the
anvi-gen-contigs-database program.
And identified single-copy genes with anvi’o’s HMM models:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 10 # 10 cores from each
#SBATCH --contiguous
#SBATCH --mem=20G #per node
#SBATCH -t 0-06:00:00
#SBATCH -p shared
#SBATCH --job-name="hmm"
#SBATCH -o odyssey_hmm.out
#SBATCH -e odyssey_hmm.err
#SBATCH --mail-type=END
indiv="AE"
# need to activate venv directly because node receiving job order doesn't play with bashrc aliases
source ~/virtual-envs/anvio-dev-github/bin/activate
# identify single copy genes
anvi-run-hmms -c $indiv-CONTIGS.db --num-threads 10
Functional annotation
To aid in binning, we annotated quickly with only eggNOG at the start:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 12 # 12 cores from each
#SBATCH --contiguous
#SBATCH --mem=8G #per node
#SBATCH -t 0-24:00:00
#SBATCH -p shared
#SBATCH --job-name="eggnog"
#SBATCH -o odyssey_anviAnnotate.out
#SBATCH -e odyssey_anviAnnotate.err
#SBATCH --mail-type=END
# 20171016 - assign gene fams by interproscan
indiv="AE"
# need to activate venv directly because node receiving job order doesn't play with bashrc aliases
source ~/virtual-envs/anvio-dev-github/bin/activate
#export sequences
anvi-get-sequences-for-gene-calls -c $indiv-CONTIGS.db -o gene-calls-aa.faa --get-aa-sequences --wrap 0
module load python/2.7.14-fasrc02
/n/helmod/apps/centos7/Core/Anaconda/5.0.1-fasrc02/x/bin/python ../../eggnog/eggnog-mapper/emapper.py -i $batch --output ${batch}_eggnog -m diamond --cpu 12
Subsequent steps (mapping, profiling, and merging() take significant time, so we generally did not parallelize this step as it could just run in the background, and only became relevant once merging was done. However, a few co-assemblies had so many genes that
Mapping short reads to assembled contigs
To obtain the coverage of each of the assembled contigs in each sample
used for co-assembly, we used
bowtie2 to map short reads
from each metagenome onto a reference database corresponding to the
contigs in AE-CONTIGS.db:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 1 # 1 cores from each
#SBATCH --contiguous
#SBATCH --mem=2G #per node
#SBATCH -t 0-05:00:00
#SBATCH -p shared
#SBATCH --job-name="bt-mas-AE"
#SBATCH -o odyssey_map_master.out
#SBATCH -e odyssey_map_master.err
#SBATCH --mail-type=END
module load samtools bowtie2 zlib xz
indiv="AE"
outdir="bt_mapped"
# make the db
#bowtie2-build $indiv.fa $indiv
# get list of sample-runs, chop into per sample batches
batches=$(ls ../filtered-dehuman-reads-all-bzip/$indiv*QUALITY_PASSED_R1.fastq.gz | sed "s,^.*/,,g; s/_.*$//" | uniq)
# make dir to stuff samples into
mkdir $outdir
# iterate for each participant
for batch in $batches; do
echo "#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 12 # 12 cores from each
#SBATCH --contiguous
#SBATCH --mem=42G #per node
#SBATCH -t 1-00:00:00
#SBATCH -p shared
#SBATCH --job-name='$batch-bt2'
#SBATCH -o odyssey_map_$batch.out
#SBATCH -e odyssey_map_$batch.err
#SBATCH --mail-type=END
module load samtools bowtie2 zlib xz
batch='$batch'
samples=\$(ls ../../filtered-dehuman-reads-all-bzip/${batch}_*QUALITY_PASSED_R1.fastq.gz | sed 's/_R1.fastq.gz//')
for samp in \$samples; do
sampOut=\$(echo \"\$samp\" | sed \"s/^.*$batch/$outdir\\/$batch/; s/-.*$//\")
R1=\${samp}_R1.fastq.gz
R2=\${samp}_R2.fastq.gz
# map sample set
bowtie2 -x $indiv -1 \$R1 -2 \$R2 --no-unal -S \$sampOut.sam --threads 12
samtools view -S -b \$sampOut.sam > \$sampOut-RAW.bam
samtools sort -@ 12 \$sampOut-RAW.bam > \$sampOut.bam
samtools index \$sampOut.bam
# clean as you go
rm \$sampOut.sam \$sampOut-RAW.bam
done
" > slurm_$batch.sh
sbatch slurm_$batch.sh #&& rm slurm_$batch.sh
sleep 10
done
Profiling
The information (coverage, nucleotide variants) in the BAM files were then imported into anvi’o as PROFILE databases:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 1 # 1 cores from each
#SBATCH --contiguous
#SBATCH --mem=1G #per node
#SBATCH -t 0-05:00:00
#SBATCH -p shared
#SBATCH --job-name="prof-m-AE"
#SBATCH -o odyssey_anviProf-master.out
#SBATCH -e odyssey_anviProf-master.err
#SBATCH --mail-type=END
indiv="AE"
assembler="idba_ud"
mapdir="bt_mapped"
mkdir profs_indiv
# get list of samples
batches=$(ls $mapdir | sed 's/_.*$//' | sort | uniq)
for batch in $batches; do
#if ls profs_indiv | grep -q -v "^$batch$"; then
echo "#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 16 # cores from each
#SBATCH --contiguous
#SBATCH --mem=42G #per node
#SBATCH -t 1-24:00:00
#SBATCH -p shared
#SBATCH --job-name=\"prof-$batch\"
#SBATCH -o odyssey_anviProf-$batch.out
#SBATCH -e odyssey_anviProf-$batch.err
#SBATCH --mail-type=END
source ~/virtual-envs/anvio-dev-github/bin/activate
module load samtools zlib xz
bams=\"$mapdir/${batch}_*bam\"
# need to merge the multiple runs into single bam to pass to anvio
samtools merge $mapdir/${batch}.bam \$bams
samtools index $mapdir/${batch}.bam
echo \"Done merging $batch bams, starting profiling...\"
#rm -rf $sample
# profile it
anvi-profile -c $indiv-CONTIGS.db --profile-SCVs -i $mapdir/${batch}.bam --num-threads 16 --write-buffer-size 500 -o profs_
indiv/$batch
echo \"Done profiling $batch\"
" > prof-$batch.sh
sbatch prof-$batch.sh # && rm prof-$batch.sh
sleep 10
#fi
done
These individual per-metagenome PROFILE databases were then merged into a single MERGED database:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 1 # 1 cores from each
#SBATCH --contiguous
#SBATCH --mem=120G #per node
#SBATCH -t 3-00:00:00
#SBATCH -p shared
#SBATCH --job-name="merge-AE"
#SBATCH -o odyssey_anviMerge.out
#SBATCH -e odyssey_anviMerge.err
#SBATCH --mail-type=END
# 20170918 merge each individually profiled sample
# need to activate venv directly because node receiving job order doesn't play with bashrc aliases
source ~/virtual-envs/anvio-dev-venv/bin/activate
indiv="AE"
# import function annotation as should be done by now
anvi-import-functions -c $indiv-CONTIGS.db -i interpro-results-formatted.tsv
# merge
anvi-merge profs_indiv/A*/PROFILE.db --enforce-hierarchical-clustering --skip-concoct-binning -o $indiv-MERGED -c $indiv-CONTIGS.db
This step took quite a bit of memory. Significant improvements to anvi’o have been made since by the anvi’o team & community, hopefully if someone is replicating this today the newfound wisdom and grace of anvi’o will make the process more gentle :)
Binning & MAGs
The entire MAG generation process was divided into multiple steps,
detailed below but summarized here to give an idea of how they fit
together. First, CONCOCT did automatic binning of contigs to generate
very rough bins. These bins were then manually refined with
anvi-refine to split bins and identify and remove contaminant contigs.
Then, we manually refined again to do a second pass of refining, just
to be safe. MAGs were then defined from this set as bins with <10%
redundancy and either ≥80% completion or ≥2 Mbp. These MAGs were merged
across co-assemblies (17 sets of MAGs from 17 independent co-assemblies
-> one set of MAGs). Redundant MAGs (duplicate MAGs reconstructed the
same bacteria in different mouths) were then identified, and removed to
identify a final set of 390 non-redundant MAGs.
Binning 1: Rough pass with CONCOCT
At the time of binning these samples, anvi’o had a now-obsolete program
anvi-cluster-with-concoct. In the spirit of transparency, this was the
script that was actually run:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 1 # 1 cores from each
#SBATCH --contiguous
#SBATCH --mem=80G #per node
#SBATCH -t 3-00:00:00
#SBATCH -p shared
#SBATCH --job-name="concoct-AE"
#SBATCH -o odyssey_anviMerge.out
#SBATCH -e odyssey_anviMerge.err
#SBATCH --mail-type=END
# 20170918 merge each individually profiled sample
# need to activate venv directly because node receiving job order doesn't play with bashrc aliases
source ~/virtual-envs/anvio-dev-venv/bin/activate
indiv="AE"
# CONCOCT clustering
anvi-cluster-with-concoct -p $indiv-MERGED/PROFILE.db -c $indiv-CONTIGS.db --num-clusters 50 -C CONCOCT
Curently, as of May 2021, the modern anvi’o program would be
anvi-cluster-contigs with the --driver CONCOCT flag.
At this point, we manually refined bins on the basis of tetranucleotide
frequency, differential coverage, and single-nucleotide variant (SNV)
frequencies across all metagenomes corresponding to that co-assembly
(4-7 metagenomes per donor or co-assembly). We did this on a bin-by-bin
process with the anvi-refine feature.
Binning 2: Manual refinement
Having done a first-pass binning (across all donors), we then went through for a second pass. To roughly prune out the obviously poor MAGs of which there was no hope of saving (e.g., MAGs that were very small or very incomplete, regardless of potential contamination), we ran the following script inspired by Tom Delmont & A. Murat Eren’s script for their TARA Oceans analysis:
#!/bin/bash
for SET in `ls | grep "A[A-Z]-idba_ud"`
do
indiv=$(echo "$SET" | sed 's/-.*$//')
anvi-rename-bins -c $SET/$indiv-CONTIGS.db \
-p $SET/$indiv-MERGED/PROFILE.db \
--collection-to-read CONCOCT \
--collection-to-write ROUGH \
--call-MAGs \
--size-for-MAG 2 \
--min-completion-for-MAG 70 \
--max-redundancy-for-MAG 100 \
--prefix ROUGH_$indiv \
--report-file $SET/$indiv-rough-bins.txt
done
This process then re-names bins with a prefix ROUGH, which we then
used to identify bins to go back and re-refine. After this second pass,
we then modified the above script to pull out all of the refined MAGs
that met better standards (less than 10% redundant and either ≥ 80%
complete or ≥ 2 million nucleotides in length, as many oral microbes
have approximately 2 Mbp genomes).
#!/bin/bash
for SET in `ls | grep "A[A-Z]-idba_ud"`
do
indiv=$(echo "$SET" | sed 's/-.*$//')
# define the bins
anvi-rename-bins -c $SET/$indiv-CONTIGS.db \
-p $SET/$indiv-MERGED/PROFILE.db \
--collection-to-read ROUGH \
--collection-to-write FINAL \
--call-MAGs \
--size-for-MAG 2 \
--min-completion-for-MAG 80 \
--max-redundancy-for-MAG 10 \
--prefix FINAL_$indiv \
--report-file $SET/$indiv-final-bins.txt
# summarize each final collection
anvi-summarize -c $SET/$indiv-CONTIGS.db \
-p $SET/$indiv-MERGED/PROFILE.db \
-C FINAL \
-o $SET/$indiv-SUMMARY
done
Merging co-assemblies
From co-assembly, the possibility exists that the same microbial populations could be binned independently in multiple mouths. However, this issue raises some secondary issues - what exactly is same in this case? 1 SNV? 100 SNVs? Some percent of genes shared/not shared? Average nucleotide identity? Similar abundance trajectories? These definition-related issues aren’t identical, nor are they hierarchical, and they all make different assumptions about microbial ecology and evolution. So we took an approach to identify ‘redundant’ MAGs that accomodates all of these important perspectives on microbial populations, which are applied in these next few steps. These steps are very much inspired by the anvi’o documentation and the Meren lab, with some minor tweaks for oral microbes.
Identifying and removing redundant MAGs
First we collected all of the MAGs and placed them into a directory:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 1 # cores from each
#SBATCH --contiguous
#SBATCH --mem=8G #per node
#SBATCH -t 0-08:00:00
#SBATCH -p shared
#SBATCH --job-name="mergeRedundant"
#SBATCH -o odyssey_mergeRedundant.out
#SBATCH -e odyssey_mergeRedundant.err
#SBATCH --mail-type=END
mkdir REDUNDANT-MAGs
mkdir REDUNDANT-MAGs/indiv_MAGs
source ~/virtual-envs/anvio-dev-venv/bin/activate
for SET in `ls | grep "A[A-Z]-idba_ud"`
do
indiv=$(echo "$SET" | sed 's/-.*$//')
MAGs=`anvi-script-get-collection-info -c $SET/$indiv-CONTIGS.db -p $SET/$indiv-MERGED/PROFILE.db -C
FINAL | grep MAG | awk '{print $1}'`
for MAG in $MAGs; do
# reformat fasta and put into redundant MAGs
anvi-script-reformat-fasta $SET/$indiv-SUMMARY/bin_by_bin/$MAG/$MAG-contigs.fa --simplify-names
--prefix $MAG -o REDUNDANT-MAGs/indiv_MAGs/$MAG.fa
done
done
cat REDUNDANT-MAGs/indiv_MAGs/*fa > REDUNDANT-MAGs/indiv_MAGs/REDUNDANT-MAGs.fa
NB: from this step forward, the subsequent analyses are applied singly to all MAGs generated from any mouth; no longer are processes done by volunteer.
And then we ran the steps above (generating CONTIGS and merged
PROFILE databases). Note that at this point, since these are merged,
we mapped ALL metagenomes, i.e., 81 metagenomes from 17 participants,
onto the combined CONTIGS database of 409 MAGs from all co-assemblies.
This metagenomic read recruitment is thus competitive, with each read
being mapped against all contigs from all MAGs before being assigned to
to the contig locus to which it maps best.
Correlation
MAGs were correlated to determine which MAGs were tightly correlated in coverage (indicative of being potentially duplicate MAGs).
#!/usr/bin/env Rscript
library(reshape2)
options(max.print=1000000000)
mean_coverage <- t(read.delim(file="REDUNDANT-MAGs-SUMMARY/bins_across_samples/mean_coverage.txt",
header = TRUE,
stringsAsFactors = FALSE,check.names = F, row.names = 1))
correlations <- melt(data = cor(mean_coverage, use="complete.obs", method="pearson"))
names(correlations) <- c('MAG_1', 'MAG_2', 'correlation')
#write.table(correlations, "REDUNDANT-MAGs-PEARSON.txt", sep="\t", quote=FALSE, col.names=NA)
Taxonomy
We also ran GTDB-Tk for taxonomy:
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 8 # 12 cores from each
#SBATCH --contiguous
#SBATCH --mem=180G #per node
#SBATCH -t 0-06:00:00
#SBATCH -p shared
#SBATCH --job-name="gtdbtk"
#SBATCH -o odyssey_gtdbtk.out
#SBATCH -e odyssey_gtdbtk.err
#SBATCH --mail-type=END
source ~/virtual-envs/anvio-dev-github/bin/activate
export GTDBTK_DATA_PATH=$PWD/release95
gtdbtk classify --align_dir gtdbtk/align --genome_dir gtdbtk/genomes --out_dir gtdbtk/classify --extension fa --cpus 8
ANI
Taxonomy was then used when calculating average nucleotide identity (ANI) between read pairs so that ANI was only compared intra-phyla for efficiency and accuracy:
#!/bin/bash
grep "MAG" REDUNDANT-MAGs-gtdbtk.txt | awk '{print $1 "\t" $4}' | sed 's/k__.*;c__//; s/;.*$//; s/k__//' > REDUNDANT-MAGs-AFFILIATIONS.txt
awk '{print $1"\t"$3}' REDUNDANT-MAGs-SUMMARY/bins_summary.txt > REDUNDANT-MAGs-LENGTH.txt
B="REDUNDANT-MAGs-SUMMARY/bin_by_bin"
echo -e "MAG\tlength\tcompletion\tredundancy" > REDUNDANT-MAGs-STATS.txt
for MAG in $(ls --color=none $B | grep MAG); do
length=`cat $B/$MAG/$MAG-total_length.txt`
completion=`cat $B/$MAG/$MAG-percent_completion.txt`
redundancy=`cat $B/$MAG/$MAG-percent_redundancy.txt`
echo -e "$MAG\t$length\t$completion\t$redundancy"
done >> REDUNDANT-MAGs-STATS.txt
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 32 # 12 cores from each
#SBATCH --contiguous
#SBATCH --mem=120G #per node
#SBATCH -t 1-12:00:00
#SBATCH -p shared
#SBATCH --job-name="ANI"
#SBATCH -o odyssey_anviANI.out
#SBATCH -e odyssey_anviANI.err
#SBATCH --mail-type=END
# need to activate venv directly because node receiving job order doesn't play with bashrc aliases
source ~/virtual-envs/anvio-dev-venv/bin/activate
module load MUMmer
anvi-compute-ani -i REDUNDANT-MAGs-internal-genomes-table.txt --method ANIm -o REDUNDANT-MAGs-ANI -T 32
Defining redundant MAGs
With this collected information (coverage correlation, ANI, taxonomy), we proceeded to define and remove redundant MAGs (retaining the higher-quality MAG). To tweak this, we used Alon Shaiber’s redundancy script, modified slightly:
#!/usr/bin/env python
import argparse
import pandas as pd
__author__ = "Alon Shaiber"
__copyright__ = "Copyright 2017, The anvio Project"
__status__ = "Development"
parser = argparse.ArgumentParser(description='A short helper script to generate a table that is ready for redundancy analysis (the next step is export-non-redundant-genomes.py)')
group0 = parser.add_argument_group('ESSENTIAL INPUTS', "You must supply the following TAB-delimited files.")
group1 = parser.add_argument_group('ESSENTIAL OUTPUT', "You must supply a name for the output file.")
group0.add_argument('-m', '--metadata', metavar='PATH', default=None, dest='metadata', help='A TAB-delimited file with metadata regarding each genome. It has to contain the following columns: name, length, completion, redundancy, domain.')
group0.add_argument('-C', '--correlation', metavar='PATH', default=None, dest='correlation', help='A TAB-delimited file containing the correlation information between the occurance of pairs of genomes. It must contain the following columns: key, item_1, item_2, correlation.')
group0.add_argument('-a', '--ani', metavar='PATH', default=None, dest='ani', help='A TAB-delimited matrix with alignment percentages pairs of genomes (for example the output from anvi-compute-ani)')
group0.add_argument('-l', '--ani_length', metavar='PATH', default=None, dest='ani_len', help='A TAB-demlimited matrix of alignment lengths, e.g. reported by anvi-compute-ani')
group1.add_argument('-t', '--taxonomy', metavar='PATH', default=None, dest='taxonomy', help='A TAB-delimited file containing the taxonomic affiliation for each MAG')
group1.add_argument('-o','--output', metavar='PATH', default='FINAL-MAGs-LIST.txt', help='Name of output file for the data table')
group1.add_argument('-A','--minANI', metavar='DECIMAL', default=0.98, help='Minimum percent ANI to be considered redundant [DEFAULT 0.98]')
group1.add_argument('-L','--minAlign', metavar='DECIMAL', default=0.75, help='Minimum fraction of the smaller MAG\'s length for considering the alignment threshold of --minANI [DEFAULT 0.75]')
group1.add_argument('-R','--minCor', metavar='DECIMAL', default=0.85, help='Minimum correlation coefficient to be considered redundant [DEFAULT 0.85]')
args = parser.parse_args()
minANI = args.minANI
minAlign = args.minAlign
minCor = args.minCor
out_filepath = args.output
metadata = pd.read_csv(args.metadata, sep='\t')
if args.taxonomy:
taxonomy = pd.read_csv(args.taxonomy, sep='\t', header=None)
taxonomy.columns = ['MAG', 'tax']
correlation = pd.read_csv(args.correlation, sep='\t')
correlation.columns.values[0] = 'key'
ani = pd.read_csv(args.ani, sep='\t')
ani = pd.melt(ani, id_vars=['key'])
ani.columns = ['MAG_1', 'MAG_2', 'ANI']
ani_len = pd.read_csv(args.ani_len, sep='\t')
redundant_MAGs = metadata['MAG'].unique()
redundant_groups = []
for MAG in redundant_MAGs:
# find MAGs > minANI
high_ANI = list(ani.loc[(ani['MAG_1'] == MAG) # for this mag
& (ani['ANI'] >= minANI) # if ANI is high enough
& (ani['MAG_2'] != MAG), # but not this mag itself
'MAG_2']) # report the other MAG
long_ANI = [x for x in high_ANI if # if alignment fraction is >= minAlign of the smaller genome
float(ani_len.loc[ani_len['key'] == MAG, x] /
metadata.loc[(metadata['MAG'] == x) | (metadata['MAG'] == MAG), 'length'].min()) >= minAlign]
correlated = [x for x in long_ANI if
float(correlation.loc[(correlation['MAG_1'] == MAG) & (correlation['MAG_2'] == x), 'correlation']) >= minCor]
if correlated:
correlated.append(MAG)
redundant_groups.append(correlated)
MAGs_to_keep = redundant_MAGs
for grp in redundant_groups:
best_MAG = metadata.loc[(metadata['MAG'].isin(grp)) &
(metadata['completion'] == metadata.loc[metadata['MAG'].isin(grp),'completion'].max()),'MAG'].item()
mags_to_replace = [x for x in grp if x is not best_MAG]
MAGs_to_keep = [best_MAG if x in mags_to_replace else x for x in MAGs_to_keep]
print(str(len(redundant_MAGs)) + " MAGs came in")
print(str(len(set(MAGs_to_keep))) + " MAGs made it out")
print("Saving output to: " + out_filepath)
with open(out_filepath, 'w') as f:
[f.write(x + '\n') for x in sorted(set(MAGs_to_keep))]
We applied this script to identify a list of non-redundant MAGs with the following command:
./anvi-define-redundant-MAG-groups.py \
-m REDUNDANT-MAGs-STATS.txt \ # metadata of MAGs so the better MAG can be retained
-C REDUNDANT-MAGs-PEARSON.txt \ # correlations
-a REDUNDANT-MAGs-ANI/ANIm_percentage_identity.txt \ # ANI percent identity
-l REDUNDANT-MAGs-ANI/ANIm_alignment_lengths.txt \ # ANI alignment lengths
-o NON-REDUNDANT-MAGs-LIST.txt
Dropping redundant MAGs
Having identified the non-redundant set of MAGs, we concatenated their
corresponding fastas into a single fasta with
cat indiv_MAGs/*fa > $indiv.fa, which then was used to create a new
anvi’o contigs database NON-REDUNDANT-MAGs-CONIGS.db using the same
script as before except for setting the indiv variable to
NON-REDUNDANT-MAGs. This merged fasta and corresponding contigs db was
then used with the exact workflow as before to recruit reads from all
samples, profile the results, and merge the profiles into
NON-REDUNDANT-MAGs-MERGED. The association between contigs and MAGs
was then reconstructed from the contig IDs with the following lines:
grep '>' NON-REDUNDANT-MAGs.fa | awk '{print $0 "\t" $0}' | sed 's/>//g; s/_[0-9]*$//' > NON-REDUNDANT-MAGs-COLLEC
TION.txt
anvi-import-collection NON-REDUNDANT-MAGs-COLLECTION.txt -c NON-REDUNDANT-MAGs-CONTIGS.db -p NON-REDUNDANT-MAGs-MER
GED/PROFILE.db -C NON_REDUNDANT_MAGs --contigs-mode
anvi-summarize -c NON-REDUNDANT-MAGs-CONTIGS.db -p NON-REDUNDANT-MAGs-MERGED/PROFILE.db -C NON_REDUNDANT_MAGs -o NO
N-REDUNDANT-MAGs-SUMMARY
Finishing touches
With this final set of MAGs, we finalized the annotation and taxonomy.
Taxonomic assignment (GTDB)
We used GTDB-Tk by Parks et al., 2015
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 32 # 12 cores from each
#SBATCH --contiguous
#SBATCH --mem=64G #per node
#SBATCH -t 1-00:00:00
#SBATCH -p shared
#SBATCH --job-name="gtdbtk-MAGs"
#SBATCH -o odyssey_gtdbtk.out
#SBATCH -e odyssey_gtdbtk.err
#SBATCH --mail-type=END
# need to activate venv directly because node receiving job order doesn't play with bashrc aliases
source ~/virtual-envs/gtdbtk/bin/activate
prefix="REDUNDANT-MAGs"
#cd indiv_MAGs
gtdbtk tree indiv_MAGs -x .fa -t 32 `pwd`/$prefix-gtdbtk-TREE
gtdbtk tree_qa `pwd`/$prefix-gtdbtk-TREE -f $prefix-gtdbtk.txt
Adding EggNOG annotations
EggNOG annotations were run as above, also InterProScan was run as described in http://dutter.github.io/projects/oral_metapan#functional-annotation
DESMAN to identify subpopulations
For these final non-redundant MAGs, we used
DESMAN to identify
‘haplotypes,’ or putatively co-inherited sets of genes inherited
together. Briefly, the method works by identifying SNVs for informative
positions from bacterial core genes (which we have from above, during
the anvi-profile and anvi-run-hmms steps), and using the frequencies
of each nucleotide variant to estimate which occur on the same
chromosome. We use the results of this approach as an operational
definition for sub-population: temporally distinct genome types within
(=sub-) a set of closely-related genomes represented by a high quality
MAG (= population).
Exporting SNVs from anvi’o
We used the set of 71 bacterial core genes identified by anvi’o. This
script loops through all MAGs in the NON-REDUNDANT set, and identifies
each mag’s gene callers for those HMM-identified core genes, and
thenanvi-gen-variability-profile outputs the SNV frequencies for
variant sites in each gene. Note that the
--min-coverage-in-each-sample flag was set to 0, as DESMAN has its
own method for filtering out informative vs. uninformative SNV sites, so
we pass it all possible variant frequencies for it to do that.
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 1 # 12 cores from each
#SBATCH --contiguous
#SBATCH --mem=12G #per node
#SBATCH -t 0-3:00:00
#SBATCH -p shared
#SBATCH --array=0-29%30
#SBATCH --job-name="varDes"
#SBATCH -o odyssey_var_for_des.out
#SBATCH -e odyssey_var_for_des.err
#SBATCH --mail-type=NONE
# set this up first
# split -l 13 -d -a 4 NON-REDUNDANT-MAGs-LIST.txt NON-REDUNDANT-MAGs-LIST.txt-
# mkdir mag_bacteria71_genes
prefix=NON-REDUNDANT-MAGs
taskNum=$(printf %04d $SLURM_ARRAY_TASK_ID) # format slurm task id to be padded with 0s like split
for MAG in $(cat NON-REDUNDANT-MAGs-LIST.txt-$taskNum); do
anvi-get-sequences-for-hmm-hits -c $prefix-CONTIGS.db -p $prefix-MERGED/PROFILE.db -C NON_REDUNDANT_MAGs -b $MAG --hmm-sources Bacteria_71 -o mag_bacteria71_genes/$MAG-genes.fna
gene_calls=$(grep '>' mag_bacteria71_genes/$MAG-genes.fna | sed 's/^.*gene_callers_id:\([0-9]*\).*$/\1/' | tr '\n' ',' | sed 's/,$//')
# --quince-mode essential for DESMAN later
anvi-gen-variability-profile -c $prefix-CONTIGS.db -p $prefix-MERGED/PROFILE.db --gene-caller-ids $gene_calls --engine NT --min-coverage-in-each-sample 0 --quince-mode -o mag_bacteria71_genes/$MAG-Bacteria_71-SNVs.txt
done
However, DESMAN requires different input formatting than anvi’o outputs (see DESMAN github), which was generated with the following short R code:
MAGs <- list.files("catalyst_idba/mag_bacteria71_genes/", pattern=".txt", full.names = F)
for (MAG in gsub("-Bacteria.*$", "", MAGs)){
#MAG <- "FINAL_AI_MAG_00002"
vars <- read.csv(paste0("catalyst_idba/mag_bacteria71_genes/", MAG, "-Bacteria_71-SNVs.txt"), sep="\t", header=T)
vars <- melt(vars, id.vars = c('corresponding_gene_call','sample_id','pos'), measure.vars = c("A","T","C","G"))
vars$sampleNuc <- paste(vars$sample_id, vars$variable, sep = "-")
vars <- dcast(vars, corresponding_gene_call + pos ~ sampleNuc, value.var = 'value', fill=0)
colnames(vars)[1:2] <- c('Contig', 'Position')
write.table(vars, paste0("catalyst_idba/mag_bacteria71_genes/", MAG, "-Bacteria_71-SNVs.csv"), quote = F, col.names = T, row.names = F, sep = ",")
}
DESMAN variant sites
This formatted input is then passed to DESMAN for the tool to apply its variant filtration method and identify informative positions:
# inside folder DESMAN
for stub in $(ls catalyst_idba/mag_bacteria71_genes/*csv | sed 's/.csv//'); do
desman/Variant_Filter.py $stub.csv -o $stub-DESMAN -p
done
The selected variants were then used for the next step.
DESMAN
DESMAN requires an estimation of the expected number of haplotypes/strains to exist within the population. This is clearly difficult to know for certain prior, so the recommended approach is to run desman for several different G’s (number of strains to deconvolve into) to generate decompositions for all reasonable/feasible numbers. From these, the posterior deviances can then be used to estimate the tradeoff point between better fit to the observed frequencies vs. overfitting.
We ran these on the cluster like so (note that the ls ... split line
was run first to set up numerically incremented batches for the slurm
array):
#!/bin/bash
#SBATCH -N 1 #1 nodes of ram
#SBATCH -n 10 # 10 cores from each
#SBATCH --contiguous
#SBATCH --mem=124G #per node
#SBATCH -t 0-32:00:00
#SBATCH -p shared
#SBATCH --array=0-363%364
#SBATCH --job-name="desman"
#SBATCH -o odyssey_desman_%a.out
#SBATCH -e odyssey_desman_%a.err
#SBATCH --mail-type=NONE
# set this up first
# ls mag_bacteria71_genes/*SNVs.csv | sed 's/.csv//' | split -l 1 -d -a 4 - mag_bacteria71_genes/desBatches/batch-
desmanDir=mag_bacteria71_genes
taskNum=$(printf %04d $SLURM_ARRAY_TASK_ID) # format slurm task id to be padded with 0s like split
MAG=$(cat $desmanDir/desBatches/batch-$taskNum)
source ~/virtual-envs/anvio-dev-github/bin/activate
for strains in 1 2 3 4 5 6 7 8 10 12; do
echo $strains # so we know what ID we're on
(desman $MAG-DESMANsel_var.csv -g $strains -i 500 -e $MAG-DESMANtran_df.csv -o $MAG-DESMAN-FINAL_${strains}_0)&
# the ()& means run and detatch as subprocesses since single threads, then wait for all to finish
done; wait
Note that only 364 MAGs entered the deconvolution step; 26 MAGs have fallen out by this point due to a lack of usable SNVs.
In reality, this DESMAN step in will crash for many due to memory or
time requirements, since this is attempting to run 364 jobs that each
run 10 tasks in parallel of varying complexity, so this script ended up
being run iteratively with varying #SBATCH --array IDs. That is, we
first ran it on all 364x10, caught the array ID’s for jobs that timed
out or crashed from memory, updated the memory/time requests and set the
--array to use only those failed runs, and repeated this until all
were done.
Identifying optimum DESMAN strains
We then concatenated the fit results from all strain numbers per MAG like so:
for stub in $(ls *-SNVs.csv | sed 's/.csv//'); do cat $stub*/fit.txt | cut -d"," -f2- > $stub-DESMAN-DEV.csv; done
Then heuristic to estimate optimal number of haplotypes to use: -f is
minimum fraction (at least 1%), -d is %step in posterior deviance,
-m is max errors
for stub in $(ls *-SNVs.csv | sed 's/.csv//'); do resolvenhap.py -d 0.4 -m 2 -f 0.01 $stub-DESMAN-FINAL; done >> desmanStrainOptima.csv
Then in R, we plotted the posterior deviance for all G’s while highlighting the predicted optimum G to make sure the heuristic produced the expected, or at least reasonable, results:
library(ggplot2)
MAG <- "FINAL_AI_MAG_00008"
desStrains <- read.csv('catalyst_idba/mag_bacteria71_genes/desmanStrainOptima.csv', col.names = c('nHap','nGoodHap','index','avg_err','base_pred_file'))
numStrains <- desStrains$nHap[grep(MAG,desStrains$base_pred_file)]
desmanStats <- read.csv(paste0("catalyst_idba/mag_bacteria71_genes/", MAG, "-Bacteria_71-SNVs-DESMAN-DEV.csv"), header=F, col.names = c("H","G","LP","Dev"))
desmanStats$Gfactor <- factor(desmanStats$G, sort(desmanStats$G))
desColors <- setNames(rep('black', length(desmanStats$G)),desmanStats$G)
desColors[as.character(numStrains)] <- 'red'
ggplot(desmanStats, aes(x = G, y = Dev)) + geom_smooth()+ geom_point(aes(color=Gfactor)) + theme_bw() +
ylab("Mean posterior deviance") + xlab("Number of strains - G") +
scale_color_manual(values = desColors, guide=F) +
theme(axis.text=element_text(size=12),axis.title=element_text(size=14,face="bold"), panel.grid = element_blank(), panel.border = element_blank(),
axis.line = element_line())
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
And for the example here, FINAL_AI_MAG_00008, the chosen number of
haplotypes (4) is indeed around the knee in the deviance curve, so 4
haplotypes does appear optimal.
Analyses and plots
The majority of the remaining analyses and figures were generated with R, documented with the code below
Libraries and basic constants
library(ggplot2)
library(ggridges)
library(reshape2)
library(plyr)
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.6.2
library(gtools)
library(grid)
library(scales)
library(vegan)
library(beyonce)
library(tidyr)
library(RColorBrewer)
theme_basic <- theme_classic() + theme(axis.ticks.x = element_blank(), axis.text = element_text(color = 'black'), plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5), strip.background = element_rect(fill = '#dcdcdc', colour = NA))
theme_basic_horizontal <- theme_basic + theme(axis.text.x = element_text(color = 'black', angle = 90, hjust = 1, vjust = 0.5))
dans_colors <- c("#FF0000","#F6A300","#0068CC","#6600AA","#AC0088","#AA33FF","#00FFFF","#00CC00","#006611","#00AC99",
"#AC6844","#FFFF00",
"#991100","#ACAC11","#a0f0aa","#FF00FF","#FF8611","#B9F6F6","#001166","#AC9A00","#994141","#ff1169",
"#0AF622","#119924","#Ac3311","#004A9A","#AcAc99","turquoise","tomato","sienna1","rosybrown",
"peachpuff","olivedrab3","mistyrose1","mediumorchid","indianred2","#114914","#660011","ivory3","deeppink","#331111")
cogValues <- c(
"K"="K: Transcription",
"L"="L: Replic., Recom., Repair",
"D"="D: Cell division/chrom. partition",
"E"="E: AA trans./metabol.",
"J"="J: Translation, ribo. struct.",
"G"="G: Carbohydrate trans./metabol.",
"C"="C: Energy production/metabol.",
"V"="V: Defense",
"P"="P: Inorg. ion trans./metabol.",
"H"="H: Coenzyme trans./metabol.",
"F"="F: Nucl. trans/metabol.",
"U"="U: Intracell. traf., sec., trans.",
"T"="T: Signal transduction",
"Q"="Q: Sec. metabolites",
"O"="O: Post-transl. modification",
"I"="I: Lipid trans./metabol.",
"M"="M: Cell envelope/wall",
"N"="N: Cell motility",
"S"="S: Function Unknown"
)
Functions used
makegtdbtktaxMat <- function(gtdbtkPath="catalyst_idba/gtdbtk/classify/gtdbtk.bac120.summary.tsv", shortform=T){
gtdbtk <- read.csv(gtdbtkPath, header = T, sep = "\t")
colnames(gtdbtk)[1] <- c("MAG")
gtdbtk <- gtdbtk[-nrow(gtdbtk),] # drop row of dashes
gtdbtk$participant <- gsub("^.*_(A[A-Z])_.*$","\\1", as.character(gtdbtk$MAG))
# clean up c/g/etc prefixes in naming
gtdbtk$classification <- gsub("[a-z]__","",as.character(gtdbtk$classification))
taxMat <- rbind.fill.matrix(lapply(strsplit(as.character(gtdbtk$classification), ";"), function(x) matrix(x, ncol = length(x), byrow=T)))
taxMat[is.na(taxMat)] <- 'Unidentified'
colnames(taxMat) <- c("Domain","Phylum","Class","Order","Family","Genus","Species")
gtdbtk <- cbind(gtdbtk,as.data.frame(taxMat))
if(shortform){
gtdbtk <- gtdbtk[,c("MAG","participant",c("Domain","Phylum","Class","Order","Family","Genus","Species"),"classification_method","note","warnings")]
}
return(gtdbtk)
}
correlateByBin <- function(infoDF, covDF, binDF, bin, minMedCov=1, maxCorThresh = 1, minCorThresh=-1, minCountThresh=1) {
covDF <- covDF[covDF$key %in% infoDF[infoDF$bin == bin,"gene_callers_id"],]
indivs <- unique(gsub("[0-9]*$","",colnames(covDF)[2:ncol(covDF)]))
corDF <- data.frame(gene=rep(covDF$key, length(indivs)), Individual=sort(rep(indivs,nrow(covDF))), bin=bin, p.value=NA, correlation=NA)
for (indiv in indivs) {
binTrace <- as.numeric(binDF[binDF$bins == bin, grep(indiv, colnames(binDF))])
if (median(binTrace) >= minMedCov){
for (gene in unique(as.character(covDF$key))){
geneTrace <- as.numeric(covDF[covDF$key == gene, grep(indiv, colnames(covDF))])
geneCor <- cor.test(geneTrace, binTrace, method = 'spearman', exact = F)
corDF[(corDF$gene == gene) & (corDF$Individual == indiv), c("p.value","correlation")] <- c(geneCor$p.value, geneCor$estimate)
}
} else {
print(paste0(indiv, " had <", minMedCov, "X coverage (median), so skipping. Bye!"))
corDF[corDF$Individual == indiv, "bin"] <- NA
}
}
corDF <- corDF[!is.na(corDF$bin),]
if (nrow(corDF) > 0){
if (minCountThresh == 'all') {
minCountThresh <- length(unique(as.character(corDF$Individual)))
}
corDFrecast <- dcast(corDF, gene ~ Individual, value.var="correlation")
corDF <- corDF[!is.na(corDF$correlation),]
return(corDF[(corDF$correlation > minCorThresh) & (corDF$correlation < maxCorThresh),])
} else {
print("NOTHING PASSES FILTERS! Maybe reconsider filters?")
return()
}
}
traceGenesByMAG <- function(infoDF, covDF, binDF, bin, statsDF, focusGenes = NULL, minGeneCov=0, tax=makegtdbtktaxMat(), facetGrid=T, focusColor=NA, focusIn=NA,
sqrt_transform=F, log10_transform=T) {
covDF <- covDF[covDF$key %in% as.character(infoDF[infoDF$bin == bin,"gene_callers_id"]),] # get coverages for this MAG
print(paste0(nrow(covDF), " genes in ", bin))
if ((!is.null(focusGenes)) & is.na(focusColor)){
covDF <- covDF[covDF$key %in% focusGenes,]
}
# get overall coverage for MAGs
binTrace <- setNames(as.numeric(binDF[binDF$bins == bin, grep("A", colnames(binDF))]), grep("A", colnames(binDF), value = T))
binTrace <- binTrace[mixedorder(names(binTrace))] # sort human (e.g. AB1 comes before AB15 before AB2)
# make long for ggplotting
covDF <- melt(covDF, id.vars = 'key')
colnames(covDF) <- c('Gene','Sample','Coverage')
# order things nicely and make participant
covDF$Gene <- factor(covDF$Gene)
covDF$Day <- gsub("A[A-Z]","", as.character(covDF$Sample))
covDF$Day <- factor(covDF$Day, mixedsort(unique(as.character(covDF$Day)))) # sort human
covDF$Participant <- gsub("(A[A-Z]).*$","\\1",as.character(covDF$Sample))
# grab genes above threshold
if (minGeneCov > 0) {
geneThresh <- ddply(covDF, ~ Gene + Participant, summarise, meanCov = mean(Coverage)) %>% filter(meanCov >= minGeneCov) %>% pull(Gene)
covDF <- covDF[covDF$Gene %in% unique(geneThresh),]
}
# make a little dataframe to facilitate overlaying the bin's coverage
binCov <- data.frame(Sample=factor(gsub("A[A-Z]","",names(binTrace)), levels = unique(gsub("A[A-Z]","",names(binTrace)))), Coverage=binTrace,
Participant = gsub("(A[A-Z]).*$","\\1",names(binTrace)))
showLegend <- F
if (length(unique(as.character(covDF$Gene))) <= 20) {
showLegend <- "legend"
print(setNames(infoDF$Pfam[match(unique(as.character(covDF$Gene)), infoDF$gene_callers_id)], unique(as.character(covDF$Gene))))
} else {
print("More than 20 genes are here, skipping legend and not dumping functions")
}
p <- ggplot(covDF, aes(x = Day, y = Coverage))
if (is.na(focusColor)) { ## case 1: no genes to focus
p <- p + geom_line(aes(group = Gene), color='grey70', size = 0.1, alpha=0.6)
} else if (is.na(focusIn)) { ## case 2: genes to focus but globally, not participant specific
geneColors <- setNames(rep('grey70', length(levels(covDF$Gene))), levels(covDF$Gene))
geneColors[as.character(focusGenes)] <- focusColor
p <- p + geom_line(aes(group = Gene, color=Gene), size = 0.1, alpha=0.6) +
scale_color_manual(values = geneColors)
} else { ## case 3: genes to focus but in specific participants
covDF$combinedGene <- paste0(as.character(covDF$Participant),"x",as.character(covDF$Gene))
geneColors <- setNames(rep('grey70', nrow(covDF)), as.character(covDF$combinedGene))
geneColors[paste0(as.character(focusIn),"x",as.character(focusGenes))] <- focusColor
p <- p + geom_line(data=covDF, aes(group = combinedGene, color=combinedGene), size = 0.1, alpha=0.6) + # reload data since updated here
scale_color_manual(values = geneColors)
}
p <- p +
geom_line(data = binCov, aes(group = Participant, x=Sample), size = 1, color = 'black') +
scale_x_discrete(expand = c(0,0)) + labs(x="Sample Day") +
labs(title=bin, subtitle = paste0(paste0(as.character(as.matrix(tax[tax$MAG == bin, c("Genus","Species")])), collapse = " "), # this is "Genus sp."
": ",round(statsDF$completion[statsDF$MAG == bin],1), "% C / ", round(statsDF$redundancy[statsDF$MAG == bin], 1), "% R, ",
round(statsDF$length[statsDF$MAG == bin]/1e6, 2), " Mb")) + guides(color=showLegend) +
theme_basic + theme(strip.background = element_blank(), panel.spacing = unit(0.1, 'lines'), axis.ticks.x = element_line(), axis.text.x = element_text(size = 4))
if (facetGrid){
p <- p + facet_grid(~Participant, space = 'free_x', scales = 'free_x')
} else {
p <- p + facet_wrap(~Participant, ncol = ceiling(length(unique(covDF$Participant))/3), scales = 'free_x')
}
if (log10_transform){
p <- p + scale_y_log10(expand = c(0,0))#, trans = 'log10', breaks = trans_breaks(trans = 'log10', inv = function(x) 10^(x + 1), n = 8))
} else if (sqrt_transform){
p <- p + scale_y_continuous(expand = c(0,0), trans = 'sqrt', breaks = trans_breaks(trans = 'sqrt', inv = function(x) x^2, n = 8))
} else {
p <- p + scale_y_continuous(expand = c(0,0))
}
print(p)
}
detrendGenesByMAG <- function(infoDF, covDF, binDF, bin, statsDF, focusGenes = NULL) {
covDF <- covDF[covDF$key %in% as.character(infoDF[infoDF$bin == bin,"gene_callers_id"]),] # get coverages for this MAG
print(paste0(nrow(covDF), " genes in ", bin))
if (!is.null(focusGenes)){
covDF <- covDF[covDF$key %in% focusGenes,]
}
binTrace <- binDF[binDF$bins == bin, grep("A", colnames(binDF))]
covDF[, grep("A", colnames(covDF))] <- covDF[, grep("A", colnames(covDF))] - matrix(binTrace) # subtract bin cov vector from each gene vector
# make long for ggplotting
covDF <- melt(covDF, id.vars = 'key')
colnames(covDF) <- c('Gene','Sample','Coverage')
# order things nicely and make participant
covDF$Gene <- factor(covDF$Gene)
covDF$Day <- gsub("A[A-Z]","", as.character(covDF$Sample))
covDF$Day <- factor(covDF$Day, mixedsort(unique(as.character(covDF$Day)))) # sort human
covDF$Participant <- gsub("(A[A-Z]).*$","\\1",as.character(covDF$Sample))
return(covDF)
}
detrendDistribution <- function(infoDF, covDF, binDF, bin, statsDF, focusGenes = NULL, sdMultiple=2, tax=makegtdbtktaxMat(), plot_it=T, focusColor='#Fa6a3a') {
covDF <- detrendGenesByMAG(infoDF = infoDF, covDF=covDF, binDF=binDF, bin=bin, focusGenes = focusGenes)
covDF <- covDF %>% group_by(Participant) %>% mutate(outlierInParticipant=Coverage > sdMultiple*sd(Coverage), outlierThresh=sdMultiple*sd(Coverage))
if (plot_it){
p <- ggplot(covDF, aes(Coverage, y=Participant)) +
#geom_histogram(aes(y=..density..), bins = 20, fill='grey50') +
#geom_density_ridges() +
stat_density_ridges(bandwidth = 4) +
#geom_density_ridges(stat='binline', bins=100) +
geom_segment(aes(y=as.numeric(factor(Participant)), yend=as.numeric(as.factor(Participant)) + 1, x=outlierThresh, xend=outlierThresh), color=focusColor, linetype=2, size=0.5) +
scale_y_discrete(expand = c(0,0)) + scale_x_continuous( n.breaks = 3) + labs(x="Coverage deviation") +
labs(title=bin, subtitle = paste0(paste0(as.character(as.matrix(tax[tax$MAG == bin, c("Genus","Species")])), collapse = " "), # this is "Genus sp."
": ",round(statsDF$completion[statsDF$MAG == bin],1), "% C / ", round(statsDF$redundancy[statsDF$MAG == bin], 1), "% R, ",
round(statsDF$length[statsDF$MAG == bin]/1e6, 2), " Mb")) + guides(color=F) +
theme_basic + theme(strip.background = element_blank(), panel.spacing = unit(0.1, 'lines'), axis.ticks.x = element_line())
print(p)
}
covDF <- covDF %>% group_by()
return(covDF)
}
plotCorrWithMAG <- function(corDF=correlateByBin(), sigThresh = 1, absCorThresh = 1, tax=makegtdbtktaxMat()){
corDF <- corDF[corDF$p.value <= sigThresh,]
corDF$gene <- factor(corDF$gene, levels = sort(unique(as.character(corDF$gene))))
p <- ggplot(corDF, aes(x = gene, y = Individual)) +
geom_tile(aes(fill = correlation)) +
scale_y_discrete(expand = c(0,0)) +
scale_fill_gradient2(low = 'blue', mid = 'lightgoldenrodyellow', high = 'red', midpoint = 0, limits = c(-1,1), na.value = 'grey10') +
theme_basic
print(p)
return(corDF$gene[(corDF$correlation <= absCorThresh) & (!is.na(corDF$gene))])
}
plotFuncCounts <- function(vec=matched$`infoDF$Pfam[match(corDF$gene, infoDF$gene_callers_id)]`, minCount = 100, logAx='identity', title="", dropNA=F) {
if (dropNA) {
vec <- vec[vec != ""]
}
vec <- factor(vec, levels = levels(vec)[rev(order(table(vec)))])
qplot(vec[vec %in% names(table(vec)[table(vec) >= minCount])]) + scale_y_continuous(expand = c(0,0), trans=logAx) +
labs(x="Function", y="Count", title = title) +
theme_basic_horizontal
}
plotDiversityTime <- function(matForDiv, title="MAG", method='simpson') {
vec <- diversity(matForDiv, index = method, MARGIN = 2)
df <- data.frame(div=vec, samp=names(vec), indiv=gsub("[0-9]*","",names(vec)))
df$samp <- gsub("-","", df$samp)
df$day <- gsub("A[A-Z]", "", as.character(df$samp))
df$day <- factor(df$day, levels = mixedsort(unique(as.character(df$day))))
p <- ggplot(df, aes(y=div, x=day, group=indiv)) +
labs(x='Sample day', #title=paste0('Diversity: ', paste0(title, collapse = " ")),
y = paste0(toupper(unlist(strsplit(method,""))[1]), paste0(unlist(strsplit(method,""))[-1], collapse = ""), " diversity")) +
geom_line(stat = 'identity', color='black') + geom_point() +
facet_grid(~indiv, scales = 'free_x', space = 'free_x') +
scale_y_continuous(expand = c(0,.01)) + theme_basic_horizontal +
theme(panel.spacing.x = unit(0.1, 'lines'), strip.background = element_blank(), axis.line.x.bottom = element_line(), axis.text.x = element_text(hjust=1))
print(p)
return(cbind(df,dataset=title))
}
correlateMagGenePairs <- function(geneCovs, geneInfo, binDetect, byIndiv=T, detectThresh=0.75, excludeOwn=F, corMeth='pearson',
MAGs=unlist(sapply(levels(binDetect$bins), function(x) if(sum(binDetect[binDetect$bins == x, -1]>= detectThresh) >= 4) x))) {
# function to loop through all pairs of genes in a MAG and correlate across indiviuals' samples in which the MAG was detected by some threshold
# Ideally this can show how much linkage is within a mouth vs between mouths by the shape of the curve
corVec <- c()
for (MAG in MAGs) {
mouths_detected_in <- table(gsub("[0-9]*","", colnames(binDetect)[which(binDetect[binDetect$bins == MAG,] >= detectThresh)]))
mouths_detected_in <- names(mouths_detected_in)[mouths_detected_in >= 4]
if (!byIndiv & (length(mouths_detected_in) > 0)){
magCor <- cor(t(geneCovs[geneCovs$key %in% geneInfo$gene_callers_id[geneInfo$bin == MAG],
grep(paste0(mouths_detected_in, collapse = "|"), colnames(geneCovs))]), method = corMeth)
magCor <- na.omit(magCor[lower.tri(magCor)])
corVec <- c(corVec, magCor)
} else {
for (mouth in mouths_detected_in) {
magCor <- cor(t(geneCovs[geneCovs$key %in% geneInfo$gene_callers_id[geneInfo$bin == MAG], grep(mouth, colnames(geneCovs))]), method = corMeth)
magCor <- na.omit(magCor[lower.tri(magCor)])
corVec <- c(corVec, magCor)
}
}
}
return(corVec)
}
correlateGenesAcrossMAGs <- function(geneCovs, geneInfo, binDetect, binCovs, byIndiv=T, detectThresh=0.75, excludeOwn=F, corMeth='pearson',
MAGs=unlist(sapply(levels(binDetect$bins), function(x) if(sum(binDetect[binDetect$bins == x, -1]>= detectThresh) >= 4) x))){
# subset to focal MAGs for speed if diong a subset
geneCovs <- geneCovs[geneCovs$key %in% geneInfo[geneInfo$bin %in% MAGs,"gene_callers_id"],]
corDF <- data.frame(gene=as.character(geneCovs$key), mouths_detected_in=NA, bin=NA, p.value=NA, correlation=NA)
corDF$bin <- as.character(corDF$bin)
corDF$mouths_detected_in <- as.character(corDF$mouths_detected_in)
rownames(corDF) <- as.character(geneCovs$key)
for (MAG in MAGs) {
# subset to grab just mouths where this MAG is at least detect_thresh in 4 samples
mouths_detected_in <- table(gsub("[0-9]*","", colnames(binDetect)[which(binDetect[binDetect$bins == MAG,] >= detectThresh)]))
mouths_detected_in <- names(mouths_detected_in)[mouths_detected_in >= 4]
print(mouths_detected_in)
binTrace <- as.numeric(binCovs[binCovs$bins == MAG, grep(paste0(mouths_detected_in, collapse = "|"), colnames(binCovs))])
# for each gene in this MAG's genes
for (gene in unique(as.character(geneCovs$key[geneCovs$key %in% geneInfo$gene_callers_id[geneInfo$bin == MAG]]))){
geneTrace <- as.numeric(geneCovs[geneCovs$key == gene, grep(paste0(mouths_detected_in, collapse = "|"), colnames(geneCovs))])
geneCor <- cor.test(geneTrace, binTrace, method = corMeth, exact = F)
corDF[gene,c("p.value","correlation")] <- c(p.value=geneCor$p.value, correlation=geneCor$estimate)
corDF$mouths_detected_in[corDF$gene == gene] <- paste0(mouths_detected_in, collapse = "|")
corDF$bin[corDF$gene == gene] <- MAG
}
}
return(corDF)
}
fetchVariabilityCDNForGenes <- function(genes=c("3555"), geneInfo){
bin <- unique(as.character(geneInfo$bin[geneInfo$gene_callers_id %in% genes]))
varData <- read.csv(paste0('catalyst_idba/codon_var/', bin, '-var-CDN.tsv'), sep = "\t")
varData <- varData[varData$corresponding_gene_call %in% genes,]
return(varData)
}
codonSynonByMAG <- function(MAG_list, geneInfo, geneCovs, binCovs, corMat=geneVsMAG, corThresh=0.2, minMedCov=2, greaterThan=F){
corMat <- corMat[!is.na(corMat$correlation),]
all_syn <- data.frame(MAG=character(), synonymity=numeric())
for (MAG in MAG_list){
if (greaterThan) {
genes_not_correlated <- as.character(corMat$gene[abs(corMat$correlation) >= corThresh & (corMat$bin == MAG)])
} else {
genes_not_correlated <- as.character(corMat$gene[abs(corMat$correlation) <= corThresh & (corMat$bin == MAG)])
}
print(paste0(length(genes_not_correlated), " genes meeting criteria (", MAG,")"))
if (!(MAG %in% gsub("-var-CDN.tsv", "", list.files("catalyst_idba/codon_var/")))) {
print(paste0(MAG, " has no associated variability file..."))
next()
}
if (length(genes_not_correlated) > 0){ # some MAGs might not have any genes meeting criteria
varinfo <- fetchVariabilityCDNForGenes(genes = genes_not_correlated, geneInfo = geneInfo)
varinfo <- varinfo[grep(unique(corMat$mouths_detected_in[corMat$bin == MAG]), varinfo$sample_id),]
genes_with_var <- varinfo$corresponding_gene_call
mag_syn <- varinfo$synonymity
gene_lengths <- varinfo$gene_length
sample_id <- varinfo$sample_id
if (length(mag_syn) == 0){
print("No codon information for any of ^these genes")
next()
}
all_syn <- rbind(all_syn, data.frame(MAG=MAG, synonymity=mag_syn, gene=genes_with_var, gene_length=gene_lengths, sample=sample_id))
}
}
return(all_syn)
}
Reading the processed data
#########################
## DATA ##
#########################
gene_info <- read.csv('catalyst_idba/NON-REDUNDANT-MAGs-gene_info.txt', sep="\t", header = T)
gene_covs <- read.csv('catalyst_idba/NON-REDUNDANT-MAGs-GENE-COVERAGES.txt', sep = "\t", header = T)
gene_detect <- read.csv('catalyst_idba/NON-REDUNDANT-MAGs-GENE-DETECTION.txt', sep = "\t", header = T)
bin_covs <- read.csv('catalyst_idba/NON-REDUNDANT-MAGs-SUMMARY/bins_across_samples/mean_coverage.txt', sep = "\t", header =T)
bin_detect <- read.csv('catalyst_idba/NON-REDUNDANT-MAGs-SUMMARY/bins_across_samples/detection.txt', sep = "\t", header =T)
bin_stats <- read.csv('catalyst_idba/REDUNDANT-MAGs-STATS.txt', sep = "\t", header = T)
eggnog <- read.csv('catalyst_idba/eggnog_annotations.tsv', sep = "\t" , header = T)
rownames(eggnog) <- eggnog$query_name
gtdbtk <- makegtdbtktaxMat()
rownames(gene_info) <- gene_info$gene_callers_id
gene_covs <- gene_covs[-46594,] # remove duplicate gene caller id from AN as it is an AO gene according to gene_info
gene_detect <- gene_detect[-46594,] # same for gene_detect
rownames(gene_covs) <- gene_covs$key
Final figures
S. Figure 1
Overall MAG taxonomy:
gtdbtk$Genus <- factor(gtdbtk$Genus, levels = names(table(gtdbtk$Genus))[order(table(gtdbtk$Genus))])
ggplot(gtdbtk, aes(x=Genus, fill=Genus, group=Genus)) +
geom_histogram(stat = 'count', position = 'identity') + labs(x="Genus", y = "MAG Count") +
scale_fill_manual(values = c(grey.colors(15,rev = T),rev(dans_colors)), guide=F) + coord_flip() +
scale_y_continuous(expand = c(0,0)) + theme_basic +
theme(panel.spacing = unit(0.1, 'line'), strip.background = element_rect(fill='#EaEaEa'), axis.ticks.x = element_line(color='black'), axis.text.y = element_text(size=6),
panel.grid.minor.x = element_line(), panel.grid.major.x = element_line())
Fig 1A
Cartoon of how the anvi’o metric ‘detection’ (horizontal coverage) works, for Figure 1A
plotDetectionCartoon <- function(gene_data, geneInfo, focalMAGs, log=F, gtdbtk=makegtdbtktaxMat(), colorGenes=NA) {
gene_data <- gene_data[gene_data$key %in% geneInfo$gene_callers_id[geneInfo$bin %in% focalMAGs],]
gene_data$MAG <- geneInfo$bin[match(gene_data$key, geneInfo$gene_callers_id)]
df <- melt(gene_data, id.vars = c("key","MAG"))
colnames(df)[colnames(df) == "variable"] <- "sample"
df$sample <- factor(df$sample, levels = mixedsort(unique(as.character(df$sample))))
df$value <- as.numeric(as.character(df$value))
df$key <- factor(df$key, levels = sort(unique(as.numeric(as.character(df$key)))))
df$MAG <- factor(df$MAG, levels = sort(unique(as.character(df$MAG))))
df <- df[order(df$MAG, df$sample, df$key),]
if (log) {
df$value <- log10(df$value + 1)
}
geneColors <- setNames(rep('black',length(levels(df$key))), levels(df$key))
geneColors[colorGenes] <- 'red'
df$detected <- "not"
df$detected[df$value > 0] <- "yes"
detectDF <- df[cumsum(rle(as.character(df$detected))$lengths),]
p <- ggplot(na.omit(df), aes(x=key, y = value)) +
geom_segment(aes(x=key, xend=key, y=0, yend=Inf, color=detected)) +
geom_col(fill='black', color='black', size=0.1)
if (log) {
p <- p + facet_grid(sample~MAG, scales='free_y')
} else {
p <- p + facet_grid(sample~MAG, scales = 'free_y')
}
p <- p + #scale_fill_manual(values = geneColors) +
scale_color_manual(values = setNames(c("#Ba52B2","#FfB6E1"),c("yes","not"))) +
labs(y = "Coverage", x = "Gene") + scale_y_continuous(n.breaks = 3) +
theme_basic + theme(legend.position = 'none') +
theme(panel.border = element_blank(), axis.text.x = element_blank(), axis.ticks.y = element_blank(), panel.spacing = unit(0, "lines"),
axis.text.y = element_text(size=5), strip.text.y = element_text(size=4))
return(p)
}
plotDetectionCartoon(gene_covs[,grep("key|AB1|AC|AI15", colnames(gene_covs))], gene_info, focalMAGs = 'FINAL_AI_MAG_00008', log = T)
Figure 1 C
plotDiversityTime(matForDiv = bin_covs[,2:ncol(bin_covs)], title = "MAGs", method = 'simpson')
## div samp indiv day dataset
## AB1 0.9808752 AB1 AB 1 MAGs
## AB15 0.9820555 AB15 AB 15 MAGs
## AB2 0.9684930 AB2 AB 2 MAGs
## AB3 0.9743385 AB3 AB 3 MAGs
## AB5 0.9708051 AB5 AB 5 MAGs
## AB6 0.9744770 AB6 AB 6 MAGs
Fig. 2BC, SFig 5
We can plot the proportion of each MAG explained by each haplotype for each metagenome (e.g., Figure 2BC) or loop through them all (Supplemental Figure 5) by uncommenting the for loop:
MAG <- "FINAL_AI_MAG_00008"
#for (MAG in gsub("-Bacteria.*$","", gsub("binDESMAN/","",desStrains$base_pred_file))) {
numStrains <- desStrains$nHap[grep(MAG,desStrains$base_pred_file)]
if (length(list.files(paste0("catalyst_idba/mag_bacteria71_genes/", MAG, "-Bacteria_71-SNVs-DESMAN-FINAL_",numStrains,"_0/"))) < 2) {
# that means the expected files aren't there, print that we're missing and skip
print(paste0("Missing files for ", MAG, " with n_hap=", numStrains))
next()
}
desmanFreq <- read.csv(paste0("catalyst_idba/mag_bacteria71_genes/", MAG, "-Bacteria_71-SNVs-DESMAN-FINAL_",numStrains,"_0/Gamma_mean.csv"), header = T,
col.names = c("sample", paste0("hap", 1:as.numeric(numStrains))))
desmanFreq <- melt(desmanFreq)
colnames(desmanFreq) <- c("sample","haplotype","frequency")
desmanFreq$indiv <- factor(gsub("[0-9]*","", desmanFreq$sample))
desmanFreq$day <- factor(gsub("A[A-Z]", "", desmanFreq$sample), levels = mixedsort(unique(gsub("A[A-Z]", "", desmanFreq$sample))))
desmanFreq$category <- 'haplotypes'
magCov <- bin_covs[bin_covs$bins == MAG,2:ncol(bin_covs)]
magDF <- data.frame(sample=names(magCov), frequency=as.numeric(magCov)/max(as.numeric(magCov)))
magDF$indiv <- factor(gsub("[0-9]*","", magDF$sample))
magDF <- magDF[magDF$indiv %in% desmanFreq$indiv,]
magDF$day <- factor(gsub("A[A-Z]", "", magDF$sample), levels = mixedsort(unique(gsub("A[A-Z]", "", as.character(magDF$sample)))))
magDF$category <- 'aMAG'
p <- ggplot(desmanFreq, aes(x=day, y = frequency)) + ggtitle(gsub("(^.*)(MAG_)(.*$)","\\2\\1\\3", gsub("_000","_", gsub("FINAL_","",MAG)))) +
geom_col(aes(fill=haplotype)) +
geom_col(data=magDF, fill='black', aes(group=indiv)) + facet_grid(category~indiv, space = 'free_x', scales='free_x') + # SINGLE
#geom_line(data=magDF, stat = 'identity', color=NA, aes(group=indiv)) + facet_grid(.~indiv, space = 'free_x', scales='free_x') + # LOOP FOR ALL
scale_y_continuous(expand = c(0,0), sec.axis = sec_axis(~. * max(as.numeric(magCov)), breaks =function(x) c(0, floor(x[2])))) +
scale_x_discrete(expand = c(0,0)) + labs(x='Sample Day', y = 'Haplotype proportion') +
scale_fill_manual(values = c(beyonce_palette(66)[seq(from = 1, to = 11, length.out = min(11,length(levels(desmanFreq$haplotype))))], rainbow(9))) +
theme_classic() + theme(axis.ticks.x = element_blank(), axis.text = element_text(color = 'black'), plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5), axis.text.x = element_text(color = 'black', angle = 90, hjust = 1, vjust = 0.5),
panel.spacing = unit(4/length(levels(desmanFreq$indiv)), 'lines'))
# for single:
g <- ggplot_gtable(ggplot_build(p))
g$heights[8] <- .2*g$heights[8] # set the top facet row to be 20% of the height
grid.draw(g)
# for loop:
#print(p)
#}
Supplemental Figure 6A
And then we can assess the genes correlated with a specific haplotype, as in Supplemental Figure 6A:
MAG <- "FINAL_AI_MAG_00008"
numStrains <- desStrains$nHap[grep(MAG,desStrains$base_pred_file)]
desmanFreq <- read.csv(paste0("catalyst_idba/mag_bacteria71_genes/", MAG, "-Bacteria_71-SNVs-DESMAN-FINAL_",numStrains,"_0/Gamma_mean.csv"), header = T,
col.names = c("sample", paste0("hap", 1:as.numeric(numStrains))))
# get the genes from the samples the haplotypes are detected in (might be different from coverage due to SNV coverage depth)
relevantGeneSamps <- gene_covs[gene_covs$key %in% gene_info$gene_callers_id[gene_info$bin == MAG], c('key', as.character(desmanFreq$sample))]
hapCor <- cor(desmanFreq[,2:ncol(desmanFreq)],t(relevantGeneSamps[,2:ncol(relevantGeneSamps)]))
colnames(hapCor) <- as.character(relevantGeneSamps$key)
hapCorFilled <- hapCor
hapCorFilled[is.na(hapCorFilled)] <- 0
geneOrder <- hclust(vegdist(t(hapCorFilled), method = 'euc'), method = 'complete')
hapWide <- melt(hapCor)
colnames(hapWide) <- c("haplotype","gene","correlation")
hapWide$gene <- factor(hapWide$gene, levels = geneOrder$labels[geneOrder$order])
# genes associated with Hap4
hap4genes <- hapCor[,(hapCor["hap2",] > 0.8) & apply(hapCor[paste0("hap",c(1,3,4)),], 2, function(x) all(x < 0.2))]
ggplot(hapWide, aes(x=gene, y = correlation, group=haplotype)) +
geom_line(aes(color=haplotype)) +
scale_y_continuous(expand = c(0,0)) + scale_x_discrete(expand = c(0,0)) +
scale_color_manual(values = c(beyonce_palette(66)[seq(from = 1, to = 11, length.out = min(11,length(levels(hapWide$haplotype))))], rainbow(9))) +
theme(axis.ticks.x = element_blank(), axis.text = element_text(color = 'black'), plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5), strip.background = element_rect(fill = '#dcdcdc', colour = NA),
axis.text.x = element_blank(), axis.line = element_line())
Supplemental Figure 6C
ai_08_gcs <- droplevels(gene_info[gene_info$bin == "FINAL_AI_MAG_00008",])
ai_08_gcs$association <- 'Other'
ai_08_gcs$association[ai_08_gcs$gene_callers_id %in% colnames(hap4genes)] <- 'Haplotype 4'
ai_08_gcs$contig <- factor(ai_08_gcs$contig, levels=names(sort(table(ai_08_gcs$contig), decreasing = F)))
ggplot(ai_08_gcs) +
geom_bar(aes(y=contig, fill=association), position='stack', width=0.7) +
scale_x_continuous(expand = c(0,0)) + labs(y=paste0("Contigs in MAG_AI_08 (n=", length(levels(ai_08_gcs$contig)),")"), x="Number of genes", fill="Associated with:") +
scale_fill_manual(values = c('#E74690','grey50')) + guides(fill=guide_legend(nrow=2)) +
theme_basic + theme(axis.text.y = element_blank(), legend.position = c(.6, .15), legend.background = element_rect(color='black', size=0.25),
axis.ticks.x = element_line(), axis.ticks.y = element_blank())
Fig. 3A and S. Fig 7
This block generates figure 3A and the corresponding supplemental
figures 7 and 8. The first plot generated is SFig. 7, which shows the
distribution of detrended coveage (relative to the parent MAG) for each
gene in MAG_AD_09 by donor. This confirms that genes above 5 standard
deviations are truly outliers. Figure 3A then plots the coverage of each
gene (thin line) in the MAG for each timepoint along with the MAG
coverage (thick line). Colored lines mark the genes with more than 5
standard deviations above the average detrended coverage by donor.
# supplemental figure 7: sd distributions for MAG_AD_09
ad09_5sd <- detrendDistribution(gene_info, gene_covs, bin_covs, bin="FINAL_AD_MAG_00009", bin_stats, plot_it = T,
focusGenes = NULL, sdMultiple = 5)
## [1] "1439 genes in FINAL_AD_MAG_00009"
# Figure 3A - gene traces
traceGenesByMAG(gene_info, gene_covs, bin_covs, bin="FINAL_AD_MAG_00009", bin_stats, facetGrid = T, log10_transform = F, sqrt_transform = T, focusColor='#Fa6a3a',
focusGenes = as.character(ad09_5sd$Gene[ad09_5sd$outlierInParticipant]),
focusIn=as.character(ad09_5sd$Participant[ad09_5sd$outlierInParticipant]))
## [1] "1439 genes in FINAL_AD_MAG_00009"
## [1] "More than 20 genes are here, skipping legend and not dumping functions"
We then generated Supplemental Table 3 (Pfam annotations for these highlighted outlier genes):
length(unique(ad09_5sd$Gene[ad09_5sd$outlierInParticipant]))
## [1] 48
setNames(as.character(gene_info[unique(as.character(ad09_5sd$Gene[ad09_5sd$outlierInParticipant])),"Pfam"]), unique(as.character(ad09_5sd$Gene[ad09_5sd$outlierInParticipant])))
s_table3 <-cbind(gene=unique(as.character(ad09_5sd$Gene[ad09_5sd$outlierInParticipant])),
Pfam = as.character(gene_info[unique(as.character(ad09_5sd$Gene[ad09_5sd$outlierInParticipant])),"Pfam"]),
eggnog[unique(as.character(ad09_5sd$Gene[ad09_5sd$outlierInParticipant])), c("COG_Category", "Preferred_name","KEGG_ko","EC"),])
print(head(s_table3))
## gene Pfam COG_Category
## NA 165278 <NA>
## 165279 165279 S
## 165031 165031 Bacterial regulatory proteins, tetR family K
## 164689 164689 TfoX N-terminal domain K
## NA.1 165280 <NA>
## 164484 164484 Glycosyl transferase family 2 S
## Preferred_name KEGG_ko EC
## NA <NA> <NA> <NA>
## 165279
## 165031
## 164689 ko:K07343
## NA.1 <NA> <NA> <NA>
## 164484 ycbB ko:K00786
Figure 3BC
First we looped through all MAGs to get the total outlier genes from the entire dataset:
mags_to_check <- unlist(sapply(levels(bin_detect$bins), function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.50) >= 15 ) x))
length(mags_to_check)
## [1] 199
for (mag in mags_to_check) {
single_outliers <- detrendDistribution(gene_info, gene_covs, bin_covs, bin=mag, bin_stats, plot_it=F,
focusGenes = NULL, sdMultiple = 5)
single_outliers <- single_outliers[single_outliers$outlierInParticipant,]
if (mag == mags_to_check[1]){
mag_outliers <- data.frame(MAG=mag, Gene=unique(as.character(single_outliers$Gene)))
} else {
mag_outliers <- rbind(mag_outliers, data.frame(MAG=mag, Gene=unique(as.character(single_outliers$Gene))))
}
}
Then this plotting function was defined to display the data:
processOutlierFunctions <- function(outliers=mag_outliers, func_category='COG_Category', minCount=3, dropBlank=T){
if (func_category %in% colnames(gene_info)){
outliers[,func_category] <- gene_info[as.character(outliers$Gene),func_category]
} else {
outliers[,func_category] <- eggnog[as.character(outliers$Gene),func_category]
}
if (dropBlank){
outliers <- outliers[!(outliers[,func_category] %in% c("",NA)),]
}
outliers <- outliers %>% mutate(count=1) %>% group_by(.dots=func_category) %>% summarise(mean=sum(count), sd=sd(count))
outliers <- as.data.frame(outliers)
outliers <- outliers[outliers$mean >= minCount,]
outliers[,func_category] <- as.character(outliers[,func_category])
outliers[outliers[,func_category] %in% names(cogValues),func_category] <- cogValues[outliers[,func_category][outliers[,func_category] %in% names(cogValues)]]
outliers[,func_category] <- factor(outliers[,func_category], levels = unique(as.character(outliers[order(outliers$mean, decreasing = F),func_category])))
p <- ggplot(outliers, aes_string(y=func_category, x = "mean")) +
geom_col(fill='#a0dff8ff') + labs(x="Count in all outliers", y=gsub("_"," ", func_category)) +
scale_x_continuous(expand = c(0,0)) +
theme_classic() + theme(axis.ticks = element_line(color='black'), panel.grid.major.x = element_line())
print(p)
return(outliers)
}
With this function Figure 3B was plotted like so:
ad09_5sd$MAG <- "FINAL_AD_MAG_00009"
x <- processOutlierFunctions(unique(ad09_5sd[ad09_5sd$outlierInParticipant,c("MAG","Gene")]), func_category = 'COG_Category', minCount = 1)
For Figure 3C, for brevity we only plotted the COG categories above 60 and then cropped the output plot.
x <- processOutlierFunctions(mag_outliers, func_category = 'Pfam', minCount =60, dropBlank=T)
## `summarise()` ungrouping output (override with `.groups` argument)
S. Fig. 8
To assess the overall coupling (or decoupling) of genes, we computed all possible pairwise correlations for MAGs meeting a detection criterion of >0.75 in varying numbers of samples:
corYesByIndiv30 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = T,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 30) x)))
corYesByIndiv20 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = T,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 20) x)))
corYesByIndiv40 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = T,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 40) x)))
geneVarStats <- function(df=corNotByIndiv40, subsampleTo=1E6){
med <- median(df)
ave <- mean(df)
sdiv <- sd(df)
ma <- max(df)
mi <- min(df)
quants <- quantile(df)
return(data.frame(median=med, average=ave, max=ma, min=mi, sd=sdiv, q1=quants[2], q3=quants[4], correlations=sample(df, subsampleTo)))
}
cors <- rbind(cbind(indiv="ByIndividual", numIndiv=40, geneVarStats(corYesByIndiv40)),
cbind(indiv="ByIndividual", numIndiv=30, geneVarStats(corYesByIndiv30)),
cbind(indiv="ByIndividual", numIndiv=20, geneVarStats(corYesByIndiv20)))
rm(corYesByIndiv40, corYesByIndiv30, corYesByIndiv20) # to free up memory
#corYesByIndiv10 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = T,
# MAGs=unlist(sapply(levels(bin_detect$bins),
# function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 10) x)))
#cors <- rbind(cors, cbind(indiv="ByIndividual", numIndiv=10, geneVarStats(corYesByIndiv10)))
#rm(corYesByIndiv10)
corNotByIndiv40 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = F,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 40) x)))
corNotByIndiv30 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = F,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 30) x)))
corNotByIndiv20 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = F,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 20) x)))
corNotByIndiv10 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = F,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 10) x)))
corNotByIndiv5 <- correlateMagGenePairs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, detectThresh = 0.75, byIndiv = F,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 5) x)))
cors <- rbind(cors,
cbind(indiv="NotByIndividual", numIndiv=40, geneVarStats(corNotByIndiv40)),
cbind(indiv="NotByIndividual", numIndiv=30, geneVarStats(corNotByIndiv30)),
cbind(indiv="NotByIndividual", numIndiv=20, geneVarStats(corNotByIndiv20)),
cbind(indiv="NotByIndividual", numIndiv=10, geneVarStats(corNotByIndiv10)),
cbind(indiv="NotByIndividual", numIndiv=5, geneVarStats(corNotByIndiv5)))
Final prep and plotting:
densiCor <- cors %>% group_by(indiv, numIndiv, q1, q3, median) %>% summarise(densY= list(density(correlations)$y), densX=list(density(correlations)$x)) %>% unnest(c(densY,densX))
densiCor$q1q3 <- (densiCor$densX >= densiCor$q1) * (densiCor$densX <= densiCor$q3)
densiCor$q1q3[densiCor$densX >= densiCor$q3] <- 3
densiCor <- densiCor %>% group_by(indiv, numIndiv) %>% mutate(medianY=densY[which.min(abs(densX-median))])
ggplot(densiCor[abs(densiCor$densX) <= 1,], aes(x=densX, y=densY, color=indiv)) +
geom_line() +
geom_ribbon(aes(ymin=0, ymax=densY, fill=interaction(indiv,factor(q1q3)))) +
geom_segment(data=densiCor[!duplicated(densiCor[,-c(6,7)]),], y=0, aes(x=median, xend=median, yend=medianY, color=indiv)) +
facet_grid(numIndiv~., switch = 'y', scales='free_y', space = 'free_y') +
scale_x_continuous(expand = c(0,0), limits = c(-1,1)) + scale_y_continuous(position = 'right', expand = c(0,0)) +
scale_color_manual(values = c('#00AaFa','#Fa6a3a44'), labels=c("Mouth-by-mouth", "Across mouths")) +
scale_fill_manual(values = c(NA,NA,'#00AaFa44','#Fa6a3a44',NA,NA), breaks = c(3,4), labels=c("Mouth-by-mouth", "Across mouths")) +
labs(x="Pearson correlation coefficient", sec.y = "Number of samples detected in\n(at least 75% of MAG's nucleotides covered)",
fill="Correlation\nscope", color="Correlation\nscope", y='Density') +
theme_bw() + theme_basic + theme(axis.ticks.x = element_line())
Figure 4
We then turned to look at nucleotide variability data to assess the frequency of single codon variants, specifically whether decoupled genes were enriched in non-synonymous codon variants (indicative of phenomena like HGT into distinct populations), and whether non-synonymous codons sweep through a population over the two-week windoe or remain stable.
geneVsMAG <- correlateGenesAcrossMAGs(geneCovs = gene_covs, binDetect = bin_detect, geneInfo = gene_info, binCovs = bin_covs, detectThresh = 0.75, byIndiv = F,
MAGs=unlist(sapply(levels(bin_detect$bins),
function(x) if(sum(bin_detect[bin_detect$bins == x, -1]>= 0.75) >= 40) x)))
quantile(abs(geneVsMAG$correlation), probs = c(0.01,0.1,0.9,0.99), na.rm = T) # USE TO GET GOOD VALUES for next step!
## 1% 10% 90% 99%
## 0.03641639 0.41827478 0.98533270 0.99720204
synUnlinked01 <- codonSynonByMAG(MAG_list = unique(as.character(geneVsMAG$bin)), corMat = geneVsMAG, gene_info, gene_covs, bin_covs, corThresh = 0.0364, greaterThan=F)
## [1] "0 genes meeting criteria (FINAL_AT_MAG_00012)"
## [1] "4 genes meeting criteria (FINAL_AJ_MAG_00005)"
## [1] "FINAL_AJ_MAG_00005 has no associated variability file..."
## [1] "14 genes meeting criteria (FINAL_AE_MAG_00012)"
## [1] "97 genes meeting criteria (FINAL_AC_MAG_00033)"
## [1] "FINAL_AC_MAG_00033 has no associated variability file..."
## [1] "43 genes meeting criteria (FINAL_AN_MAG_00021)"
## [1] "FINAL_AN_MAG_00021 has no associated variability file..."
## [1] "20 genes meeting criteria (FINAL_AD_MAG_00009)"
## [1] "0 genes meeting criteria (FINAL_AC_MAG_00020)"
## [1] "14 genes meeting criteria (FINAL_AE_MAG_00015)"
## [1] "1 genes meeting criteria (FINAL_AC_MAG_00006)"
## [1] "No codon information for any of ^these genes"
## [1] "63 genes meeting criteria (FINAL_AQ_MAG_00003)"
## [1] "2 genes meeting criteria (FINAL_AD_MAG_00003)"
## [1] "No codon information for any of ^these genes"
## [1] "9 genes meeting criteria (FINAL_AJ_MAG_00003)"
## [1] "18 genes meeting criteria (FINAL_AI_MAG_00002)"
## [1] "FINAL_AI_MAG_00002 has no associated variability file..."
## [1] "32 genes meeting criteria (FINAL_AR_MAG_00012)"
## [1] "FINAL_AR_MAG_00012 has no associated variability file..."
## [1] "73 genes meeting criteria (FINAL_AE_MAG_00007)"
## [1] "FINAL_AE_MAG_00007 has no associated variability file..."
## [1] "8 genes meeting criteria (FINAL_AJ_MAG_00008)"
## [1] "7 genes meeting criteria (FINAL_AC_MAG_00034)"
## [1] "FINAL_AC_MAG_00034 has no associated variability file..."
## [1] "4 genes meeting criteria (FINAL_AD_MAG_00002)"
## [1] "0 genes meeting criteria (FINAL_AC_MAG_00008)"
## [1] "7 genes meeting criteria (FINAL_AC_MAG_00028)"
## [1] "10 genes meeting criteria (FINAL_AD_MAG_00021)"
## [1] "FINAL_AD_MAG_00021 has no associated variability file..."
## [1] "0 genes meeting criteria (FINAL_AD_MAG_00015)"
## [1] "22 genes meeting criteria (FINAL_AH_MAG_00007)"
## [1] "12 genes meeting criteria (FINAL_AQ_MAG_00010)"
## [1] "3 genes meeting criteria (FINAL_AD_MAG_00004)"
## [1] "FINAL_AD_MAG_00004 has no associated variability file..."
synUnlinked10 <- codonSynonByMAG(MAG_list = unique(as.character(geneVsMAG$bin)), corMat = geneVsMAG, gene_info, gene_covs, bin_covs, corThresh = 0.418, greaterThan=F)
## [1] "41 genes meeting criteria (FINAL_AT_MAG_00012)"
## [1] "116 genes meeting criteria (FINAL_AJ_MAG_00005)"
## [1] "FINAL_AJ_MAG_00005 has no associated variability file..."
## [1] "111 genes meeting criteria (FINAL_AE_MAG_00012)"
## [1] "302 genes meeting criteria (FINAL_AC_MAG_00033)"
## [1] "FINAL_AC_MAG_00033 has no associated variability file..."
## [1] "377 genes meeting criteria (FINAL_AN_MAG_00021)"
## [1] "FINAL_AN_MAG_00021 has no associated variability file..."
## [1] "195 genes meeting criteria (FINAL_AD_MAG_00009)"
## [1] "103 genes meeting criteria (FINAL_AC_MAG_00020)"
## [1] "249 genes meeting criteria (FINAL_AE_MAG_00015)"
## [1] "30 genes meeting criteria (FINAL_AC_MAG_00006)"
## [1] "331 genes meeting criteria (FINAL_AQ_MAG_00003)"
## [1] "123 genes meeting criteria (FINAL_AD_MAG_00003)"
## [1] "233 genes meeting criteria (FINAL_AJ_MAG_00003)"
## [1] "351 genes meeting criteria (FINAL_AI_MAG_00002)"
## [1] "FINAL_AI_MAG_00002 has no associated variability file..."
## [1] "196 genes meeting criteria (FINAL_AR_MAG_00012)"
## [1] "FINAL_AR_MAG_00012 has no associated variability file..."
## [1] "503 genes meeting criteria (FINAL_AE_MAG_00007)"
## [1] "FINAL_AE_MAG_00007 has no associated variability file..."
## [1] "140 genes meeting criteria (FINAL_AJ_MAG_00008)"
## [1] "52 genes meeting criteria (FINAL_AC_MAG_00034)"
## [1] "FINAL_AC_MAG_00034 has no associated variability file..."
## [1] "37 genes meeting criteria (FINAL_AD_MAG_00002)"
## [1] "284 genes meeting criteria (FINAL_AC_MAG_00008)"
## [1] "196 genes meeting criteria (FINAL_AC_MAG_00028)"
## [1] "276 genes meeting criteria (FINAL_AD_MAG_00021)"
## [1] "FINAL_AD_MAG_00021 has no associated variability file..."
## [1] "5 genes meeting criteria (FINAL_AD_MAG_00015)"
## [1] "138 genes meeting criteria (FINAL_AH_MAG_00007)"
## [1] "76 genes meeting criteria (FINAL_AQ_MAG_00010)"
## [1] "159 genes meeting criteria (FINAL_AD_MAG_00004)"
## [1] "FINAL_AD_MAG_00004 has no associated variability file..."
synLinked90 <- codonSynonByMAG(MAG_list = unique(as.character(geneVsMAG$bin)), corMat = geneVsMAG, gene_info, gene_covs, bin_covs, corThresh = 0.985, greaterThan=T)
## [1] "1206 genes meeting criteria (FINAL_AT_MAG_00012)"
## [1] "118 genes meeting criteria (FINAL_AJ_MAG_00005)"
## [1] "FINAL_AJ_MAG_00005 has no associated variability file..."
## [1] "7 genes meeting criteria (FINAL_AE_MAG_00012)"
## [1] "79 genes meeting criteria (FINAL_AC_MAG_00033)"
## [1] "FINAL_AC_MAG_00033 has no associated variability file..."
## [1] "911 genes meeting criteria (FINAL_AN_MAG_00021)"
## [1] "FINAL_AN_MAG_00021 has no associated variability file..."
## [1] "83 genes meeting criteria (FINAL_AD_MAG_00009)"
## [1] "108 genes meeting criteria (FINAL_AC_MAG_00020)"
## [1] "No codon information for any of ^these genes"
## [1] "2 genes meeting criteria (FINAL_AE_MAG_00015)"
## [1] "No codon information for any of ^these genes"
## [1] "1 genes meeting criteria (FINAL_AC_MAG_00006)"
## [1] "No codon information for any of ^these genes"
## [1] "1 genes meeting criteria (FINAL_AQ_MAG_00003)"
## [1] "3 genes meeting criteria (FINAL_AD_MAG_00003)"
## [1] "2 genes meeting criteria (FINAL_AJ_MAG_00003)"
## [1] "7 genes meeting criteria (FINAL_AI_MAG_00002)"
## [1] "FINAL_AI_MAG_00002 has no associated variability file..."
## [1] "62 genes meeting criteria (FINAL_AR_MAG_00012)"
## [1] "FINAL_AR_MAG_00012 has no associated variability file..."
## [1] "0 genes meeting criteria (FINAL_AE_MAG_00007)"
## [1] "FINAL_AE_MAG_00007 has no associated variability file..."
## [1] "368 genes meeting criteria (FINAL_AJ_MAG_00008)"
## [1] "166 genes meeting criteria (FINAL_AC_MAG_00034)"
## [1] "FINAL_AC_MAG_00034 has no associated variability file..."
## [1] "437 genes meeting criteria (FINAL_AD_MAG_00002)"
## [1] "0 genes meeting criteria (FINAL_AC_MAG_00008)"
## [1] "2 genes meeting criteria (FINAL_AC_MAG_00028)"
## [1] "300 genes meeting criteria (FINAL_AD_MAG_00021)"
## [1] "FINAL_AD_MAG_00021 has no associated variability file..."
## [1] "235 genes meeting criteria (FINAL_AD_MAG_00015)"
## [1] "93 genes meeting criteria (FINAL_AH_MAG_00007)"
## [1] "478 genes meeting criteria (FINAL_AQ_MAG_00010)"
## [1] "115 genes meeting criteria (FINAL_AD_MAG_00004)"
## [1] "FINAL_AD_MAG_00004 has no associated variability file..."
synLinked99 <- codonSynonByMAG(MAG_list = unique(as.character(geneVsMAG$bin)), corMat = geneVsMAG, gene_info, gene_covs, bin_covs, corThresh = 0.997, greaterThan=T)
## [1] "293 genes meeting criteria (FINAL_AT_MAG_00012)"
## [1] "4 genes meeting criteria (FINAL_AJ_MAG_00005)"
## [1] "FINAL_AJ_MAG_00005 has no associated variability file..."
## [1] "0 genes meeting criteria (FINAL_AE_MAG_00012)"
## [1] "0 genes meeting criteria (FINAL_AC_MAG_00033)"
## [1] "FINAL_AC_MAG_00033 has no associated variability file..."
## [1] "61 genes meeting criteria (FINAL_AN_MAG_00021)"
## [1] "FINAL_AN_MAG_00021 has no associated variability file..."
## [1] "0 genes meeting criteria (FINAL_AD_MAG_00009)"
## [1] "0 genes meeting criteria (FINAL_AC_MAG_00020)"
## [1] "0 genes meeting criteria (FINAL_AE_MAG_00015)"
## [1] "0 genes meeting criteria (FINAL_AC_MAG_00006)"
## [1] "0 genes meeting criteria (FINAL_AQ_MAG_00003)"
## [1] "0 genes meeting criteria (FINAL_AD_MAG_00003)"
## [1] "0 genes meeting criteria (FINAL_AJ_MAG_00003)"
## [1] "0 genes meeting criteria (FINAL_AI_MAG_00002)"
## [1] "FINAL_AI_MAG_00002 has no associated variability file..."
## [1] "0 genes meeting criteria (FINAL_AR_MAG_00012)"
## [1] "FINAL_AR_MAG_00012 has no associated variability file..."
## [1] "0 genes meeting criteria (FINAL_AE_MAG_00007)"
## [1] "FINAL_AE_MAG_00007 has no associated variability file..."
## [1] "0 genes meeting criteria (FINAL_AJ_MAG_00008)"
## [1] "1 genes meeting criteria (FINAL_AC_MAG_00034)"
## [1] "FINAL_AC_MAG_00034 has no associated variability file..."
## [1] "43 genes meeting criteria (FINAL_AD_MAG_00002)"
## [1] "0 genes meeting criteria (FINAL_AC_MAG_00008)"
## [1] "0 genes meeting criteria (FINAL_AC_MAG_00028)"
## [1] "0 genes meeting criteria (FINAL_AD_MAG_00021)"
## [1] "FINAL_AD_MAG_00021 has no associated variability file..."
## [1] "0 genes meeting criteria (FINAL_AD_MAG_00015)"
## [1] "1 genes meeting criteria (FINAL_AH_MAG_00007)"
## [1] "No codon information for any of ^these genes"
## [1] "169 genes meeting criteria (FINAL_AQ_MAG_00010)"
## [1] "0 genes meeting criteria (FINAL_AD_MAG_00004)"
## [1] "FINAL_AD_MAG_00004 has no associated variability file..."
plot_syn <- rbind(cbind(synUnlinked10, type="Bottom 10%"), cbind(synLinked90, type="Top 10%"),
cbind(synUnlinked01, type="Bottom 1%"), cbind(synLinked99, type='Top 1%'))
plot_syn <- plot_syn[!is.na(plot_syn$synonymity),]
plot_syn$gene <- factor(plot_syn$gene)
ggplot(plot_syn, aes(x=type, y = synonymity, fill=type, color=type)) +
geom_point(data=group_by(plot_syn, type) %>% sample_n(1000), size=0.4, alpha=0.5, position=position_jitter(width=0.3, height=0)) +
geom_violin(alpha=0.6) +
labs(fill="Percentile", x = "Percentile of gene pairwise correlations", y= "Synonymity") +
scale_fill_manual(values = brewer.pal(name='BrBG', n=4)[c(1,4,2,3)]) + scale_color_manual(values = brewer.pal(name='BrBG', n=4)[c(1,4,2,3)], guide=F) +
theme_basic
wilcox.test(x=plot_syn$synonymity[plot_syn$type == "Bottom 10%"], y=plot_syn$synonymity[plot_syn$type == "Top 10%"])
##
## Wilcoxon rank sum test with continuity correction
##
## data: plot_syn$synonymity[plot_syn$type == "Bottom 10%"] and plot_syn$synonymity[plot_syn$type == "Top 10%"]
## W = 3.6391e+11, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(x=plot_syn$synonymity[plot_syn$type == "Bottom 1%"], y=plot_syn$synonymity[plot_syn$type == "Top 1%"])
##
## Wilcoxon rank sum test with continuity correction
##
## data: plot_syn$synonymity[plot_syn$type == "Bottom 1%"] and plot_syn$synonymity[plot_syn$type == "Top 1%"]
## W = 4053762338, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
Some functions specific for Fig 4B:
fetchVariabilityCDNForMAG <- function(MAG="FINAL_AD_MAG_00009", focalGenes, minCov=10) {
codonVar <- read.csv(paste0("catalyst_idba/codon_var/", MAG, "-var-CDN.tsv"), sep="\t")
codonVar <- codonVar[codonVar$coverage >= minCov,]
codonVar <- codonVar[codonVar$corresponding_gene_call %in% focalGenes,]
codonVar[,grep("^...$",colnames(codonVar))] <- t(apply(codonVar[,grep("^...$",colnames(codonVar))], 1, function(x) x/sum(x)))
codonVar <- melt(codonVar, measure.vars = grep("^...$",colnames(codonVar)), variable.name = "codon", value.name = "codon_count")
codonVar$Day <- gsub("A[A-Z]","", as.character(codonVar$sample_id))
codonVar$Day <- factor(codonVar$Day, levels = mixedsort(unique(as.character(codonVar$Day)))) # sort human
codonVar$Participant <- gsub("(A[A-Z]).*$","\\1",as.character(codonVar$sample_id))
return(codonVar)
}
traceGenesByMAGwithSyn <- function(infoDF, covDF, binDF, bin, statsDF, focusGenes = "164483", codonPos="208", minGeneCov=0, tax=makegtdbtktaxMat(), proportional=T, facetGrid=T,
sqrt_transform=F, log10_transform=F) {
if (is.null(focusGenes)){
covDF <- covDF[covDF$key %in% infoDF[infoDF$bin == bin,"gene_callers_id"],] # get coverages for this MAG
print(paste0(nrow(covDF), " genes in ", bin))
} else {
covDF <- covDF[covDF$key %in% focusGenes,]
}
# get overall coverage for MAGs
binTrace <- setNames(as.numeric(binDF[binDF$bins == bin, grep("A", colnames(binDF))]), grep("A", colnames(binDF), value = T))
binTrace <- binTrace[mixedorder(names(binTrace))] # sort human (e.g. AB1 comes before AB15 before AB2)
codonVar <- fetchVariabilityCDNForGenes(genes = focusGenes, geneInfo = infoDF)
if (proportional){
codonVar[,grep("^...$",colnames(codonVar))] <- t(apply(codonVar[,grep("^...$",colnames(codonVar))], 1, function(x) x/sum(x)))
}
codonVar <- melt(codonVar, measure.vars = grep("^...$",colnames(codonVar)), variable.name = "codon", value.name = "codon_count")
# make long for ggplotting
covDF <- melt(covDF, id.vars = 'key')
colnames(covDF) <- c('Gene','Sample','Coverage')
# order things nicely and make participant
covDF$Gene <- factor(covDF$Gene)
covDF$Day <- gsub("A[A-Z]","", as.character(covDF$Sample))
covDF$Day <- factor(covDF$Day, levels = mixedsort(unique(as.character(covDF$Day)))) # sort human
covDF$Participant <- gsub("(A[A-Z]).*$","\\1",as.character(covDF$Sample))
codonVar$Day <- gsub("A[A-Z]","", as.character(codonVar$sample_id))
codonVar$Day <- factor(codonVar$Day, levels = mixedsort(unique(as.character(codonVar$Day)))) # sort human
codonVar$Participant <- gsub("(A[A-Z]).*$","\\1",as.character(codonVar$sample_id))
# grab genes above threshold
geneThresh <- ddply(covDF, ~ Gene + Participant, summarise, meanCov = mean(Coverage)) %>% filter(meanCov >= minGeneCov) %>% pull(Gene)
covDF <- covDF[covDF$Gene %in% unique(geneThresh),]
# make a little dataframe to facilitate overlaying the bin's coverage
binCov <- data.frame(Sample=factor(gsub("A[A-Z]","",names(binTrace)), levels = unique(gsub("A[A-Z]","",names(binTrace)))), Coverage=binTrace,
Participant = gsub("(A[A-Z]).*$","\\1",names(binTrace)))
codonVar <- codonVar[codonVar$codon_number %in% c(codonPos),]
codonVar <- codonVar[codonVar$codon_count > 0,]
covDF$codon_count <- covDF$Coverage
binCov$codon_count <- binCov$Coverage
# dummy category for y facet
covDF$category <- "higher"
binCov$category <- "higher"
codonVar$category <- "lower"
p <- ggplot(codonVar, aes(x = Day, y = codon_count)) +
facet_grid(category~Participant, scales='free', space='free_x') +
geom_line(data=covDF[covDF$Participant %in% codonVar$Participant,], aes(y=codon_count, x=Day, group=Gene), color='grey50', size = 0.8) +
geom_line(data = binCov[binCov$Participant %in% codonVar$Participant,], aes(group = Participant, x=Sample), size = 1, color = 'black') +
geom_col(aes(fill=codon), color='black') +
geom_label(data=codonVar %>% group_by(Day, Participant, synonymity, category) %>% summarise(codon_count=sum(codon_count)),
aes(label=round(synonymity,2), y=codon_count -0.25), size=2.5, label.padding = unit(0.15,'lines')) +
scale_y_continuous(expand = c(0,0)) + scale_x_discrete(expand = c(0,0.5)) + labs(x="Sample Day", y='Variant frequency', fill='Codon') + scale_fill_discrete(drop=T) +
labs(title=bin, subtitle = paste0(paste0(as.character(as.matrix(tax[tax$MAG == bin, c("Genus","Species")])), collapse = " "), # this is "Genus sp."
": ",round(statsDF$completion[statsDF$MAG == bin],1), "% C / ", round(statsDF$redundancy[statsDF$MAG == bin], 1), "% R, ",
round(statsDF$length[statsDF$MAG == bin]/1e6, 2), " Mb")) + guides(color=F) +
theme_basic + theme(strip.background = element_blank(), strip.text.y = element_blank(), strip.background.x = element_blank(), panel.spacing = unit(0.2, 'lines'))
print(p)
}
These functions were used to produce Figure 4B
traceGenesByMAGwithSyn(gene_info, gene_covs, bin_covs, bin="FINAL_AD_MAG_00009",statsDF = bin_stats, sqrt_transform = T, focusGenes = "164483", codonPos="208")
## `summarise()` regrouping output by 'Day', 'Participant', 'synonymity' (override with `.groups` argument)
And the two callouts in Figure 4C:
traceGenesByMAGwithSyn(gene_info, gene_covs, bin_covs, bin="FINAL_AE_MAG_00012",statsDF = bin_stats, sqrt_transform = T, focusGenes = "205006", codonPos="580")
## `summarise()` regrouping output by 'Day', 'Participant', 'synonymity' (override with `.groups` argument)
traceGenesByMAGwithSyn(gene_info, gene_covs, bin_covs, bin="FINAL_AQ_MAG_00003",statsDF = bin_stats, sqrt_transform = T, focusGenes = "581326", codonPos="59")
## `summarise()` regrouping output by 'Day', 'Participant', 'synonymity' (override with `.groups` argument)
Figure 4C:
magList <- unique(as.character(synUnlinked01$MAG))
for (mag in magList) {
df <- fetchVariabilityCDNForMAG(MAG=mag, focalGenes = unique(as.character(synUnlinked01$gene[synUnlinked01$MAG == mag])), minCov=10)
df <- df %>% group_by(unique_pos_identifier, corresponding_gene_call, Participant, codon) %>% filter(n()>=4 & any(codon_count > 0))
if (nrow(df) < 4) {
next()
}
df <- df %>% summarise(lmCoef = lm(codon_count~as.numeric(Day))$coefficients[2], meanSyn = mean(synonymity, na.rm = T))
df <- data.frame(df)
df$MAG <- mag
if (mag == magList[1]) {
codonTrend <- df
} else {
codonTrend <- rbind(codonTrend, data.frame(df))
}
}
codonTrend$quantile <- cut(codonTrend$lmCoef, c(-Inf, quantile(codonTrend$lmCoef, c(0.01,.99)), Inf), labels = c('2%','98%','2%'))
codonTrend$MAG <- factor(gsub("(^.*)(MAG_)(.*$)","\\2\\1\\3", gsub("_000","_", gsub("FINAL_","",as.character(codonTrend$MAG)))))
ggplot(codonTrend, aes(x=factor(corresponding_gene_call), y = lmCoef)) +
facet_grid(.~MAG, scales = 'free_x', space = 'free_x') +
geom_hline(yintercept = quantile(codonTrend$lmCoef, 0.99), lty=2) +
geom_hline(yintercept = quantile(codonTrend$lmCoef, 0.01), lty=2) +
geom_point(aes(color=quantile), position = position_jitter(height = 0, width = 0.2), alpha=0.4, pch=16) +
geom_violin(fill=NA, alpha=0.4, color='#00000066') +
labs(x="Genes", y="Slope of codon proportion over time in a mouth") +
scale_color_manual(values = c('darkred','steelblue3')) +
theme_classic() + theme(axis.text.x = element_blank())
Figure 4D:
ggplot(codonTrend, aes(x=lmCoef, y=meanSyn, color=quantile)) +
geom_point(alpha=0.4, pch=16) +
labs(x="Slope of codon proportion over time in a mouth", y = "Average synonymity") +
scale_color_manual(values = c('darkred','steelblue3')) +
theme_classic() + theme(panel.spacing = unit(0.2, 'line'))
