Basic sequence analysis and manipulation

Practical 2: Genomics, systems biology and multi-omics

Learning objectives

Parse FASTA records; validate sequence symbols; calculate called-base GC content; construct a reverse complement.

Biological question and data

How do file structure and unknown bases affect a sequence summary?

FASTA inputs contain only A, C, G, T and N; these teaching records have no quality scores. N is unknown and is excluded from the GC denominator.

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. Parse explicit records

Use the local teaching parser, which rejects missing headers, duplicate IDs and invalid bases.

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"))
sequences <- read_teaching_fasta(c(">demo1","ACGTNN","GCGC",">demo2","TTTACG"))
sequences
       demo1        demo2 
"ACGTNNGCGC"     "TTTACG"
Show R code 
stopifnot(nchar(sequences["demo1"])==10)

Wrapped sequence lines are joined within their records. Sequence length includes unknown bases.

2. Calculate called-base summaries

Count unknowns separately. An all-unknown sequence has no estimable called-base GC proportion.

Show R code 
summary_one <- function(s) {
 x <- strsplit(s,"")[[1]]; called <- x%in%c("A","C","G","T")
 c(length=length(x),called=sum(called),unknown=sum(!called),
 GC=if(any(called))mean(x[called]%in%c("G","C")) else NA_real_)
}
result <- t(vapply(sequences,summary_one,numeric(4)))
result
      length called unknown        GC
demo1     10      8       2 0.7500000
demo2      6      6       0 0.3333333
Show R code 
stopifnot(abs(result["demo1","GC"]-6/8)<1e-12,
 is.na(summary_one("NNN")["GC"]))

demo1 has ten characters but eight called bases; six called bases are G or C.

3. Reverse-complement the sequence

Reverse the base order as well as complementing each symbol. Keep unknown N as N.

Show R code 
reverse_complement <- function(s) {
 x <- rev(strsplit(s,"")[[1]])
 paste(unname(c(A="T",C="G",G="C",T="A",N="N")[x]),collapse="")
}
rc <- vapply(sequences,reverse_complement,character(1))
rc
       demo1        demo2 
"GCGCNNACGT"     "CGTAAA"
Show R code 
stopifnot(identical(unname(vapply(rc,reverse_complement,character(1))),unname(sequences)))

Reverse complementing twice restores the original sequence, a useful orientation check.

4. Make the denominator visible

Unknown bases change what can be calculated. Compare GC fractions using all supplied characters and only called bases in a hypothetical sequence.

Show R code 
seq_chars <- strsplit('GCGCNNNNAT','')[[1]]
called <- seq_chars %in% c('A','C','G','T')
c(all_character_GC=mean(seq_chars %in% c('G','C')),
 called_base_GC=mean(seq_chars[called] %in% c('G','C')),
 called_fraction=mean(called))
all_character_GC   called_base_GC  called_fraction 
       0.4000000        0.6666667        0.6000000

These summaries answer different questions. Excluding unknown positions gives a called-base composition; it does not recover the composition of unobserved positions. Report the unknown fraction alongside it.

5. Preserve orientation and record context

Reverse complementation changes orientation, while reversing alone does not complement bases. A double reverse complement should recover a valid input sequence under the supported alphabet. Before comparing records, check whether each represents a complete sequence or a region, what the header identifies and whether ambiguous symbols are supported by the parser. Restricting to A/C/G/T/N is a teaching contract, not the complete vocabulary of all biological sequence files.

Interpretation challenge

Two sequences have identical called-base GC fractions but one has 40% unknown bases. Are their measurement qualities equivalent?

No. The fraction describes only called positions. Unknown coverage, sequence length and the observation process must accompany the composition statistic.

Interpretation and limits

Called-base GC is a descriptive proportion, not sequencing quality or coding potential. Unknown bases and selected regions can make whole-record comparisons misleading. This parser is intentionally small, not a general FASTA library.

Exercises

  1. Recalculate GC including N in the denominator; explain the change.
  2. Add a duplicate header.
  3. Why does complementing without reversing give a different object?
  1. It decreases to 6/10, treating unknowns as if they were called non-GC bases.
  2. The parser rejects the ambiguous identifier.
  3. A reverse complement changes strand orientation; complement alone preserves the original character order.

Continue learning

B04 Programming and B05 Genomics.

Continue with sequence similarity. Return to the practical index or the unified map.