BIO331 – Lab 07: Bioinformatics and Working with Strings of Biological Data

Bioinformatics

Author

Dee Ruttenberg

Published

March 18, 2026

Today we’re going to apply our work on R to Bioinformatic data sets. The typical output of a DNA sequencing experiment is a fastq file, consisting of often millions of configs which need to be stitched together (do you remember the two ways we do that?). Fastq files are different from fasta files as they contain quality scores – this is useful to see how reliable your fastq data is. The first tool we are going to use is Fastqc, to interpret the results of our fastq files.

Quality Control

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
Bioconductor version '3.19' is out-of-date; the current release version '3.23'
  is available with R version '4.6'; see https://bioconductor.org/install
BiocManager::install("Rqc")
Bioconductor version 3.19 (BiocManager 1.30.27), R 4.4.3 (2025-02-28)
Warning: package(s) not installed when version(s) same as or greater than current; use
  `force = TRUE` to re-install: 'Rqc'
Installation paths not writeable, unable to update packages
  path: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
  packages:
    aplot, bitops, bookdown, caTools, class, cluster, colorspace, combinat,
    cpp11, devtools, dichromat, diffobj, doBy, ellipsis, emmeans, estimability,
    expm, factoextra, FactoMineR, Formula, fracdiff, furrr, future, geiger,
    gert, ggDNAvis, ggfun, ggplot2, ggpubr, ggsci, ggthemes, gh, gridExtra,
    Hmisc, htmlTable, httr2, igraph, jose, KernSmooth, lattice, lazyeval,
    listenv, litedown, lme4, MASS, Matrix, multcompView, mvtnorm, nlme, nnet,
    optimParallel, pak, parallelly, pkgdown, pkgload, Rcpp, RcppArmadillo,
    RcppThread, RCurl, renv, restfulr, roxygen2, rpart, rsconnect, RSQLite,
    rstatix, s2, segmented, seqinr, seqmagick, sessioninfo, sf, shiny,
    snowflakeauth, sourcetools, sp, spatial, survival, swagger, tidytree,
    urlchecker, usethis, vegan, webutils, XML, yulab.utils, zip, zoo
Old packages: 'backports', 'bit64', 'broom', 'bslib', 'callr', 'cli', 'clipr',
  'curl', 'data.table', 'dbplyr', 'dplyr', 'fs', 'glue', 'httr', 'knitr',
  'magrittr', 'openssl', 'processx', 'ps', 'purrr', 'ragg', 'readxl', 'rlang',
  'rmarkdown', 'rstudioapi', 'S7', 'selectr', 'stringi', 'tinytex', 'vctrs',
  'vroom', 'withr', 'xfun', 'xml2'

We will use a specific small dataset used specifically for teaching bioinformatics. Rqc will output an .html file (essentially a web page). For this class, I will ask that you understand specifically the results of ‘Average Quality’, ‘Cycle-Specific Average Quality’, ‘Cycle-specific Quality Distribution’, and ‘Cycle-Specific Base Call Proportion’.

library(Rqc)
Loading required package: BiocParallel
Loading required package: ShortRead
Loading required package: BiocGenerics

Attaching package: 'BiocGenerics'
The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs
The following objects are masked from 'package:base':

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
    Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
    tapply, union, unique, unsplit, which.max, which.min
Loading required package: Biostrings
Loading required package: S4Vectors
Loading required package: stats4

Attaching package: 'S4Vectors'
The following object is masked from 'package:utils':

    findMatches
The following objects are masked from 'package:base':

    expand.grid, I, unname
Loading required package: IRanges
Loading required package: XVector
Loading required package: GenomeInfoDb

Attaching package: 'Biostrings'
The following object is masked from 'package:base':

    strsplit
Loading required package: Rsamtools
Loading required package: GenomicRanges
Loading required package: GenomicAlignments
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: 'MatrixGenerics'
The following objects are masked from 'package:matrixStats':

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Attaching package: 'Biobase'
The following object is masked from 'package:MatrixGenerics':

    rowMedians
The following objects are masked from 'package:matrixStats':

    anyMissing, rowMedians
Loading required package: ggplot2
folder <- system.file(package="ShortRead", "extdata/E-MTAB-1147")
rqc(path = folder, pattern = ".fastq.gz")
'/private/var/folders/v6/3txfs1g125n_x8nnndk4t60r0000gn/T/RtmpGMhlwz/rqc_report.html' has been created.

Bioconductor

Bioconductor and Biostrings are especially powerful tools to examine fasta files in R. In this section, we’ll learn some tools to manipulate them.

library(Biostrings)
seqs <- readDNAStringSet("Full_cmv_259sequences_Aln.fasta")

As Fasta files are labeled, you can identify specific sequences, either by loading a specfic label or identifying a keyword:

my_seq <- seqs["Europe_United Kingdom_GQ221974.1_Human_herpesvirus_5_strain_3157_complete_genome"]
names_with_uk <- grep("United Kingdom", names(seqs), ignore.case = TRUE)
UK_seqs <- seqs[names_with_uk]
UK_subset <- subseq(UK_seqs, start = 20000, end = 25000)

You can, using R, code an aligned genome as a matrix, in order to identify snps.

mat <- as.matrix(seqs)
#mat <- as.matrix(UK_seqs)
snp_positions <- which(apply(mat, 2, function(col) {
  bases <- col[col %in% c("A","C","G","T")]
  length(unique(bases)) > 1
}))
snp_matrix <- mat[, snp_positions]
snp_list <- lapply(snp_positions, function(pos) {
  data.frame(
    Sequence = rownames(mat),
    Position = pos,
    Base = mat[, pos]
  )
})
names(snp_list) <- paste0("Pos_", snp_positions)

#Alignment

The final skill we need to learn is how to align .fasta files to a reference genome. The data we have today is aligned so we’ll have to arbitrarily “unalign” it. We will use our first sequence as the reference.

UK_seqs_nogap <- DNAStringSet(gsub("-", "", as.character(UK_seqs))) 
ref <- seqs[6]

The region of interest from hcmv will be our query site.

library(Biostrings)
query <- readDNAStringSet("hcmv15-hcmv1b-1.fasta")
reversecompquery <- reverseComplement(query)
ref <- seqs["Europe_Germany_JX512204.1_Human_herpesvirus_5_strain_HAN16_complete_genome"]
#ref <- readDNAStringSet("hcmvgenome.fasta")
query <- DNAString(toupper(as.character(reversecompquery)))
#ref <- DNAStringSet(gsub("-", "", as.character(ref))) 

hits <- matchPattern(query, ref[[1]], max.mismatch=20)

We can now look at our aligned data!

hits

Homework Question 1: Interpretting QC data.

  1. Look at the ‘Cycle-specific Base Call Proportion’ rqc file. What is the ratio of A/C/G/T in cycles 1-10. What is the ratio of A/C/G/T in cycles 65-75? Based on your knowledge of DNA-sequencing, propose a hypothesis for why the cycles have a different ratio.

  2. Look at the ‘Cycle specific Quality Distribution’. What is the typical quality of cycles 1-10? What is the typical quality of cycles 65-75? Based on your knowledge of DNA-sequencing, propose a hypothesis for why the cycles have a different quality

Homework Question 2: Loading a File, Reverse Complements

For one of the research projects, we will be studying a specific protein coding site of hcmv associated with the transition of hcmv from latent to active. Align the following two queries to any genome in the reference file. Which one successfully aligns? What does that mean about the site of interest?

query <- readDNAStringSet("hcmv15-hcmv1b-1.fasta")
reversecompquery <- reverseComplement(query)

Homework Question 3: Isolating a Site

Isolate the ‘hcmv15-hcmv1b-1.fasta’ site in all 259 seqs. Your final result should be a DNAString set of 259 sequences about 4-5 kb long.