How can I find the GC-content of a FASTA file using a Bash script?

I want to find the GC-content from a FASTA format format file using a Bash script.
The GC-content is basically (number of (G + C)) / (number of (A + T + G + C)).

I am trying to use the wc command. But I was not able to get an answer.

After going through documentation and videos, I came up with a solution.

filename=$@                     # Collecting all the filenames as parameters
for f in $filename              # Looping over files
do
    echo " $f is being processed..."
    gc=( $( grep -v ">" < "$f" | grep -io 'g\|c'< "$f" | wc -l))    # Reading lines that don’t start with < using -v. grep -io matches to either g or c and outputs each match on single line. wc -l counts the number of lines or indirectly the number of g and c. This is stored in a variable.
    total=( $( grep -v ">" < "$f" | tr -d '\s\r' | wc -c))          # Spaces, tabs, new line are removed from the file using tr. Then the number of characters are counted by wc -c
    percent=( $( echo "scale=2;100*$gc/$total" |bc -l))             # bc -l is used to get the answer in float format. scale=2 mentions the number of decimal points.
    echo " The GC content of $f is: "$percent"%"
    echo
done

I am learning bioinformatics.

Do not reinvent the wheel. For common bioinformatics tasks, use open-source tools that are specifically designed for these tasks, are well-tested, widely used, and handle edge cases. For example, use EMBOSS infoseq utility. EMBOSS can be easily installed, for example using conda.

Example:

Install EMBOSS package (do once):

conda create --name emboss emboss --channel iuc

Activate the conda environment and use EMBOSS infoseq, here to print the sequence name, length and percent GC:

source activate emboss

cat your_sequence_file_name.fasta | infoseq -auto -only -name -length -pgc stdin

source deactivate

This prints into standard output something like this:

Name           Length %GC
seq_foo        119    60.50
seq_bar        104    39.42
seq_baz        191    46.60
...

This should work:

#!/usr/bin/env sh
# Adapted from https://www.biostars.org/p/17680

# Fail on error
set -o errexit
# Disable undefined variable reference
set -o nounset

# ================
# CONFIGURATION
# ================
# Fasta file path
FASTA_FILE="file.fasta"
# Number of digits after decimal point
N_DIGITS=3

# ================
# LOGGER
# ================
# Fatal log message
fatal() {
  printf '[FATAL] %s\n' "$@" >&2
  exit 1
}

# Info log message
info() {
  printf '[INFO ] %s\n' "$@"
}

# ================
# MAIN
# ================
{
  # Check command 'bc' exist
  command -v bc > /dev/null 2>&1 || fatal "Command 'bc' not found"
  # Check file exist
  [ -f "$FASTA_FILE" ] || fatal "File '$FASTA_FILE' not found"

  # Count number of sequences
  _n_sequences=$(grep --count '^>' "$FASTA_FILE")
  info "Analyzing $_n_sequences sequences"
  [ "$_n_sequences" -ne 0 ] || fatal "No sequences found"

  # Remove sequence wrapping
  _fasta_file_content=$(
    sed 's/\(^>.*$\)/#\1#/' "$FASTA_FILE" \
      | tr --delete "\r\n" \
      | sed 's/$/#/' \
      | tr "#" "\n" \
      | sed '/^$/d'
  )

  # Vars
  _sequence=
  _a_count_total=0
  _c_count_total=0
  _g_count_total=0
  _t_count_total=0

  # Read line by line
  while IFS= read -r _line; do
    # Check if header
    if printf '%s\n' "$_line" | grep --quiet '^>'; then
      # Save sequence and continue
      _sequence=${_line#?}
      continue
    fi

    # Count
    _a_count=$(printf '%s\n' "$_line" | tr --delete --complement 'A' | wc --bytes)
    _c_count=$(printf '%s\n' "$_line" | tr --delete --complement 'C' | wc --bytes)
    _g_count=$(printf '%s\n' "$_line" | tr --delete --complement 'G' | wc --bytes)
    _t_count=$(printf '%s\n' "$_line" | tr --delete --complement 'T' | wc --bytes)

    # Add current count to total
    _a_count_total=$((_a_count_total + _a_count))
    _c_count_total=$((_c_count_total + _c_count))
    _g_count_total=$((_g_count_total + _g_count))
    _t_count_total=$((_t_count_total + _t_count))

    # Calculate GC content
    _gc=$(
      printf 'scale = %d; a = %d; c = %d; g = %d; t = %d; (g + c) / (a + c + g + t)\n' \
        "$N_DIGITS" "$_a_count" "$_c_count" "$_g_count" "$_t_count" \
        | bc
    )
    # Add 0 before decimal point
    _gc="$(printf "%.${N_DIGITS}f\n" "$_gc")"

    info "Sequence '$_sequence' GC content: $_gc"
  done << EOF
$_fasta_file_content
EOF

  # Total data
  info "Adenine total count: $_a_count_total"
  info "Cytosine total count: $_c_count_total"
  info "Guanine total count: $_g_count_total"
  info "Thymine total count: $_t_count_total"

  # Calculate total GC content
  _gc=$(
    printf 'scale = %d; a = %d; c = %d; g = %d; t = %d; (g + c) / (a + c + g + t)\n' \
      "$N_DIGITS" "$_a_count_total" "$_c_count_total" "$_g_count_total" "$_t_count_total" \
      | bc
  )
  # Add 0 before decimal point
  _gc="$(printf "%.${N_DIGITS}f\n" "$_gc")"
  info "GC content: $_gc"
}

The "Count number of sequences" and "Remove sequence wrapping" codes are adapted from https://www.biostars.org/p/17680

The script uses only basic commands except for bc to do the precision calculation (See bc installation).

You can configure the script by modifying the variables in the CONFIGURATION section.

Because you haven't indicated which one you want, the GC content is calculated for both each sequence and the overall. Therefore, get rid of anything that isn't necessary

Despite my lack of bioinformatics background, the script successfully parses and analyzes a fasta file.