The goal of alignstatplot is to provide a simple tool for performing sequence alignment, descriptive analysis, and innovative visualization of the generated results. alignstatplot provides some new analyses that may aid researchers in identifying genic regions associated with gene variability and can be combined with other analysis tools for studying diversity. In this documentation, we will demonstrate how to use it and how it can be combined with other well-known R packages.
The tool is available online :
https://bioinformatics.um6p.ma/AlignStatPlot/
Prefer Python, or want a point-and-click app instead of R code? alignstatplot-py is a companion Python package that reimplements the same categories of analysis (sequence stats, nucleotide diversity, distance/phylogenetic trees, PCA-based SNP clustering, and the circular/dashboard plots) on NumPy/pandas/ SciPy/scikit-learn/Biopython/Matplotlib. It ships as:
- a
pip-installable library + CLI (alignstatplot run --fasta ... --anno ...), - a desktop GUI (
alignstatplot gui, PySide6/Qt) for uploading a FASTA file and browsing plots, and - a local web GUI (
alignstatplot web) with a drag-and-drop browser frontend.
This R package remains the reference implementation and is still actively maintained; see alignstatplot-py's README for how closely its results match this package's numerically (same statistics and plot types, not bit-for-bit identical, since it doesn't depend on R/ClustalW/FactoMineR/ape).
This is an updated version focused on making the package lighter and easier to install, and giving users more control over plot appearance:
- Fewer, lighter dependencies. Dropped the Bioconductor dependency (
ggtree,BiocManager) along withtidyverse,plyr,reshape2,pals,factoextra, andphytools— all reimplemented internally on top of packages already required. The package now installs from CRAN alone, no Bioconductor configuration needed. - More control over plots. Plotting functions (
distanceHeatmap,nucTableHeatmap,plotTreeWithRuler,SNPClusterPlotPCAMap, and others) now acceptcolors/palette, font-size, and theme-style arguments (with the previous look kept as the default), and return their plot object instead of only drawing/saving as a side effect. - Input validation. Core analysis functions (
getSeqInfo,alignment2Table,AlignmentStatsPerSeq,nucFrequency,seqTableToBinary,getDistanceMatrixTabel,getTree,SNPCluster, and others) now reject malformed input with an informative error at the point of the call, instead of failing deep inside a downstream plotting step. - Regression test suite. Added a
testthat-based suite covering the analysis pipeline.
See REFACTOR_PLAN.md in the repository for the full, step-by-step account of this work.
You can install the development version of alignstatplot like so:
library(devtools)
install_github("https://github.com/AlsammanAlsamman/alignstatplot")rm(list = ls())
library(alignstatplot)
library(ggplot2) #For saving FiguresAll the analyses can be performed using one function
SeqFile<-system.file("extdata","Example_Small.fasta",package = "alignstatplot")
AnnoFile<-system.file("extdata","Example_Small.anno",package = "alignstatplot")
MaxCluster = 4 #Numbe of SNP clusters
MinimumClusterLength = 3 #Length of the SNP cluster to be considered
AlignMethod<- "ClustalW" #"ClustalW", "ClustalOmega", and "Muscle"
OutFolder<-"output" #Output folder
Verbose<-T
MaxMissPer=0.2 #Maximum missing allowed in sequence alignment
consZoomFactor <- 3 #Zooming factor for circle alignment visualization for the consensus compared to the sequences
cex.SeqLabels<-0.5 #Sequence labels font sizes
# alignstatplot(SeqFile,
# AlignMethod,
# AnnoFile,
# OutFolder,
# MaxMissPer,
# MaxCluster,
# MinimumClusterLength,
# Verbose)Alignstatplot can be used to plot sequence alignment file. The function can import sequence alignment in FASTA or clustalw
AlignFile<-system.file("extdata","Example_Sequences_Aligned.aln",package = "alignstatplot")
#AlignFile<-system.file("extdata","Example_Sequences_Aligned.fasta",package = "alignstatplot")
SeqFormat<-"clustalw" #clustalw or fasta #input format
SeqFontSize<-0.5 # font label size
#pdf(paste0(OutFolder,"/","AlignPlot.pdf"),width = 15,height = 15)
plotAlignCircle(AlignFile,SeqFormat = "clustalw",SeqFontSize = 0.5)
#dev.off()Example files can be used from the package.
SeqFile<-system.file("extdata","Example_Small.fasta",package = "alignstatplot")
#Simple Fasta File
#>Gene
#>ACTCGCTCGACTACGACCAT
AnnoFile<-system.file("extdata","Example_Small.anno",package = "alignstatplot")
#Annotation file, could contain any kind of information regarding gene structure
#It contains GeneName Region Start End direction forward or reverse
#gene1 intron 1 14 forward
#gene1 exon 232 474 forward
#Sequence alignment program to be used
AlignMethod<- "ClustalW" #"ClustalW", "ClustalOmega", and "Muscle"
#Output Folder
OutFolder<-"output"
# Prepare environment
#dir.create(file.path("output"), showWarnings = FALSE)Performing sequence alignment using the chosen tool.
#Align DNA sequences using clustalw
myClustalWAlignment <- seqAlign(SeqFile,AlignMethod)
#> use default substitution matrix
#extract sequence information
SeqInfo<-getSeqInfo(SeqFile)
#Convert sequence to list
SeqAligned<-alignment2Fasta(myClustalWAlignment,SeqInfo)
#Save aligned sequence in fasta format
alignment2Fasta(myClustalWAlignment,SeqInfo,paste(OutFolder,"Alignment.fasta",sep = "/"))Alignment statistics contains different aspects regarding sequence alignment results.
StatsTable<-AlignmentStatsPerSeq(SeqInfo,SeqAligned)
#colnames(StatsTable)
# [1] "SeqNo" "Sequence.Name" "Sequence.Length" "A" "T" "C"
# [7] "G" "Gap" "Gap.Percentage" "GC" "GC.Percentage"
head(StatsTable,5)
#> SeqNo Sequence.Name Sequence.Length A T C G Gap Gap.Percentage GC
#> 1 1 gene1 969 268 194 202 305 72 6.92% 507
#> 2 2 gene2 973 257 197 215 304 68 6.53% 519
#> 3 3 gene3 910 244 176 196 294 131 12.58% 490
#> 4 4 gene4 955 240 179 225 311 86 8.26% 536
#> 5 5 gene5 931 233 176 216 306 110 10.57% 522
#> GC.Percentage
#> 1 52.32%
#> 2 53.34%
#> 3 53.85%
#> 4 56.13%
#> 5 56.07%Plot Sequence alignment with consensus and links , Warining more than 20 genes it will not recognizable
#You can use drawConsWithNoGenes instead iof you have more than 20
#pdf(paste0(OutFolder,"/","SeqAlignmentCircleWithLinks.pdf"))
drawConsWithGenes(SeqInfo,SeqAligned)
#dev.off()geneLabels/consensusLabel let you override the displayed names (the
underlying SeqInfo$Name is still used for data linking), and each text
element on the plot has its own font-size control: cex.SeqLabels for the
gene name labels, cex.ConsLabel for the consensus label (defaults to
cex.SeqLabels), and cex.RulerLabels/cex.ConsRulerLabels for the
sequence and consensus base-pair ruler tick labels respectively.
drawConsWithGenes(SeqInfo,SeqAligned,
geneLabels = paste0("Gene ", seq_len(nrow(SeqInfo))),
consensusLabel = "Consensus (all genes)",
cex.SeqLabels = 0.6, cex.ConsLabel = 0.9,
cex.RulerLabels = 0.4, cex.ConsRulerLabels = 0.6)
Plot Sequence alignment with consensus and No links
#pdf(paste0(OutFolder,"/","SeqAlignmentCircleWithNoLinks.pdf"),width = 15,height = 15)
drawConsWithNoGenes(SeqInfo,SeqAligned,cex.SeqLabels = 1)
#dev.off()geneLabels overrides the displayed sequence names; cex.SeqLabels and
cex.bpLabels independently control the sequence-name and base-pair-ruler
text sizes.
drawConsWithNoGenes(SeqInfo,SeqAligned,
geneLabels = paste0("Gene ", seq_len(nrow(SeqInfo))),
cex.SeqLabels = 1, cex.bpLabels = 0.5)
It can handle tenth of genes and still showing the alignment differences.
knitr::include_graphics(paste0("Example_Out","/SeqAlignmentCircleWithNoLinks.png"))
#Calculate Distance Table (Similarity distance)
DistTable<-getDistanceMatrixTabel(SeqInfo,myClustalWAlignment)
#Save Distance Table
#write.table(DistTable,paste0(OutFolder,"/Distance_Table.tsv"),sep="\t")
#Calculate Phylogenetic Tree
myTree<-getTree(DistTable)
#Save Phylogenetic Tree as Text to Folder
#writeNexus(myTree,paste0(OutFolder,"/Tree.nexus"))
#Save Tree Summary"))
#sink(file=paste0(OutFolder,"/TreeSummary.txt"))
summary(myTree)
#>
#> Phylogenetic tree: myTree
#>
#> Number of tips: 10
#> Number of nodes: 8
#> Branch lengths:
#> mean: 0.0984519
#> variance: 0.00642544
#> distribution summary:
#> Min. 1st Qu. Median 3rd Qu. Max.
#> 0.006071412 0.034203569 0.054481207 0.174942542 0.233735569
#> No root edge.
#> Tip labels: gene1
#> gene2
#> gene3
#> gene4
#> gene5
#> gene6
#> gene7
#> gene8
#> gene9
#> gene10
#> No node labels.
#sink()# The similarity matrix as a heatmap. fontsizescale controls the row/column
# label size (0, the default, picks a size automatically based on the number
# of sequences).
distanceHeatmap(DistTable)
#The phylogenetic tree
plotTreeWithRuler(SeqInfo,myClustalWAlignment)
#The phylogenetic tree in combine with gene structure
plotTreeWithGenes(SeqInfo,myClustalWAlignment,AnnoFile)
#The phylogenetic tree in combine with simialrity matrix
#fontsizescale controls the tip/label size the same way as above.
plotSimilarityMatrixWithTree(SeqInfo,myClustalWAlignment)
#> NULL
The sequence alignment data are used for gene clustering using principle component analysis
plotPCA(SeqInfo,myClustalWAlignment,ncluster = 4)
#> Warning: argument frame is deprecated; please use ellipse instead.
#> Warning: argument frame.type is deprecated; please use ellipse.type instead.
# ggsave(plot=p,paste0(OutFolder,"/","Genes_PCA.pdf"),device = "pdf",
# width = 10, height = 8, dpi = 150, units = "in")#Convert alignment To Table
SeqAlignedTable<-alignment2Table(SeqInfo,SeqAligned)
#Filter Table for missing and Mono
SeqAlignedTableFiltered<-nucTableFilter(SeqAlignedTable,MaxMissPer = 0.2,removeMono = T)
#Calculate Frequency
NucCount<-nucFrequency(SeqAlignedTableFiltered)
NucCount[,1:5]
#> N9 N13 N17 N29 N31
#> A 0.1 0.0 0.0 0.1 0.3
#> C 0.0 0.7 0.4 0.0 0.0
#> T 0.0 0.1 0.4 0.0 0.0
#> G 0.7 0.0 0.0 0.9 0.7
#> N 0.2 0.2 0.2 0.0 0.0Heatmap and bar plots of the nucleotides frequency across the sequence alignment
#plot sample
#nucTableHeatmap(SeqAlignedTableFiltered[,1:50],cex.SeqLabels = 5,cex.NucLabels = 5)
#Chart of the nucleotides frequency across the sequence alignment
#nucFrequencyPlot(NucCount[,1:50],F)
#combine heatmap and charts
#plot first 50 nucleotide variations
nucTableFreqHeatmap(SeqAlignedTableFiltered[,1:50])
#If the sequence is very long it can be split and saved to a folder
# dir.create(file.path(paste(OutFolder,"HeatMap_Freq",sep="/")), showWarnings = FALSE)
# # If The sequence was very long this can be used
# PlotList<-nucTableFreqHeatmapSplit(SeqAlignedTableFiltered,50)
# outPdfDir<-paste(OutFolder,"/","HeatMap_Freq","/",sep="")
# # #Save Plots
# saveSeqPlotList(PlotList,outPdfDir)#The figure is big we can include in the documentation using this
knitr::include_graphics(paste0("Example_Out","/NucLogo.png"))
Binary format is very important step in several data analysis, where most of the algorithms are not accepting allele-based format such population structure analysis. There are several ways to convert nucleotide formatted data (sequences) to binary.
Choosing reference is helping in searching for only variations that differentiate between two cases (control/patient or resistant/sensitive).
#One reference - based conversion
RefsNames<-c("gene1")
#Using One Reference
SeqBinaryTableOneRef<-seqTableToBinary(SeqAlignedTableFiltered,RefsNames,RemoveNonRefNuc = T)
SeqBinaryTableOneRef[1:5,1:5]
#> N9 N13 N17 N29 N31
#> gene1 1 1 1 1 1
#> gene2 0 1 0 1 0
#> gene3 1 1 1 1 0
#> gene4 1 1 0 0 0
#> gene5 1 1 0 1 0In such cases, more than one reference can be used, with shared nucleotide variations across reference genotypes used to create a shared reference genotype that can be used to convert all data. This is a stricter conversion in which only a few shared variations are used; this can be used to search for specific interesting stable variations across the control cases.
RefsNames<-c("gene1","gene2","gene3")
SeqBinaryTableMoreRef<-seqTableToBinary(SeqAlignedTableFiltered,RefsNames,RemoveNonRefNuc = T)
SeqBinaryTableMoreRef[1:5,1:5]
#> N13 N29 N36 N60 N84
#> gene1 1 1 1 1 1
#> gene2 1 1 1 1 1
#> gene3 1 1 1 1 1
#> gene4 1 0 1 1 1
#> gene5 1 1 1 1 1The fisrt step will be to remove non-biallailic variations (tri-, tetra-, ..) and assign 0 or 1 for alleles. This is the default way for data conversion.
#report bialleles locations.
biallelicNuc<-getBiallelicByFreq(NucCount)
#Create a reference genotype.
GenotypeRef<-getRefGenotypeForbiallelic(NucCount,biallelicNuc)
#Convert data using this genotype as a reference (restricted to the
#biallelic columns GenotypeRef actually has a reference allele for)
SeqBinaryTableByFreqRef<-seqTableToBinary(SeqAlignedTableFiltered[,biallelicNuc],GenotypeRef,RemoveNonRefNuc = T,RefsNames = F)
SeqBinaryTableByFreqRef[1:5,1:5]
#> N9 N13 N17 N29 N31
#> gene1 1 1 1 1 0
#> gene2 0 1 0 1 1
#> gene3 1 1 1 1 1
#> gene4 1 1 0 0 1
#> gene5 1 1 0 1 1#SNP Clustering Nucleotide variation clustering is a PCA-based analysis that clusters SNP data across genes. It aids in identifying the variation that occurs across genes, where some alleles tend to appear together.
#Allele Based Variations
Cluster<-SNPCluster(SeqAlignedTableFiltered)Phylogenetic tree constructed showing the clustering association across the loci (SNPs)
SNPClusterPlot3DTree(Cluster,60)
##
One Dimension Tree 2D phylogenetic tree
SNPClusterPlot1DTree(Cluster)
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
##
PCA Plot
SNPClusterPlotPCAMap(Cluster)
##
Drawing genes with cluster Map Draw the SNP location and assign the
cluster to the genes. The function takes a sequence, a sequence
alignment object, and the number of base pairs to consider as a cluster.
For example, if three nucleotides “ATC” belong to cluster 1 and the
cluster’s minimum bp is 4, it will not be drawn on gene. As a result,
the cluster should contain more than four base pairs close to each other
to be drawn, and it could be regarded as a motif segment across the
genes. SNP clustering could be used to look for consistent patterns that
are unaffected by differences in gene length.
SNPClusterPlot(SeqInfo,SeqAligned,Cluster,MaxCluster = 3,MinimumClusterLength = 3)
Add the phylogenetic tree to the plot.
# With phylogenetic tree
SNPClusterPlotWithTree(SeqInfo,myClustalWAlignment,Cluster,MaxCluster = 3, MinimumClusterLength = 3)
With hundreds of SNPs, the dendrogram above gets crowded well before it gets
uninformative – most SNPs barely move the clustering at all. SNPClusterImpact()
scores every SNP by how much it actually drives the retained PCA/MCA dimensions the
clustering above is built on (its contribution, weighted by how much variance each
dimension explains), so low-impact SNPs can be identified and filtered out before
re-clustering.
ImpactTable<-SNPClusterImpact(Cluster)
head(ImpactTable)
#> SNP Impact Cos2 Cluster
#> 930 N930 0.3123553 0.05119092 3
#> 931 N931 0.3093967 0.04986123 3
#> 921 N921 0.3049479 0.05054261 2
#> 926 N926 0.2995618 0.05110840 3
#> 923 N923 0.2991817 0.05072784 2
#> 913 N913 0.2892038 0.04894040 3plotSNPClusterImpact() sorts SNPs by impact and shows both the score itself and the
cumulative share of total impact captured, so a cutoff can be chosen by eye instead of
guessed – the dashed line marks filterHighImpactSNPs()’s default cutoff (the elbow of
the curve), and can be overridden with an explicit topN or minImpact.
plotSNPClusterImpact(ImpactTable)
Re-clustering on just the kept SNPs gives a much less crowded, equally structured dendrogram:
SeqAlignedTableFilteredHQ<-filterHighImpactSNPs(SeqAlignedTableFiltered, ImpactTable)
dim(SeqAlignedTableFilteredHQ)
#> [1] 10 40ClusterHQ<-SNPCluster(SeqAlignedTableFilteredHQ)SNPClusterPlot1DTree(ClusterHQ)
Beyond the alignment and clustering tools above, alignstatplot provides classic
population-genetics summary statistics and a set of publication-style plots that all
share one consistent theme (theme_alignstatplot()).
Segregating sites, nucleotide diversity (π), Watterson's θ, and the transition/transversion ratio, computed directly from the aligned sequence table.
nucDiversityStats(SeqAlignedTableFiltered)
#> NumSequences AlignmentLength SegregatingSites NucleotideDiversity
#> 1 10 931 931 0.7305645
#> WattersonsTheta TsTvRatio
#> 1 0.3534858 0.8181818positionConservation() reports the majority allele, conservation, entropy, and
missingness at every alignment position; plotConservationTrack() and
plotVariantDensity() plot it as a track along the alignment.
PositionStats<-positionConservation(SeqAlignedTableFiltered)
plotConservationTrack(PositionStats)
plotVariantDensity(SeqAlignedTableFiltered, windowSize = 50)
A distributional complement to the distanceHeatmap() above.
plotIdentityDistribution(DistTable)
plotSeqStatsSummary(StatsTable)
plotBaseComposition(NucCount)
regionStats() aggregates positionConservation()'s output within each annotated
region (e.g. exon/intron), so conservation and variant rate can be compared across
region types directly.
RegionStats<-regionStats(PositionStats, AnnoFile)
plotRegionStats(RegionStats)
A single patchwork-combined figure suitable as a supplementary summary figure.
plotSummaryDashboard(SeqInfo, SeqAligned, DistTable)
library(poppr)
#> Loading required package: adegenet
#> Loading required package: ade4
#> Registered S3 method overwritten by 'vegan':
#> method from
#> rev.hclust dendextend
#>
#> /// adegenet 2.1.7 is loaded ////////////
#>
#> > overview: '?adegenet'
#> > tutorials/doc/questions: 'adegenetWeb()'
#> > bug reports/feature requests: adegenetIssues()
#> Registered S3 method overwritten by 'pegas':
#> method from
#> print.amova ade4
#> This is poppr version 2.9.3. To get started, type package?poppr
#> OMP parallel support: available
SeqBinaryTableByFreqRef[1:10,1:10]
#> N9 N13 N17 N29 N31 N35 N36 N37 N38 N39
#> gene1 1 1 1 1 0 - 0 1 0 1
#> gene2 0 1 0 1 1 1 0 0 0 0
#> gene3 1 1 1 1 1 0 0 0 1 1
#> gene4 1 1 0 0 1 0 0 0 0 1
#> gene5 1 1 0 1 1 0 0 0 0 1
#> gene6 1 1 0 1 1 0 0 0 0 1
#> gene7 - - - 1 1 - 0 1 1 -
#> gene8 1 0 1 1 0 1 0 0 1 0
#> gene9 1 1 1 1 1 0 0 1 0 1
#> gene10 - - - 1 0 0 0 0 0 1
length(SeqInfo$Name)
#> [1] 10
rownames(SeqBinaryTableByFreqRef)
#> [1] "gene1" "gene2" "gene3" "gene4" "gene5" "gene6" "gene7" "gene8"
#> [9] "gene9" "gene10"
monpop <- df2genind(SeqBinaryTableByFreqRef,
ploidy = 1,
ind.names = SeqInfo$Name,
sep="")
summary(monpop)
#>
#> // Number of individuals: 10
#> // Group sizes: 10
#> // Number of alleles per locus: 3 3 3 2 2 3 1 2 2 3 2 3 3 2 3 3 3 3 2 1 2 2 2 2 1 2 3 2 2 2 3 1 1 1 1 1 1 2 2 1 2 2 1 2 2 2 2 2 2 2 2 2 1 2 2 1 2 2 2 2 1 2 2 1 1 2 1 2 1 1 2 2 2 2 2 1 2 2 2 3 2 2 2 2 2 3 2 1 2 2 2 2 2 1 3 2 1 2 2 2 2 3 2 3 3 2 2 2 1 2 1 1 2 2 2 2 1 2 2 2 2 2 2 1 2 1 1 2 1 2 1 2 1 2 2 1 2 1 2 1 1 1 2 2 1 2 2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 1 2 1 2 1 2 2 1 2 2 2 1 1 2 2 2 2 1 2 2 3 3 2 2 3 3 2 2 3 2 3 2 3 2 3 3 3 2 3 2 2 2 1 1 1 2 1 2 2 1 2 2 1 1 2 1 1 2 2 2 1 2 2 1 1
#> // Number of alleles per group: 428
#> // Percentage of missing data: 0 %
#> // Observed heterozygosity: 0
#knitr::knit_exit()If you use this package, please cite:
Alsamman, A.M., El Allali, A., Mokhtar, M.M., Al-Sham'aa, K., Nassar, A.E., Mousa, K.H. and Kehel, Z., 2023. AlignStatPlot: An R package and online tool for robust sequence alignment statistics and innovative visualization of big data. PLoS ONE, 18(9), p.e0291204. https://doi.org/10.1371/journal.pone.0291204
Or in R:
citation("alignstatplot")