Alignment and conserved-region analysis

Practical 5: Genomics, systems biology and multi-omics

Learning objectives

Inspect a global-alignment score; distinguish gaps from nucleotide differences; calculate conservation across an explicitly supplied multiple alignment.

Biological question and data

How do alignment and missingness decisions change the apparent conservation of a region?

The local helper uses Needleman–Wunsch global alignment with fixed match, mismatch and linear-gap scores. The multiple alignment is supplied, not inferred by this pairwise helper.

No biological download is required. Numerical inputs are hypothetical, not observed animal data. The runnable examples use base R. Count and abundance units are specified where used.

1. Align two sequences under stated scores

Global alignment compares the entire supplied sequences. Different scoring choices can change an optimal alignment.

Show R code 
helper_dir <- if (dir.exists("R")) "R" else "../R"
stopifnot(file.exists(file.path(helper_dir,"sequence_helpers.R")))
source(file.path(helper_dir,"sequence_helpers.R"))
a <- "ACGTACGT"; b <- "ACGACGT"
fit <- global_alignment(a,b)
fit
$a
[1] "ACGTACGT"

$b
[1] "ACG-ACGT"

$score
[1] 12
Show R code 
stopifnot(gsub("-","",fit$a)==a,gsub("-","",fit$b)==b,
 nchar(fit$a)==nchar(fit$b),fit$score==12)
self <- global_alignment(a,a)
stopifnot(self$score==2*nchar(a),self$a==a,self$b==a)

Seven matches and one gap produce a score of 12 under +2/−1/−2 scoring. The score is not an evolutionary distance.

2. Inspect an existing multiple alignment

A column is conservatively called conserved only when every sequence supplies a nucleotide and all those nucleotides match.

Show R code 
aligned <- c(s1="ACGTACGT",s2="ACG-ACGT",s3="ACGTTCGN")
M <- do.call(rbind,strsplit(aligned,""))
called <- M%in%c("A","C","G","T"); dim(called) <- dim(M)
fully_called <- colSums(called)==nrow(M)
conserved <- apply(M,2,function(z) all(z%in%c("A","C","G","T")) && length(unique(z))==1)
data.frame(position=seq_len(ncol(M)),fully_called,conserved)
  position fully_called conserved
1        1         TRUE      TRUE
2        2         TRUE      TRUE
3        3         TRUE      TRUE
4        4        FALSE     FALSE
5        5         TRUE     FALSE
6        6         TRUE      TRUE
7        7         TRUE      TRUE
8        8        FALSE     FALSE
Show R code 
stopifnot(sum(fully_called)==6,sum(conserved)==5)

A gap and an unknown base are excluded from fully called columns; a genuine nucleotide difference breaks conservation in another column.

3. Compare denominators

Calculate both the fraction of all alignment columns and the fraction of fully called columns that are conserved.

Show R code 
c(conserved_of_all=sum(conserved)/ncol(M),
 conserved_of_fully_called=sum(conserved)/sum(fully_called))
         conserved_of_all conserved_of_fully_called 
                0.6250000                 0.8333333

These fractions answer different questions. Report the denominator and missingness rather than choosing the larger value.

4. Inspect scoring sensitivity

Different gap costs can favor a mismatch or an insertion/deletion explanation. Compare the score contributions of two hypothetical candidate alignments under two rules. These are supplied candidate counts, not new optimal alignments.

Show R code 
candidates <- data.frame(alignment=c('no_gap','two_gaps'),
 matches=c(6,7),mismatches=c(2,0),gap_columns=c(0,2))
score <- function(gap_cost) with(candidates,2*matches-mismatches-gap_cost*gap_columns)
data.frame(candidates,score_gap1=score(1),score_gap4=score(4))
  alignment matches mismatches gap_columns score_gap1 score_gap4
1    no_gap       6          2           0         10         10
2  two_gaps       7          0           2         12          6

The preferred supplied candidate can change with the cost. The existing alignment helper uses a simple gap rule; affine-gap methods distinguish opening and extending gaps and address a different scoring model.

5. Connect an alignment to downstream claims

Conservation depends on which positions and sequences are comparable. A gap can reflect evolutionary change, incomplete sequence or alignment uncertainty. Filtering uncertain columns changes the denominator for identity and distance. For downstream phylogeny or structure mapping, preserve the relationship between aligned positions and original residue identifiers.

Assess whether sequences are comparable homologous regions before calculating evolutionary distances. A mathematically optimal score is optimal under the supplied scoring rule, not a certificate of true history.

Interpretation challenge

Why should the highest-scoring alignment not automatically be described as the true evolutionary alignment?

The score encodes limited assumptions. Repeats, rearrangements, sequence errors and evolutionary processes can admit multiple plausible alignments. Biological context and sensitivity are needed.

Interpretation and limits

The example does not identify an actual functional element. Conservation depends on alignment quality, phylogenetic diversity and evolutionary processes. Production multiple aligners and evolutionary conservation models need their own method choices and diagnostics.

Exercises

  1. Change the gap score and inspect the result.
  2. Why is N not evidence for a conserved base?
  3. Why does conservation alone not prove function?
  1. Alternative alignments can become optimal, particularly where mismatch and gap alternatives compete.
  2. Its nucleotide is unknown.
  3. Context, species sampling, evolutionary constraints and biological evidence determine the interpretation.

Continue learning

B05 Genomics.

Continue with phylogenetic reasoning. Return to the practical index or the unified map.