Note

Following includes codes used to generate an 79 way cactus alignment of birds. It uses bash and R codes interchangeably along with chunks that should be submitted as jobs using a resource manager in a supercomputing cluster (here we use slurm and snakemake). Please make sure you are running the codes in the correct interface.

1 Preprocessing

We selected 78 bird genomes and 1 genome for outgroup and aligned them using progressive cactus. The selection of taxa follows birds having short tarsi that were identified using bayou. See the file bayou.Rmd file for steps on running the bayou analysis. The 79 genomes and their relevant information is shown in the following table.

Order Family Latin_name Genome_Accession Sex Scaffold_Num Scaffold_N50_bp Scaffold_L50 Total_Len_bp File_ID Notes
Crocodylia Gavialidae Gavialis_gangeticus GCA_001723915.1 Female 81 96076944 8 2640792433 Gavialis_gangeticus_GCA001723915.fa.gz
Struthioniformes Struthionidae Struthio_camelus ASM69896v1_HiC Unknown 5981 83444685 5 1225496327 Struthio_camelus_ASM69896.fa.gz https://www.dropbox.com/s/1g03spa7xo9aaq9/ASM69896v1_HiC.fasta.gz?dl=0
Rheiformes Rheidae Rhea_americana rheAme1 Unknown 1395 80688299 4 1161242837 Rhea_americana_rheAme1.fa.gz https://www.dropbox.com/s/c52npppt0m7e2dv/rheAme1_HiC.fasta.gz?dl=0
Tinamiformes Tinamidae Eudromia_elegans GCA013398775.1 Female 2603 6275819 50 958982151 Eudromia_elegans_GCA013398775.fa.gz
Casuariiformes Casuariidae Dromaius_novaehollandiae droNov1 Unknown 1863 82460724 5 1192718721 Dromaius_novaehollandiae_droNov1.fa.gz https://www.dropbox.com/s/bge457t6l9a5uqp/droNov1_HiC.fasta.gz?dl=0
Anseriformes Anatidae Anas_platyrhynchos GCA_000355885.1 Female 1665 77626585 5 1211992756 Anas_platyrhynchos_GCA900411745.fa.gz
Galliformes Phasianidae Gallus_gallus4 GCA_000002315.1 Female 16846 12877381 23 1046915324 Gallus_gallus_GCA000002315.fa.gz
Galliformes Phasianidae Gallus_gallus7 GCF_016699485.2 Female 214 90861225 4 1053332251 Gallus_gallus_GCF016699485.fa.gz
Phoenicopteriformes Phoenicopteridae Phoenicopterus_ruber GCA_009819775.1 Female 90 85472901 5 1246226932 Phoenicopterus_ruber_GCA009819775.fa.gz
Columbiformes Columbidae Columba_livia GCA_001887795.1 Male 91 24540102 11 1018016946 Columba_livia_GCA001887795.fa.gz
Musophagiformes Musophagidae Tauraco_erythrolophus GCA_009769465.1 Female 140 85562082 6 1252271553 Tauraco_erythrolophus_GCA009769465.fa.gz
Cuculiformes Cuculidae Cuculus_canorus GCA_017976375.1 Female 151 78325852 5 1180156273 Cuculus_canorus_GCA017976375.fa.gz
Caprimulgiformes Caprimulgidae Caprimulgus_europaeus GCA_907165065.1 Female 121 82614289 6 1177791212 Caprimulgus_europaeus_GCA907165065.fa.gz
Apodiformes Apodidae Apus_apus GCA_020740795.1 Female 109 115508478 4 1100361224 Apus_apus_GCA020740795.fa.gz
Apodiformes Trochilidae Calypte_anna GCA_003957555.2 Female 159 74081004 4 1059687259 Calypte_anna_GCA003957555.fa.gz
Gruiformes Rallidae Porphyrio_hochstetteri GCA_020800305.1 Female 174 71566193 5 1270294353 Porphyrio_hochstetteri_GCA020800305.fa.gz
Charadriiformes Charadriidae Pluvialis_apricaria GCA_017639485.1 Female 107 87527607 5 1247767512 Pluvialis_apricaria_GCA017639485.fa.gz
Charadriiformes Laridae Sterna_hirundo GCA_009819605.1 Female 123 85487239 5 1229972541 Sterna_hirundo_GCA009819605.fa.gz
Sphenisciformes Spheniscidae Aptenodytes_patagonicus GCA_010087175.1 Unknown 15968 2912628 122 1252291219 Aptenodytes_patagonicus_GCA010087175.fa.gz
Sphenisciformes Spheniscidae Aptenodytes_forsteri GCA_000699145.1 Male 29970 5071598 71 1257483768 Aptenodytes_forsteri_GCA000699145.fa.gz
Sphenisciformes Spheniscidae Pygoscelis_adeliae GCA_000699105.1 Male 90746 5047175 73 1226103150 Pygoscelis_adeliae_GCA000699105.fa.gz
Sphenisciformes Spheniscidae Pygoscelis_papua GCA_010090195.1 Unknown 18750 2789995 121 1306937499 Pygoscelis_papua_GCA010090195.fa.gz
Sphenisciformes Spheniscidae Pygoscelis_antarcticus GCA_010078415.1 Unknown 16691 6271244 57 1260447050 Pygoscelis_antarcticus_GCA010078415.fa.gz
Sphenisciformes Spheniscidae Megadyptes_antipodes GCA_010078485.1 Unknown 17379 23315117 17 1314943623 Megadyptes_antipodes_GCA010078485.fa.gz
Sphenisciformes Spheniscidae Eudyptula_novaehollandiae GCA_010078495.1 Unknown 15946 29280209 16 1353666924 Eudyptula_novaehollandiae_GCA010078495.fa.gz
Sphenisciformes Spheniscidae Eudyptula_albosignata GCA_010080465.1 Unknown 18591 21866543 17 1370511393 Eudyptula_albosignata_GCA010080465.fa.gz
Sphenisciformes Spheniscidae Eudyptula_minor GCA_010080355.1 Unknown 22396 21127646 19 1462615782 Eudyptula_minor_GCA010080355.fa.gz
Sphenisciformes Spheniscidae Spheniscus_demersus GCA_010077935.1 Unknown 15338 15386364 25 1275043525 Spheniscus_demersus_GCA010077935.fa.gz
Sphenisciformes Spheniscidae Spheniscus_humboldti GCA_010076325.1 Unknown 11265 6229819 66 1239252862 Spheniscus_humboldti_GCA010076325.fa.gz
Sphenisciformes Spheniscidae Spheniscus_magellanicus GCA_010076225.1 Unknown 13422 12679469 33 1260229131 Spheniscus_magellanicus_GCA010076225.fa.gz
Sphenisciformes Spheniscidae Eudyptes_pachyrhynchus GCA_010085335.1 Unknown 52188 8912950 39 1295505775 Eudyptes_pachyrhynchus_GCA010085335.fa.gz
Sphenisciformes Spheniscidae Eudyptes_sclateri GCA_010078445.1 Unknown 3484 1921244 174 1211402317 Eudyptes_sclateri_GCA010078445.fa.gz
Sphenisciformes Spheniscidae Eudyptes_chrysolophus GCA_010084205.1 Unknown 15753 13794837 29 1365695164 Eudyptes_chrysolophus_GCA010084205.fa.gz
Sphenisciformes Spheniscidae Eudyptes_schlegeli GCA_010080425.1 Unknown 17083 1894246 198 1308180354 Eudyptes_schlegeli_GCA010080425.fa.gz
Sphenisciformes Spheniscidae Eudyptes_chrysocome GCA_010085355.1 Unknown 3371 1949323 169 1231062710 Eudyptes_chrysocome_GCA010085355.fa.gz
Sphenisciformes Spheniscidae Eudyptes_filholi GCA_010085365.1 Unknown 2312 6429221 29 1223976468 Eudyptes_filholi_GCA010085365.fa.gz
Sphenisciformes Spheniscidae Eudyptes_moseleyi GCA_010082375.1 Unknown 20381 2252023 158 1301647863 Eudyptes_moseleyi_GCA010082375.fa.gz
Procellariiformes Procellariidae Macronectes_giganteus Mgig Male 861 27376298 14 1247958164 Macronectes_giganteus_Mgig.fa.gz http://kusglab.org/downloads/Mgig_genome_final.fasta.gz
Pelecaniformes Threskiornithidae Theristicus_caerulescens GCA_020745775.1 Female 88 101368359 5 1203372763 Theristicus_caerulescens_GCA020745775.fa.gz
Accipitriformes Accipitridae Aquila_chrysaetos GCA_900496995.4 Female 144 46934974 9 1233704830 Aquila_chrysaetos_GCA900496995.fa.gz
Strigiformes Strigidae Athene_cunicularia GCA_003259725.1 Male 445 42147404 6 1157069330 Athene_cunicularia_GCA003259725.fa.gz
Trogoniformes Trogonidae Trogon_surrucura GCA_020746105.1 Female 182 82617956 5 1165749515 Trogon_surrucura_GCA020746105.fa.gz
Bucerotiformes Bucerotidae Bucorvus_abyssinicus GCA_009769605.1 Female 112 40274051 9 1132597561 Bucorvus_abyssinicus_GCA009769605.fa.gz
Coraciiformes Alcedinidae Halcyon_senegalensis GCA013397595.1 Male 3682 2092632 146 1137737933 Halcyon_senegalensis_GCA013397595.fa.gz
Coraciiformes Alcedinidae Chloroceryle_aenea GCA013399075.1 Female 67943 323453 918 1067751506 Chloroceryle_aenea_GCA013399075.fa.gz
Coraciiformes Alcedinidae Actenoides_hombroni GCA_024336885.1 Unknown 10367 10995457 14 1152258273 Actenoides_hombroni_GCA024336885.fa.gz
Coraciiformes Meropidae Merops_nubicus GCA_009819595.1 Female 109 49132790 9 1149356202 Merops_nubicus_GCA009819595.fa.gz
Piciformes Lybiidae Pogoniulus_pusillus GCA_015220805.1 Female 253 33994490 13 1272358903 Pogoniulus_pusillus_GCA015220805.fa.gz
Piciformes Picidae Colaptes_auratus GCA_015227895.2 Unknown 2369 43947546 11 1378187165 Colaptes_auratus_GCA015227895.fa.gz
Falconiformes Falconidae Falco_peregrinus GCA_024432695.1 Female 642 26062782 14 1153798512 Falco_peregrinus_GCA024432695.fa.gz Replacing GCA023634155
Psittaciformes Psittaculidae Melopsittacus_undulatus GCA_012275295.1 Female 864 104092594 5 1171617451 Melopsittacus_undulatus_GCA012275295.fa.gz
Passeriformes Acanthisittidae Acanthisitta_chloris GCA_016904835.1 Female 206 40827287 8 1077387582 Acanthisitta_chloris_GCA016904835.fa.gz
Passeriformes Tyrannidae Myiozetetes_cayanensis GCA_022539395.1 Female 1692 63679618 6 1075628138 Myiozetetes_cayanensis_GCA022539395.fa.gz
Passeriformes Pipridae Chiroxiphia_lanceolata GCA_009829145.1 Female 93 75608808 5 1089631598 Chiroxiphia_lanceolata_GCA009829145.fa.gz
Passeriformes Maluridae Malurus_cyaneus GCA_009741485.1 Female 4330 68113514 6 1078891882 Malurus_cyaneus_GCA009741485.fa.gz
Passeriformes Meliphagidae Lichenostomus_cassidix GCA_008360975.2 Female 894 63800663 6 1102944376 Lichenostomus_cassidix_GCA008360975.fa.gz
Passeriformes Laniidae Lanius_collurio GCA_020086605.1 Male 39581 23147665 16 1205288378 Lanius_collurio_GCA020086605.fa.gz
Passeriformes Corvidae Corvus_cornix GCA_000738735.5 Male 47 73732264 5 1030881721 Corvus_cornix_GCA000738735.fa.gz
Passeriformes Alaudidae Eremophila_alpestris GCA_009792885.1 Male 2714 9421811 30 1041035471 Eremophila_alpestris_GCA009792885.fa.gz
Passeriformes Paridae Parus_major GCA_001522545.2 Male 1675 71365269 5 1020310768 Parus_major_GCA001522545.fa.gz
Passeriformes Acrocephalidae Acrocephalus_scirpaceus GCA_910950805.1 Unknown 200 74438198 5 1075083815 Acrocephalus_scirpaceus_GCA910950805.fa.gz
Passeriformes Hirundinidae Progne_subis GCA_022316685.1 Unknown 2896 6129949 44 1165951862 Progne_subis_GCA022316685.fa.gz
Passeriformes Hirundinidae Riparia_riparia GCA_020917445.1 Male 19999 13690488 28 1093155110 Riparia_riparia_GCA020917445.fa.gz
Passeriformes Hirundinidae Hirundo_rustica GCA_015227805.3 Female 617 76187387 5 1105955550 Hirundo_rustica_GCA015227805.fa.gz
Passeriformes Pycnonotidae Brachypodius_melanocephalos GCA01339861.1 Female 3820 17537833 14 906727256 Brachypodius_melanocephalos_GCA013398615.fa.gz added_62993985_bp_of_potential_Z_but_no_resequencing_data_for_that_region
Passeriformes Pycnonotidae Pycnonotus_jocosus GCA013400435.1 Female 142341 214665 1144 1035289814 Pycnonotus_jocosus_GCA013400435.fa.gz
Passeriformes Phylloscopidae Phylloscopus_trochilus GCA_016584745.1 Male 39596 12447930 24 1168570690 Phylloscopus_trochilus_GCA016584745.fa.gz
Passeriformes Sylviidae Sylvia_atricapilla GCA_009819655.1 Unknown 189 72978054 5 1066786587 Sylvia_atricapilla_GCA009819655.fa.gz
Passeriformes Certhiidae Certhia_americana GCA_018697195.1 Unknown 8309 64356598 6 1113356413 Certhia_americana_GCA018697195.fa.gz
Passeriformes Sturnidae Sturnus_vulgaris GCA_023376015.1 Female 1344 72244370 5 1043825371 Sturnus_vulgaris_GCA023376015.fa.gz
Passeriformes Muscicapidae Ficedula_albicollis GCA_000247815.2 Male 21428 64724594 6 1118343587 Ficedula_albicollis_GCA000247815.fa.gz
Passeriformes Estrildidae Taeniopygia_guttata GCA_003957565.4 Male 199 70982421 6 1056271262 Taeniopygia_guttata_GCA003957565.fa.gz
Passeriformes Passeridae Passer_domesticus GCA_001700915.1 Female 2301 68733010 6 1042720703 Passer_domesticus_GCA001700915.fa.gz
Passeriformes Motacillidae Motacilla_alba GCA_015832195.1 Male 619 72386170 6 1072670728 Motacilla_alba_GCA015832195.fa.gz
Passeriformes Fringillidae Serinus_canaria GCA_022539315.1 Female 684 72083177 5 1050336073 Serinus_canaria_GCA022539315.fa.gz
Passeriformes Emberizidae Emberiza_elegans GCA_022818055.1 Male 380 73992973 6 1222647951 Emberiza_elegans_GCA022818055.fa.gz
Passeriformes Icteridae Molothrus_ater GCA_012460135.1 Female 447 52124711 7 1087312585 Molothrus_ater_GCA012460135.fa.gz
Passeriformes Parulidae Setophaga_coronata GCA_001746935.2 Male 2109 71386623 5 1022651334 Setophaga_coronata_GCA001746935.fa.gz
Passeriformes Thraupidae Diglossa_brunneiventris GCA_019023105.1 Female 600 67281049 6 1079473785 Diglossa_brunneiventris_GCA019023105.fa.gz

1.1 Preparing genome data

Get the appropriate files from NCBI or other sources (provided in Notes column). Rename all files as listed in the File_ID column of the table.

Important: Gunzip all the files, the cactus pipeline does note work with gzipped or bgzipped files.

Clean the headers of the scaffold names. We got rid of any extra elements in the scaffold line except for the NCBI Accession Number. This can be easily accomplished by using the following command. Replace the genome.fa with the name of the genomes. You can loop over all the files to run the command on all the fasta files

#bash
for FILE in ./*.fa; do
  sed -i 's/\(>[A-Z]*[0-9]*\.[0-9]\).*/\1/' FILE
done

Use samtools faidx to generate index files for each genome. Small scaffolds may cause issues with aligning using cactus. So we remove all scaffolds less than 1000 bp.

#R
#load all fai index files
index_files <- list.files("./", pattern = "*.fa.gz.fai", full.names = T)
for (file in index_files){
  file_open <- read.table(file, sep="\t", header = F)
  species <- gsub(".fa.gz.fai", "", basename(file), fixed = T)
  file_keep_scafs <- subset(file_open, V2 >= 1000)
  #samtools based format scaffold:start-end
  file_keep_scafs$V6 <- paste(file_keep_scafs$V1, ":1-", file_keep_scafs$V2, sep = "")
  write.table(file_keep_scafs$V6, paste(species,".txt", sep = ""), 
              quote = F, row.names = F, col.names = F)
}
#bash
for FILE in ./*.txt; do
  species=$(basename $FILE .txt)
  genomes=$(basename $(find . -type f -name $species*.fa -print) .fa)
  echo $genomes
  samtools faidx ${genomes}.fa -r ${species}.txt > ${genomes}_1000.fa
  mv ${genomes}_1000.fa ${genomes}.fa
done

1.2 Note: Genome soft masking

The cactus snakemake pipeline will get stuck if genomes are not masked. Most genomes from NCBI are already masked but still some genomes might give you issues. The snakemake cactus pipeline does initially mask the genome but even that does not work sometimes. I would recommend running RepeatModeler and RepeatMasker on genomes that fail or stall. We had to do this for the Macronectes giganteus genome.

#bash
source activate repeatmodeler

BuildDatabase -name Mgig -engine ncbi Macronectes_giganteus_Mgig.fa
RepeatModeler -engine ncbi -pa 20 -database Mgig
RepeatMasker -species Aves -parallel 20 -xsmall -s Macronectes_giganteus_Mgig.fa
RepeatMasker -lib Mgig-families.fa -parallel 20 -xsmall -s Macronectes_giganteus_Mgig.fa.masked

1.3 Getting the tree

The cactus snakemake pipeline needs a tree to run. This tree need not be the correct species tree, just needs to be a guide tree. We used the Genetic data only trees (6670 species), Hackett backbone Jetz tree for this purpose (download here: https://birdtree.org/downloads/)

We subset the requisite tips and updated the names of the taxa to match the taxon names in our samples.

library(ape)
library(phytools)
big_tree <- read.tree("b10k.tree")
index_table <- readxl::read_xlsx("genome_list.xlsx")

#Listed taxa have different names in tree versus our sample sheet
species_all <- c(intersect(big_tree$tip.label, index_table$Latin_name), "Lichenostomus_melanops", 
                 "Corvus_corone", "Brachypodius_atriceps", "Emberiza_buchanani", "Gallus_gallus")

#Eudyptula minor has been split into four species for which we have genomes
#We created a polytomy at node with Eudpytula minor to include the four other species
all_tree <- keep.tip(big_tree, species_all)
all_tree <- bind.tip(all_tree, "Eudyptula_novaehollandiae", where = which(all_tree$tip.label=="Eudyptula_minor"))
all_tree <- bind.tip(all_tree, "Eudyptula_albosignata", where = which(all_tree$tip.label=="Eudyptula_minor"))
all_tree <- bind.tip(all_tree, "Eudyptes_schlegeli", where = which(all_tree$tip.label=="Eudyptes_chrysolophus"))

#Eudyptula chrysocome has been split into two species for which we have genomes
#We created a polytomy at node with Eudpytula chryscome to include the two other species
all_tree <- bind.tip(all_tree, "Eudyptes_filholi", where = which(all_tree$tip.label=="Eudyptes_chrysocome"))

#The following species have had their names changed or been split but we only need one taxa
#These nodes were renamed
all_tree$tip.label[which(all_tree$tip.label == "Lichenostomus_melanops")] <- "Lichenostomus_cassidix"
all_tree$tip.label[which(all_tree$tip.label == "Corvus_corone")] <- "Corvus_cornix"
all_tree$tip.label[which(all_tree$tip.label == "Brachypodius_atriceps")] <- "Brachypodius_melanocephalos"
all_tree$tip.label[which(all_tree$tip.label == "Emberiza_buchanani")] <- "Emberiza_elegans"

#Since we included two gallus gallus genomes in the cactus file, we need to duplicate
#and rename that node
all_tree$tip.label[which(all_tree$tip.label == "Gallus_gallus")] <- "Gallus_gallus4"
all_tree <- bind.tip(all_tree, "Gallus_gallus7", where = which(all_tree$tip.label=="Gallus_gallus4"))

#Add an outgroup node for the gharial genome
all_tree <- bind.tip(all_tree, "Gavialis_gangeticus", where = "root")

#plot(all_tree, root.edge = T)
write.tree(all_tree, file="all_tree.tree")
library(ape)
library(ggplot2)
library(ggtree)
library(ggtreeExtra)
tree <- read.tree("all_tree.tree")
index_table <- readxl::read_xlsx("genome_list.xlsx")
#plot(tree)

p <- ggtree(tree, layout="circular", size=0.1)+
  geom_tiplab(size=4)+
  geom_fruit(data=index_table, 
             geom=geom_bar, 
             mapping = aes(y=Latin_name, x=Scaffold_N50_bp), 
             stat="identity",
             width=0.5,
             pwidth = 0.3,
             orientation="y",
             offset=1.2,
             grid.params = list(vline = TRUE),
             axis.params = list(axis="x",
                               text.size=1))
p

1.4 Creating the cactus map file

This file has a simple format, the first line includes the newick tree. The subsequent lines include two columns, first column will include the name of the node (aka species name) and the second column will include the location of the file.

2 Running snakemake

Check this website for scripts and instruction to run this snakemake code:
https://github.com/gwct/GenomeAnnotation/tree/master/genome_alignment

The following code was run on Harvard Supercomputing Cluster and ran for ~90 days.

#!/bin/bash
#SBATCH -c 1
#SBATCH -t 99-00:00
#SBATCH -p holy-smokes
#SBATCH --nodelist=holy7c0908
#SBATCH --mem=5000
#SBATCH -o ./logs/cac_%j.out
#SBATCH -e ./logs/cac_%j.err

cd /n/holyscratch01/informatics/sshakya/genomes_cactus/

source activate PhyloAcc

singularity exec --cleanenv ~/cactus_v2.1.1-gpu.sif cactus-prepare /n/holyscratch01/informatics/sshakya/genomes_cactus/cactus_snakemake/new_bird_aln.txt --outDir /n/holyscratch01/informatics/sshakya/genomes_cactus/cactus_snakemake/snakemake_out --jobStore /n/holyscratch01/informatics/sshakya/genomes_cactus/cactus_snakemake/snakemake_temp --gpu

snakemake -p -s ~/GenomeAnnotation/genome_alignment/cactus-gpu-snakemake/cactus_gpu.smk --configfile /n/holyscratch01/informatics/sshakya/genomes_cactus/all_birds_config.yaml --profile ~/GenomeAnnotation/genome_alignment/cactus-gpu-snakemake/profiles/slurm_profile/