Bioinformatic pipeline for “Genome skimming resolves Neotropical population structure and mito-nuclear discordance in the common vampire bat (Desmodus rotundus)” by Lexi E. Frank, Marissa Milstein, Christopher Shaffer, Philip Suse, Steven Suse, Achaicus Chekema, Evan J. Kipp, Roxanne J. Larsen, Tiffany Wolf, Peter A. Larsen and Christopher Faulk.
Bioinformatic methods adapted from “A beginner’s guide to low-coverage whole genome sequencing for population genomics” tutorial and ANGSD guides.
DNA was extracted from Desmodus rotundus kidney and liver tissues with the Qiagen MagAttract Kit. Nanopore library prep was performed following manufacturer’s instructions with the SQK-NBD114.24 kit. Sequencing was performed on a P2Solo to a coverage <1X.
Combine all bams generated per barcode into one file using samtools cat. For this example, barcode01 refers to one sample/individual.
# concatenate bams
samtools cat -o allbarcode01.bam barcode01*.bamMap bams to reference genome with Dorado aligner. Dorado aligner utilizes minimap2. Index files with Samtools.
# map reads to reference
dorado aligner reference.fna allbarcode01.bam | samtools sort -@ 32 -o barcode01_mapped.bam
# index the bam
samtools index barcode01_mapped.bamMake a reference genome index using Samtools faidx.
# make reference index
samtools faidx reference.fnaGenerate summary stats using Cramino like number of mapped reads, mean coverage, read N50, etc.
# generate summary stats
cramino -t 32 allbarcode01_mapped.bam > allbarcode01_cramino.txtUse Mitohifi with singularity to assemble mitchondrial genomes guided by a reference.
# pull mitohifi
sudo docker pull ghcr.io/marcelauliano/mitohifi:master
# install singularity
wget https://github.com/sylabs/singularity/releases/download/v4.3.4/singularity-ce_4.3.4-jammy_amd64.deb
# run mitohifi on singularity
singularity shell --bind /home/larsenlab/Desktop/desmodus_skim_seq/MitoHiFi/:/MitoHiFi docker://ghcr.io/marcelauliano/mitohifi:master mitohifi.py -h
cd MitoHiFi
# find the correct reference
findMitoReference.py --species "Desmodus rotundus" --outfolder . --min_length 14000
# make fastq from Bam to use as input into mitohifi
samtools fastq allbarcode01_mapped.bam > allbarcode01.fq
# run mitohifi on the fastq file, the file final_mitogenome.fasta is the final consensus sequence
# -o 2 for vertebrate mitochondrial code for annotation
#- f is the reference fasta downloaded above, -g is the reference genbank file downloaded above
mitohifi.py -r allbarcode01.fq -f NC_022423.1.fasta -g NC_022423.1.gb -t 32 -o 2Use Flye to de novo assembly, if mitohifi does not assemble a mitogenome.
# assembly
flye --nano-corr --meta allbarcode01.fastq --genome-size 16k --out-dir Assembly_barcode01Generate phylogenetic trees from mitogenome data using Mafft aligner and iqtree tree generator. Visualize the tree in Figtree.
# put all fasta mitogenomes together into one file with references and/or outgroups.
cat *.fasta > all_mito.fasta
# Align with MAFFT
# conda install bioconda::mafft
mafft --auto --adjustdirection all_mito.fasta > all_mito.aln
# Make tree with iqtree
iqtree2 -s all_mito.aln -nt 24 -B 1000
#visualize in FigtreeMakes consensus fasta for COI and CYTB using samtools, minimap2, and medaka consensus.
# Convert bam to fastq
samtools fastq allbarcode01.bam > barcode01.fastq
# Make COI consensus
# Map reads, COI_ref.fna is reference file of coi sequences
minimap2 -t 4 -ax map-ont -g 1000 COI_ref.fasta barcode01.fastq > barcode01_COI_unsorted.sam
# Sort, filter, and extract mapped reads
samtools view -@ 4 -bS barcode01_COI_unsorted.sam | samtools sort -o barcode01_COI_sorted.bam
samtools view -@ 4 -b -F 4 barcode01_COI_sorted.bam | samtools fastq - > barcode01_COI_extracted.fastq
# Run Medaka consensus
medaka_consensus -i barcode01_COI_extracted.fastq -d COI_ref.fasta -o medaka_COI_out -m "r1041_e82_400bps_sup_v4.3.0" -t 4
# Make CYTB consensus
# Map reads, cytb_ref.fasta is a reference file of cytb sequences
minimap2 -t 4 -ax map-ont -g 1000 cytb_ref.fasta barcode01.fastq > barcode01_cytb_unsorted.sam
# Sort, filter, and extract mapped reads
samtools view -@ 4 -bS barcode01_cytb_unsorted.sam | samtools sort -o barcode01_cytb_sorted.bam
samtools view -@ 4 -b -F 4 barcode01_cytb_sorted.bam | samtools fastq - > barcode01_cytb_extracted.fastq
# Run Medaka consensus
medaka_consensus -i barcode01_cytb_extracted.fastq -d cytb_ref.fasta -o medaka_cytb_out -m "r1041_e82_400bps_sup_v4.3.0" -t 4Generate phylogenetic trees from COI/CYTB data using Mafft aligner and iqtree tree generator. Visualize the tree in Figtree.
# Combine all individual COI fasta files into one unaligned file plus outgroup, then align with MAFFT
cat *_COI.fasta > allCOI.fasta
mafft --auto --adjustdirection allCOI.fasta > aligned_COI.fasta
# Combine all individual cytb fasta files into one unaligned file, then align with MAFFT
cat *_cytb.fasta > allcytb.fasta
mafft --auto --adjustdirection allcytb.fasta > aligned_cytb.fasta
# Run the COI tree
iqtree2 -s aligned_COI_qual.fasta -m MFP -B 1000 -T AUTO
# Run the cytb tree
iqtree2 -s aligned_cytb_qual.fasta -m MFP -B 1000 -T AUTO
#visualize in figtreeANGSD is for analyzing next generation sequencing data. Use genotype likelihoods to take uncertainty into account instead of basing the analysis on called genotypes for low coverage data.
#download htslib
git clone --recurse-submodules https://github.com/samtools/htslib.git;
#download ANGSD
git clone https://github.com/angsd/angsd.git;
#install htslib
cd htslib
make
#install ANGSD
cd ../angsd
make HTSSRC=../htslibMake a list of the bam files to be analyzed together.
# make list of bam files
ls <path to bams>/*.bam > bam.filelist# -GL 1 is for the Samtools model, which is potentially better for low coverage data and faster
./angsd -GL 1 -out genolike -nThreads 20 -doGlf 2 -doMajorMinor 1 -SNP_pval 1e-6 -doMaf 1 -minInd [minimum number of individuals needing at a loci] -bam bam.filelist
# NGSadmixture
# K is the number of groups, test out K=1 up to K=5 to start, -P = number of threads
angsd/misc/NGSadmix -likes genolike.beagle.gz -K 2 -P 20 -o all_admix2 -minMaf 0.05
# Visualize in R on command line or RStudio
R
# Get ID and pop info for each individual from a flat file with sample names and localities.
pop<-read.table("population_info.txt",as.is=T)
# Read inferred admixture proportions file (repeat this step and below for each K)
q<-read.table("all_admix2.qopt")
# save the plot
png(filename = "k2_all_barplot.png", width = 2400, height = 1400, res = 300)
# Plot them (ordered by population)
ord<-order(pop[,1])
# Set margins, increase the first number for more space
par(mar=c(10,4,1,1))
# more options for plot, las = 2 makes names vertical, 1 is horizontal
bp <- barplot(t(q)[,ord],
# col=2:10,
col=c("#882255","#88CCEE","#117733","#44AA99","#AA4499","#DDCC77","#CC6677","#332288"),
space=0,
border=NA,
xlab="",
ylab="Inferred Admixture Proportions (K=2)",
las=1)
# add unique group names vertically
text(x = tapply(1:nrow(pop), pop[ord,1], mean),
y = -0.05,
labels = unique(pop[ord,1]),
xpd=T,
srt = 90,
adj = 1,
cex = 0.8,
font=2)
# add lines to separate the populations
abline(v=cumsum(sapply(unique(pop[ord,1]), function(x){sum(pop[ord,1]==x)})),col=1,lwd=1.2)
dev.off()Make a PCA in ANGSD.
# install
conda install bioconda::pcangsd
pcangsd -h
# Genotype likelihoods using 32 threads, if pca does not converge then increase iterations
pcangsd -b all_bams/angsdput.beagle.gz -o des_mito_final_pca -t 32 --iter 2000
#Plot the results in R or RStudio
#open R
R
# import data and metadata
pop<-read.table("sample_info.txt")
C <- as.matrix(read.table("sample_pca.cov"))
e <- eigen(C)
# Calculate the percentage of variance
# e$values contains the eigenvalues from your eigen(C) call
var_explained <- e$values / sum(e$values) * 100
# Round to 2 decimal places for the labels
pc1_lab <- paste0("PC1 (", round(var_explained[1], 2), "%)")
pc2_lab <- paste0("PC2 (", round(var_explained[2], 2), "%)")
# Create a named vector: Assign specific colors to specific names
# Replace these examples with your actual group names from sample_info
group_colors <- c(
"country1" = "#882255",
"country2" = "#88CCEE",
"country3" = "#117733",
"country4" = "#7CEBD8",
"country5" = "#E182D1",
"country6" = "#DDCC77",
"country7" = "#CC6677",
"country8" = "#332288",
"country9" = "#999933"
)
# Map the colors to actual data rows
cols <- group_colors[as.character(pop[,1])]
# Adjust margins to put the legend outside of the plot
pdf("sample_pca.pdf", width = 10, height = 7)
par(mar = c(5, 4, 4, 12)) # To add extra space on the right (12)
# Plot
plot(e$vectors[,1], e$vectors[,2],
col=cols,
pch = 19,
cex = 1.8, # make dots larger
xlab= pc1_lab,
ylab=pc2_lab,
main="PCA of Allele Frequencies"
# run this section to zoom in on certain parts of the pca
# main = "PCA (zoomed in) ",
# xlim = c(-0.05, 0.05), # Adjust these numbers to zoom in/out
# ylim = c(-0.05, 0.05) # Adjust these numbers to zoom in/out
)
# Place legend in margins outside of the pca plot
legend(x = "topright",
inset = c(-0.3,0),
legend = names(group_colors),
col = group_colors,
pch = 19,
bty = "n",
cex = 0.7,
xpd = T)
dev.off()Pull out mitochondrial reads from .bams and run admixture/pca analysis with only mitochondrial data.
# filter for mito reads, sort, and index, NC_022423.1 is the name of the mitochodria in the reference genome
samtools view -b allbarcode01_mapped.bam NC_022423.1 | samtools sort -o barcode01_mito_sorted.bam
#index the reads
samtools index barcode01_mito_sorted.bam
# check that it worked, should show zeros for everything (or almost everything) and many for mitochondrial (NC022423.1 in this case)
samtools idxstats barcode01_mito_sorted.bam
# Go back and make the PCA and structure plots with these filesCalculate SAF and SFS. Here it is done for only chromosome 1.
#first generate .saf file, NC_071387.1 is chromosome 1
# This should be done by population
# here we use the reference in place of an ancestral state to polarize the data for a folded analysis
# -minInd = number of inidividuals needing data at a loci to be included, -uniqueOnly 1 tells ANGSD to discard any read that aligns to more than one location in the genome
angsd -bam pop1_bam.filelist -doSaf 1 -uniqueOnly 1 -minInd 10 -out pop1_folded -anc reference.fna -GL 1 -P 20 -r NC_071387.1
# 1D SFS: EM optimization with 4 threads
# add more iterations (-maxIter) to help model converge, use -tole 1e-6 for strict convergence
angsd/misc/realSFS pop1_folded.saf.idx -maxIter 4000 -P 8 -fold 1 -tole 1e-6 > pop1_folded.sfs
# look at the output
cat pop1_folded.sfs
#in R
sfs<-scan("pop1_folded.sfs")
barplot(sfs[-1])
# 2D SFS: need this for pairwise Fst calculation
angsd/misc/realSFS pop1_folded.saf.idx pop2_folded.saf.idx -fold 1 -P 8 -tole 1e-6 -maxIter 1000 > pop1_pop2.2dsfsUse the 2D SFS from the previous step to calculate Fst.
# Index the Fst values, combines SAF files and 2D SFS prior to calculate the per-site Fst probabilities
angsd/misc/realSFS fst index pop1_folded.saf.idx pop2_folded.saf.idx -sfs pop1_pop2.2dsfs -fstout pop1_pop2 -fold 1
# Global Fst
../angsd/misc/realSFS fst stats pop1_pop2.fst.idx
# example output
#nObs is number of observations or genomic sites passing filtering in both populations
#unweighted FST: ratio for every site individually, and then takes the average of those ratios
#weighted FST: the global FST by summing up all the numerators across the genome and dividing them by the sum of all the denominators. Use and report this number.
-> Assuming idxname:pop1_pop2.fst.idx
-> Assuming .fst.gz file: pop1_pop2.fst.gz
-> FST.Unweight[nObs:8136]:0.030263 Fst.Weight:0.316762
0.030263 0.316762
# Windowed Fst
# Example: 50kb windows with a 10kb step
realSFS fst stats2 pop1_pop2.fst.idx -win 50000 -step 10000 > pop1_pop2.txt# make the thetas.idx file
angsd/misc/realSFS saf2theta pop1_folded.saf.idx -sfs pop1_folded.sfs -outname pop1
# estimate the stats
# Global estimate, output "pop1.thetas.idx.pestPG"
thetaStat do_stat pop1.thetas.idx
## view the output
cat pop1.thetas.idx.pestPG
# Sliding window estimate
angsd/misc/thetaStat do_stat pop1.thetas.idx -win 50000 -step 10000 -outnames pop1_thetas_windows