盛夏不等式: Computing GC Content (ID: GC)

Problem

給定至多含 10 條 DNA 序列之 fasta 檔，求 GC 比最高者及其 GC 比。

Given: At most 10 DNA strings in FASTA format (of length at most 1 kbp each).

Return: The ID of the string having the highest GC-content, followed by the GC-content of that string. Rosalind allows for a default error of 0.001 in all decimal answers unless otherwise stated; please see the note on absolute error below.

Background

GC-content

GC-content 是核酸中，Guanine（鳥糞嘌呤）和 Cytosine（胞嘧啶）的占比。由於鹼基層疊和氫鍵數量等因素， A-T 配對不若 G-C 配對穩定 [1]。因此，GC-content 可以粗略地反映核酸的穩定性並且應用於以下場合，

分子演化學：由於 (1) 鹼基突變為彼此的機率不一 (Substitution Biases)，而且 (2) GC-content 與序列的穩定性有關，序列的 GC-content 可能反映了該物種演化的資訊。
評估樣本汙染：依照經驗法則，同一物種的基因序列之 GC-content 呈鐘形分布 [2]，若定序資料中 GC-content 呈多峰分布，可能意味著樣本中含有其他物種的序列汙染，或是存在某些大量表現的基因影響了分布。
預測黏合溫度（annealing temperature）：黏合溫度是 P１CR 時能讓 primer 與目標序列黏合的溫度。由於 primer 與標的序列配對的穩定度越高，黏合溫度也越高，所以 GC-content 可作為推估黏合溫度的參數。

FASTA format

FASTA format 用於紀錄核酸或胜肽的序列。在 FASTA format 中，一條序列的紀錄包含 (1) 以 ">" 起頭的描述列，紀錄序列的名稱、辨識碼或來源等資訊。(2) 由鹼基或胺基酸字符組成的序列字串（鹼基和胺基酸字符的 IUPAC 命名，可參考維基百科的簡介）。典型的 FASTA format 如下，由於每列的字串長度依慣例不超過 80 個字元（這可能跟方便在終端介面閱讀有關），資料處理時要留意較長的序列往往被換行符分隔。

> OTU1
ATCGATCGATCGATCGATCAGATATATGATGATAGATCCCCATGATAGATCGATCGATCGATCGATCAGATATATGATGA
ATCGACGT
> OTU2
ATCGATCGGGCGATCGATCAGATATATGACGCGAAATCCCCATGATAG

以 FASTA format 儲存的文字檔案常有 .fa、.fna、.fasta 等副檔名。

Solution

此題可分為三個部分：讀取 fasta 檔案、計算 GC 比、求 GC 比最高的序列，對應的資料處技巧分別為字串擷取、字符統計和取極值。其中，讀取 fasta 檔用於任何需要處理序列的步驟，例如篩選 OTUs/ASVs 或註解分類資訊。由於 fasta 檔案的描述列以 ">" 開頭，可藉由 ">" 符號擷取序列編號及描述，再拼接夾在描述列之間的字串，以取得完整的序列。

取得字串後，可應用 Counting 所學的技術，計算序列當中 G 和 C 的數量之和，再除以序列的長度（即鹼基總數）即可獲得 GC 比。最後透過連續的比較，取得 GC 比最大的序列及其 GC 比。

python

在 python 中可用 count() 計算序列中 G 和 C 的數量，再用 len 統計序列長度。隨後計算 GC 比例時要留意，在 python 2 當中，若除數與被除數皆為整數，除法運算子 "/" 預設會無條件捨去結果的小數（即把結果視同整數）。因此，若要取得含小數的近似解，需先將除數與被除數其一轉變為浮點數再計算。不過除法的這項特性已在 python 3 更動，若使用 "/" 會回傳具浮點數的近似解，而使用 "//" 才會獲得無條件捨去小數的結果。至於求餘數或長除法等其他除法功能，可參考比較 Python 與除法相關的運算子與函式 – /、//、% 與 divmod

def count_GC (seq):
    return 100 * float(seq.count("G") + seq.count("C")) / len(seq)

dictionary 的 key 和 value 可用於記錄序列的名稱和字串。首先以 open() 開啟 fasta 檔案，接著建立空的 dictionary，再一列一列讀取 fasta 檔案的內容。接著，凡是碰到以 ">" 開頭者，即擷取 ">" 與 "\n"（換行號）之間的字串作為序列辨識碼。由於鹼基字串位於 fasta 檔以 ">" 開頭的描述列之間，所以合併讀取到下一個 ">" 之前的內容，即可獲得該序列的鹼基字串。

此處用到的 funtion 分別為，

startswith(str)：判斷字串是否以指定字符串為首，此處用於判斷讀取的列是否為描述列
lstrip(str)：從字串左側移除指定字符串，此處用於移除描述列的 ">" 符號
replace(old, new): 以 new（新字符串）取代字串中的 old（舊字符串），此處用於移除 "\n"

# Input the path to fasta file
def read_fasta (path):
    f = open(path,  "r")
    seqs = {}
    for line in f:
        if line.startswith(">"):
            ID = line.lstrip(">").replace("\n", "")
            seqs[ID] = "" 
        else:
            seqs[ID] = seqs[ID] + line.replace("\n", "")       
    return seqs

求 GC 比最大值的方法是先假設最大的 GC 比為 0，接著比較各序列的 GC 比是否大於假設的最大 GC 比。若是，則更新假設；若否，則維持假設。一旦所有序列的 GC 比都比較過，即可獲得此檔案的最大 GC 比。

seqs = read_fasta("seqs.fasta")
maxGC = 0
for s in seqs:
    if maxGC < count_GC(seqs[s]):
        ID = s
        maxGC = count_GC(seqs[s])     
print ID
print maxGC

R

R 有其獨特的向量式運算，亦即諸如四則運算或字符查找等 function 可直接套用在向量的每個項目。由於向量運算的效率往往高於迴圈，所以此處擷取 seqinr 套件的 read.fasta() 的程式碼（seqinr 是專門用於處理序列資料的套件），嘗試以向量式運算求出 GC 比最大的序列。

首先，以 readLines() 讀取 fasta 檔案的所有內容並存為字串向量。接著，使用 grep() 搜尋以 ">" 開頭的描述列之索引號，並以 gsub 移除 ">" 取回序列辨識碼。最後利用 paste() 連接位於描述列的索引號之間的鹼基字串。

read_fasta <- function (path) {
  lines <- readLines(path)
  ind <- grep(">", lines)
  start <- ind + 1
  end <- c((ind - 1)[-1], length(lines))
  seqs <- lapply(1:length(ind), function(i) paste(lines[start[i]:end[i]], collapse = ""))
  names(seqs) <- gsub(">", "", lines[ind])
  return(seqs)
}

而計算 GC 比方面，我則模仿 python 的 count()，建立了計算字串中給定字符出現頻率的 function。簡言之，透過 gregexpr() 識別給定字符在字串中出現的位置，再計算回傳位置的總數得出字符出現頻率。由於 gregexpr() 以 list() 回傳結果，所以要以 lengths() 而非 length() 計算回傳值的長度。若要使用 length()，需要以 unlist() 把 list 型式的回傳結果轉變為向量再計算。

length: 計算向量長度
lengths: 計算 list 中，每個元素的長度

count <- function (str, char) {
  return(lengths(gregexpr(char, str)))
  # alternative
  # return(length(unlist(gregexpr(char, str))))
}
count_GC <- function (seq) {
  return(100 * (count(seq, "G") + count(seq, "C")) / nchar(seq))
}

最後計算不同序列的 GC 比，以 which.max() 找出 GC 比最大者的索引值取出答案。

seqs <- read_fasta("./test.fasta")
GC <- lapply(seqs, count_GC)
print(GC[which.max(GC)])

References

Yakovchuk, Protozanova, and Frank-Kamenetskii. (2006). Base-stacking and base-pairing contributions into thermal stability of the DNA double helix. Nucleic acids research, 34(2), 564-574.
Romiguier et al. (2010). Contrasting GC-content dynamics across 33 mammalian genomes: relationship with life-history traits and chromosome sizes. Genome research, 20(8), 1001-1009.

盛夏不等式

Computing GC Content (ID: GC)