Parsing Information Out of Sequencing Data

Sequencing data is a critical component in various fields such as genetics, genomics, and molecular biology. The raw sequencing data consists of a series of nucleotide sequences that are read in order to determine the genetic material of an organism. However, this raw data can be overwhelming and difficult to analyze manually.

One common approach to manage and analyze large amounts of sequencing data is by converting it into a text format. This allows for easier manipulation and analysis using computational tools. In this article, we will discuss how to parse information out of sequencing data in R, specifically focusing on calculating nucleotide diversity at each position.

Introduction

Nucleotide diversity is a measure of genetic variation within a population or species. It can be used to understand the evolutionary history of an organism and its adaptability to different environments. Calculating nucleotide diversity involves counting the frequency of each nucleotide (A, C, G, T) at each position in a sequence.

In this article, we will explore how to calculate nucleotide diversity using R and the seqinr package.

Prerequisites

To follow along with this article, you should have:

• A basic understanding of R programming language
• The seqinr package installed and loaded in your R environment
• A fasta file containing the sequences you want to analyze (you can create one using a tool like FastQC)

Reading Fasta Files

To start, we need to read our fasta file into R. The seqinr package provides a function called read.fasta() that allows us to do this.

# Load the necessary packages
library(tidyverse)
library(seqinr)

# Read the fasta file
fasta_files <- list.files("/path/to/your/directory", pattern = ".fa", full.names = TRUE)
seqs <- lapply(fasta_files, function(x){
  seqs <- read.fasta(x)
})

Calculating Nucleotide Diversity

To calculate the nucleotide diversity at each position, we can use two functions from the seqinr package: getFrag() and count(). The getFrag() function allows us to subset a sequence into smaller fragments based on a start and end position. The count() function then calculates the frequency of each nucleotide within these fragments.

# Obtain fragment lengths - assuming all sequences are the same length!
l <- length(seqs[[1]])

# Use lapply to estimate frequency for each position
p <- lapply(1:l, function(i){
  # Obtain the nucleotide for the current position
  pos_seq <- getFrag(seqs, i, i)
  
  # Get the frequency of each nucleotide
  pos_freq <- count(pos_seq, wordsize = 1, freq = TRUE)
  
  # Convert to data.frame and add information about the nucleotide position
  pos_freq <- as.data.frame() %>%
    rename(nuc = Var1, freq = Freq) %>%
    mutate(pos = i)
}, seqs = seqs)

# Bind all tables together
diversity <- bind_rows(p) %>% 
  filter(freq > 0) %>% 
  group_by(pos) %>% 
  summarise(shannon_entropy = -sum(freq * log2(freq)),
            het = 1 - sum(freq^2), n_nuc = n())

Visualizing the Results

Once we have calculated the nucleotide diversity for each position, we can visualize our results using ggplot2. This allows us to create a plot that shows how the diversity changes across different positions.

# Create a scatterplot of Shannon Entropy vs. Position (bp)
ggplot(diversity, aes(x = pos, y = shannon_entropy)) +
  geom_line() +
  geom_point(aes(colour = factor(n_nuc))) +
  labs(x = "Position (bp)", y = "Shannon Entropy", 
       colour = "Number of\nnucleotides")

Applying the Code to Multiple Fasta Files

If you have multiple fasta files that you want to analyze, you can use the lapply() function in R to apply the code above to each file.

# Find all the fasta files of interest
fasta_files <- list.files("/path/to/your/directory", 
                          pattern = ".fa", full.names = TRUE)

# Use lapply to apply the code above to each file
my_diversities <- lapply(fasta_files, function(f){
  # Read the fasta file
  seqs <- read.fasta(f)
  
  # Obtain fragment lengths - assuming all sequences are the same length!
  l <- length(seqs[[1]])
  
  # Use the `lapply` function to estimate frequency for each position
  p <- lapply(1:l, function(i){
    # Obtain the nucleotide for the current position
    pos_seq <- getFrag(seqs, i, i)
    
    # Get the frequency of each nucleotide
    pos_freq <- count(pos_seq, wordsize = 1, freq = TRUE)
    
    # Convert to data.frame and add information about the nucleotide position
    pos_freq <- as.data.frame() %>%
      rename(nuc = Var1, freq = Freq) %>%
      mutate(pos = i)
  }, seqs = seqs)
  
  # Bind all tables together
  diversity <- bind_rows(p) %>% 
    filter(freq > 0) %>% 
    group_by(pos) %>% 
    summarise(shannon_entropy = -sum(freq * log2(freq)),
              het = 1 - sum(freq^2), n_nuc = n())
})

# Name the list elements and bind them together
names(my_diversities) <- basename(fasta_files)
my_diversities <- bind_rows(my_diversities, .id = "file_name")

This will give you a table of diversities for each file. You can then use ggplot2 to visualize it, similarly to what I did above, but perhaps using facets to separate the diversity from each file into different panels.

By following this article, you should now have a better understanding of how to parse information out of sequencing data in R and calculate nucleotide diversity at each position.

Last modified on 2024-03-27