GenFlow: A Bioinformatic Workflow for PhyloGenomic Analysis

cat Data/short.txt
GCF_000633175.1
GCF_000820325.1
GCF_016805405.1
GCA_963892115.1
GCF_002850695.3
GCF_022631195.1
GCF_000014805.1
GCF_029025785.1
GCF_002285935.1
GCF_023920205.1
GCF_000820005.1
GCF_028355655.1

GenFlow -g short.txt -f INISA09F.fasta,INISA10F.fasta,INISA16F.fasta -t 8 -G 0.8 -F 0.8

#Download ncbi tool
curl -o datasets 'https://ftp.ncbi.nlm.nih.gov/pub/datasets/command-line/v2/linux-amd64/datasets'
chmod +x datasets

#Create a folder to save some logs
mkdir logs

#Download your references genomes and place them in a new folder named "Genomes"
./datasets download genome accession --inputfile short.txt >> logs/download.log 2>&1
unzip ncbi*
rm *.zip
mkdir Genomes
mv ncbi_dataset/data/GC*/*.fna Genomes/
rm -r ncbi_dataset md5sum.txt README.md datasets dataformat
cp *.fasta Genomes
cd Genomes

#Download esearch tool and set the path
sh -c "$(curl -fsSL https://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/install-edirect.sh)"
sh -c "$(wget -q https://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/install-edirect.sh -O -)" > ../logs/esearch.log 2>&1
echo "export PATH=\$HOME/edirect:\$PATH" >> $HOME/.bash_profile >> ../logs/esearch.log 2>&1
export PATH=${HOME}/edirect:${PATH} >> ../logs/esearch.log 2>&1

#Replace name with taxa
for f in GC* ; do term=$(echo $f | cut -f1,2 -d'_') ; esearch -db assembly -query $term | esummary | \
xtract -pattern DocumentSummary -sep ' ' -element Organism,Sub_value,AssemblyAccession | sed 's/ /_/g' | sed 's/://g' | \
sed 's+/+_+g' | sed 's/,/_/g' |  sed 's/[.]//g' | sed 's/-/_/g' | sed 's/_([^)(]*)//' | sed 's/=//g' ; done > NEW
ls -1 *.fna > OLD
paste OLD NEW|while read OLD NEW;do mv ${OLD} ${NEW};done
for i in ls *GCF*; do mv $i $i\.fasta; done
rm OLD NEW

#Reformat fasta files to be compatible with Anvio. New files will have the extension ".fa"
for i in `ls -1 *.fasta | sed 's/.fasta//'`
do
anvi-script-reformat-fasta $i\.fasta \
                           -o $i\.fa \
                           -l 1000 \
                           --simplify-names --seq-type NT >> ../logs/reformat.log 2>&1
    anvi-gen-contigs-database -f $i\.fa -o $i\.db >> ../logs/database.log 2>&1
    anvi-run-hmms -c $i\.db -T 32 >> ../logs/hmms.log 2>&1
done
#Some bash codes to create a requires file called to perform pangenomics (yes, we are doing that)
ls -1 *.db > path.txt
sed -i '1s/^/contigs_db_path\n/' path.txt
ls -1 *.db | sed 's/.db//' > name.txt
sed -i '1s/^/name\n/' name.txt
paste name.txt path.txt > external-genomes.txt
rm name.txt path.txt
#Save the genomes in a new Anvio database and perform pangenomics
anvi-gen-genomes-storage -e external-genomes.txt -o Filo-GENOMES.db >> \
../logs/storage.log 2>&1
anvi-pan-genome -g Filo-GENOMES.db --project-name Filo --num-threads 32 >> \
../logs/pangenome.log 2>&1
#Create the new file with all SCCGs with a geometric and functional index over 0.8
NG=$(ls *.fasta | wc -l)
anvi-get-sequences-for-gene-clusters -g Filo-GENOMES.db -p Filo/Filo-PAN.db \
-o core-sequences.fasta \
--max-num-genes-from-each-genome 1 --min-num-genomes-gene-cluster-occurs "${NG}" \
--concatenate-gene-clusters --min-geometric-homogeneity-index 0.8 \
--min-functional-homogeneity-index 0.8 >> ../logs/core.log 2>&1

mafft --retree 1 --thread 32 --maxiterate 0 core-sequences.fasta > core-sequences-aligned.fasta
mkdir ../results
anvi-gen-phylogenomic-tree -f core-sequences-aligned.fasta \
                           -o ../results/phylogenomic-tree.txt

average_nucleotide_identity.py -i fa -o ../results/pyANI --labels labels.txt --classes classes.txt -g --gmethod seaborn --gformat svg -v -l ba_ANI.log >> ../logs/pyANI.log

library(pheatmap)
library(viridis)
library(ggplot2)
library(viridis)
library(reshape2)
library(grid)
data <- read.table("ANIm_percentage_identity.tab", sep = "\t", header = T, row.names = 1)

# Convert the data to a matrix 
data_matrix <- as.matrix(data)

# Multiply all values by 100 for percentage
data_matrix <- data_matrix * 100

# Melt the data to long format
data_long <- melt(data_matrix)

# Custom formatting to remove decimal places for whole numbers
formatted_numbers <- apply(data_matrix, c(1, 2), function(x) {
  # Format numbers without decimal points
  if (x %% 1 == 0) {
    return(as.character(as.integer(x)))  # Convert to integer and then to character
  } else {
    return(as.character(round(x, 1)))  # Round to one decimal place for non-whole numbers
  }
})

# Create the heatmap
heatmap_plot <- pheatmap(
  data_matrix,
  cluster_rows = TRUE,           # Cluster the rows (to maintain order)
  cluster_cols = TRUE,           # Cluster the columns
  clustering_method = "complete", # Complete linkage clustering
  clustering_distance_rows = "euclidean", # Euclidean distance for rows
  clustering_distance_cols = "euclidean", # Euclidean distance for columns
  color = colorRampPalette(c("white", "orange", "red"))(100), # Gradient from white to orange to red
  display_numbers = formatted_numbers, # Use custom formatted numbers
  fontsize_number = 9,           # Increase font size of values to 8
  cellwidth = 20,                # Increase cell width for better spacing
  cellheight = 20,               # Increase cell height for better spacing
  scale = "none",                # No scaling of ANI values
  border_color = NA,             # Remove borders around cells
  angle_col = 90,                # Rotate x-axis labels 90 degrees
  legend = TRUE,                 # Title for the legend
  treeheight_row = 200,          # Increase row dendrogram height for larger branches
  treeheight_col = 200,          # Increase column dendrogram height for larger branches
  show_colnames = TRUE,          # Show x-axis labels
  show_rownames = TRUE,          # Show y-axis labels
  margin = c(10, 10, 10, 10),    # Increase margins (bottom, left, top, right)
  number_color = "black" 
)

# Save the heatmap to a PDF with appropriate dimensions
pdf("heatmap.pdf", width = 18, height = 18)  # Adjust size as needed

# Print the heatmap to the PDF device
print(heatmap_plot)

# Close the PDF device
dev.off()