This repository is made for eBook entitled “Bioinformatics Recipes for Plant Genomics: Data, Code, and Workflows” in Bio-101
This is an example workflow to construct gene tree and perform molecular dating using BEAST.
-
Running environment:
- The workflow was constructed based on Linux CentOS 7.7.1908 with OpenJDK Runtime Environment v1.8.0_242.
-
Required software and versions:
To download a pre-built version of GATK
cd program/
wget https://github.com/broadinstitute/gatk/releases/download/4.0.8.1/gatk-4.0.8.1.zip
unzip gatk-4.0.8.1.zip
cd gatk-4.0.8.1/
./gatk # check
cd ../../
To download and install BEAGLE and BEAST
# BEAGLE
cd lib/
wget https://github.com/beagle-dev/beagle-lib/archive/refs/tags/v3.1.1.tar.gz
tar -zxvf v3.1.1.tar.gz
cd beagle-lib-3.1.1/
./autogen.sh
./configure --prefix=$PWD
make install
ll lib # check
cd ../../
BEAGLE_PATH=$PWD/lib/beagle-lib-3.1.1/lib/ # save the absolute path for BEAST
# BEAST
cd program/
wget https://github.com/beast-dev/beast-mcmc/releases/download/v1.10.4/BEASTv1.10.4.tgz
tar -zxvf BEASTv1.10.4.tgz
cd BEASTv1.10.4/bin
export LD_LIBRARY_PATH=$BEAGLE_PATH:$LD_LIBRARY_PATH
./beast -beagle_info # check
cd ../../../
The input BAM files bam/*.bam
are compressed binary versions of SAM files that represent aligned sequences of the samples to the reference genome.
BAM/SAM is commonly used to represent sequence alignment of next-generation sequencing data. Specification of the format can be found here.
You can generate BAM files for your samples with most alignment tools such as bwa or bowtie.
We will start from these files and generate site-by-site GVCF files for making gene sequence files.
The file input/sample_info.txt
a tab-delimited file containing the names and group assignments of the samples. The group assignment will be used for plotting.
The first ten rows of the file:
ANN0955 Ann-0
ANN0803 Ann-0
ANN0979 Ann-0
HT170_SD2W-14 Ann-0
ANN1272 Ann-0
HT219_SD2W-38 Ann-2
ANN1374 Ann-2
ANN1146 Ann-2
ANN1159 Ann-2
ANN0870 Ann-2
The first 5 samples are from species Ann
with genotype 0/0 of the focal chromosomal inversion (group 0), and the next 5 samples are from species Ann
with genotype 1/1 (group 2). We also include 2 samples from each of other related (sub)species and the outgroup.
The file input/genes.locations.txt
is a tab-delimited file with the following columns: gene name (strings without spaces), chromosome, start and end (in bp, based on the reference genome). It will be used to extract gene sequences from BAM files.
The first five rows of the example data:
gene_1 HanXRQChr17 27130172 27133822
gene_2 HanXRQChr17 191345873 191353428
gene_3 HanXRQChr17 191494056 191496274
gene_4 HanXRQChr17 193856630 193871419
gene_5 HanXRQChr17 193971597 193974985
gatk="program/gatk-4.0.8.1/gatk"
ref="input/HanXRQr1.0-20151230.fa"
BAM_DIR="input/bam"
# Index the reference for GATK
samtools faidx $ref
$gatk CreateSequenceDictionary --java-options -Xmx2G -R $ref
# Index the BAM files for GATK
for bam in `ls $BAM_DIR/*.bam`
do
samtools index $bam
done
ll input/bam/*.bam.bai
# Generate gene FASTAs from BAMs
while read gene_name chr start end
do
if [ ! -s "cache/$gene_name.fasta" ]
then
echo -e "$chr\t$start\t$end" > cache/tmp.$gene_name.bed
cut -f1 input/sample_info.txt | parallel -j 10 bash script/make_fasta.sh $BAM_DIR $ref $gene_name
for f in cache/tmp.$gene_name.*.fa; do (cat "${f}"; echo); done > cache/$gene_name.fasta
rm cache/tmp.*.bed cache/tmp.*.*.g.vcf* cache/tmp.*.*.fa
fi
done < input/genes.locations.txt
This step first produces a g.vcf
for each gene for each sample from the BAM file using gatk HaplotypeCaller
. The parameter -ERC BP_RESOLUTION
outputs the genotype calls site by site, which will ease the conversion of g.vcf
to FASTA format using the custom script script/gvcf2fasta_nogaps.pl
.
A FASTA will be generated in cache/
for each gene with sequences of all samples. Check if the numbers of sequences are correct:
ls cache/gene_*.fasta | wc -l # 82 genes in total
wc -l cache/gene_*.fasta # should be 40 (20 sequences in each FASTA) for all genes
We filter the genes based on individual and average missing rates across samples and generate a gene list for phylogenetic analysis
cut -f1 input/genes.locations.txt | perl script/filter_phylo_genes.pl > cache/phylo_genes.list
Check the numbers of genes that pass the filter:
wc -l cache/phylo_genes.list # 40 out of 82 genes kept
# FASTA list
cat cache/phylo_genes.list | sed 's/$/&.fasta/g' | sed 's/^/cache\//g' > cache/gene_fastas.list
# Generate the XML file for BEAST v1 using custom script
perl script/fastas2BEAST1xml.pl \
--fastaList cache/gene_fastas.list \
--chainLength 1000000 \
--screenEcho 1000 \
--fileLog 100 \
--out output/out
Here we set the number of steps of MCMC (--chainLength
) to be 1M and set the parameter values to be displayed on the screen (--screenEcho
) every 1k steps and recorded in the log file (--fileLog
) every 100 steps.
This will generate output/out.xml
, which BEAST will use to run the analysis.
BEAGLE_PATH=$PWD/lib/beagle-lib-3.1.1/lib/ # specify the path to the BEAGLE library
export LD_LIBRARY_PATH=$BEAGLE_PATH:$LD_LIBRARY_PATH
program/BEASTv1.10.4/bin/beast -threads 30 output/out.xml
This is the main command of BEAST. It takes 1-2h with 30 threads. This step will generate output/out.log
, output/out.trees
and output/out.ops
.
We assume that the phylogenetic MCMC in BEAST reach convergence after 10k runs. We then create a maximum clade credibilty (MCC) tree summarizing the posterior tree distribution using the program TreeAnnotator within the BEAST package.
program/BEASTv1.10.4/bin/treeannotator -burnin 100000 -heights median output/out.trees output/out.mcc.tre
This will generate the file output/out.mcc.tre
that contains a consensus tree as well as other information of phylogenetic analysis in BEAST.
Rscript script/plot_beast.R
It is a free and open source software, licensed under GPLv3.