如何使用 bash 脚本查找 fasta 文件的 gc 内容？

Question

我正在学习生物信息学。我想使用 Bash 脚本从 fasta 文件中查找 GC 内容。 GC 内容基本上是（（g + c）的数量/（（a + t + g + c）的数量）。感谢任何帮助。

我正在尝试使用 wc 命令。但我无法得到答案。

Answer 1

这应该有效：

#!/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
' "$@" >&2
  exit 1
}

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

# ================
# 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/(^>.*$)/##/' "$FASTA_FILE" 
      | tr --delete "
" 
      | sed 's/$/#/' 
      | tr "#" "
" 
      | 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
' "$_line" | grep --quiet '^>'; then
      # Save sequence and continue
      _sequence=${_line#?}
      continue
    fi

    # Count
    _a_count=$(printf '%s
' "$_line" | tr --delete --complement 'A' | wc --bytes)
    _c_count=$(printf '%s
' "$_line" | tr --delete --complement 'C' | wc --bytes)
    _g_count=$(printf '%s
' "$_line" | tr --delete --complement 'G' | wc --bytes)
    _t_count=$(printf '%s
' "$_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_DIGITS" "$_a_count" "$_c_count" "$_g_count" "$_t_count" 
        | bc
    )
    # Add 0 before decimal point
    _gc="$(printf "%.${N_DIGITS}f
" "$_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_DIGITS" "$_a_count_total" "$_c_count_total" "$_g_count_total" "$_t_count_total" 
      | bc
  )
  # Add 0 before decimal point
  _gc="$(printf "%.${N_DIGITS}f
" "$_gc")"
  info "GC content: $_gc"
}

这 ”计算序列数“ 和 ”删除序列包装"代码改编自https://www.biostars.org/p/17680

该脚本仅使用除bc 之外的基本命令来进行精度计算（参见bc installation）。

您可以通过修改CONFIGURATION 部分中的变量来配置脚本。

因为您没有指明您想要哪一个，所以会计算每个序列和整体的 GC 含量。因此，摆脱任何不必要的东西:)

尽管我缺乏生物信息学背景，但该脚本成功地解析和分析了一个 fasta 文件。