Function to calculate the GC content variation in a sequence

Question

I came across this BMC Genomics paper: Analysis of intra-genomic GC content homogeneity within prokaryotes

And I implemented some Python functions to make this available as part of a personal project. The functions are these:

import numpy as np
from collections import defaultdict
def get_sequence_chunks(sequence, window, step, as_overlap=False):
 """GC Content in a DNA/RNA sub-sequence length k. In
 overlap windows of length k, other wise non overlapping.
 
 Inputs:
 sequence - a string representing a DNA sequence. 
 as_overlap - boolean that represents if overlap is needed.
 k - a integer representing the lengths of overlapping bases.
 Default is 20.
 
 Outputs:
 subseq - a string representing a slice of the sequence with or without
 overlapping characters.
 """
 # make sequence upper case and getting the length of it
 sequence, seq_len = sequence.upper(), len(sequence)
 # non overlap sequence length
 non_overlap = range(0, seq_len - window + 1, step)
 # overlap sequence length
 overlap = range(0, seq_len - window + 1)
 # overlap is needed
 if as_overlap:
 # iterates to the overlap region
 for i in overlap:
 # creates the substring 
 subseq = sequence[i:i + window]
 yield subseq
 # if non overlap is choosed
 else:
 # iterates to the non overlap region
 for j in non_overlap:
 # creates the sub-string to count the gc_content
 subseq = sequence[j:j + window]
 yield subseq
def gc_content(seq, percent=True):
 """
 Calculate G+C content, return percentage (as float between 0 and 100).
 Copes mixed case sequences, and with the ambiguous nucleotide S (G or C)
 when counting the G and C content. The percentage is calculated against
 the full length.
 
 Inputs:
 sequence - a string representing a sequence (DNA/RNA/Protein)
 
 Outputs:
 gc - a float number or a percent value representing the gc content of the sequence. 
 
 """
 seq_len = len(seq) 
 gc = sum(seq.count(x) for x in ["G", "C", "g", "c", "S", "s"])
 if percent:
 try:
 return gc * 100.0 / seq_len
 except ZeroDivisionError:
 return 0.0 
 else:
 return gc / seq_len 
def get_gc_content_slide_window(sequence, window, step):
 """
 Calculate teh GC (G+C) content along a sequence. Returns a dictionary -like
 object mapping the sequence window to the gc value (floats).
 Returns 0 for windows without any G/C by handling zero division errors, and 
 does NOT look at any ambiguous nucleotides.
 
 Inputs:
 sequence - a string representing a sequence (DNA/RNA/Protein)
 window_size - a integer representing the length of sequence to be analyzed
 at each step in the full sequence.
 step - a integer representing the length of oerlapping sequence allowed.
 
 Outputs:
 gc- a dictionary-like object mapping the window size sequence to it gc values 
 as floating numbers
 """
 gc = defaultdict(float)
 for s in get_sequence_chunks(sequence, window, step, as_overlap=False):
 gc[s] = gc.get(s, 0.0) + gc_content(s)
 return gc 
def difference_gc(mean_gc, gc_dict):
 """
 Calculates the difference between the mean GC content of window i, and the
 mean chromosomal GC content, as Di = GC i − GC
 """
 # get the keys and default values in the dictionary
 d_i = defaultdict(float, [(k, 0.0) for k in gc_dict.keys()])
 # iterate through all keys, gc calculated values for the
 # genome chunk
 for chunk, gc_cnt in gc_dict.items():
 # get the difference between the chromosomal mean and the chunks
 dif = gc_cnt - mean_gc
 # add the difference to the appropriate chunk
 d_i[chunk] = dif
 return d_i
def get_chromosomal_gc_variantion(difference_dic):
 """
 Calculates the chromosomal GC variation defined as the log-transformed average 
 of the absolute value of the difference between the mean GC content of each 
 non-overlapping sliding window i and mean chromosomal GC content. 
 chromosomal_gc_variantion = log(1/N*sum(|Di|), where N is the maximum number of 
 non-overlapping window size in bp in a sliding windows.
 """
 # get the number of chunks
 n = len(difference_dic.keys())
 arr = np.array(list(difference_dic.values()))
 var = np.log((1/n) * sum(abs(arr)))
 return var

I worked with e. coli genome: Escherichia coli strain RM13322 chromosome, complete genome

gc = gc_content(ec, percent=True) = 50.76819424506625

It is similar to found here: Genome Assembly and Annotation report (26286) ~ 50.73.

The slide window gave it non overlapping results like, where ec = escherichia coli genome, 100 size window, and step = 100:

get_gc_content_slide_window(ec, 100, 100)
defaultdict(float,
 {'GTTGGCATGCCATGGTGATGTTTGTTAGGAAAGCAAAGATGGCAAAACTGCTGGGGGTTTTGTGGTTGAGTATGCCAATATAATTAATAGATTAAAGAGT': 39.0,
 'TAGTTGTGAAGAAAATATGGATAAACAGGACGACGAATGCTTTCACCGATAAGGACAACTTTCCATAACTCAGTAAATATAGTGCAGAGTTCACCCTGTC': 39.0,
 'AAACGGTTTCTTTTGCAGGAAAGGAATATGAGTTAAAGGTCATTGATGAAAAAACGCCTATTCTTTTTCAGTGGTTTGAACCTAATCCTGAACGATATAA': 34.0,
 'GAAAGATGAGGTTCCAATAGTTAATACTAAGCAGCATCCCTATTTAGATAATGTCACAAATGCGGCAAGGATAGAGAGTGATCGTATGATAGGTATTTTT': 36.0,
 'GTTGATGGCGATTTTTCAGTCAACCAAAAGACTGCTTTTTCAAAATTGGAACGAGATTTTGAAAATGTAATGATAATCTATCGGGAAGATGTTGACTTCA': 34.0,
 'GTATGTATGACAGAAAACTATCAGATATTTATCATGATATTATATGTGAACAAAGGTTACGAACTGAAGACAAAAGAGATGAATACTTGTTGAATCTGTT': 29.0,...}

The difference function (slide window gc content - gc content all sequence), where ec_gc = gc content of e.coli genome and ec_non = gc content in each subsequence of size 100 non overlapped:

difference_gc(ec_gc, ec_non)
defaultdict(float,
 {'GTTGGCATGCCATGGTGATGTTTGTTAGGAAAGCAAAGATGGCAAAACTGCTGGGGGTTTTGTGGTTGAGTATGCCAATATAATTAATAGATTAAAGAGT': -11.768194245066248,
 'TAGTTGTGAAGAAAATATGGATAAACAGGACGACGAATGCTTTCACCGATAAGGACAACTTTCCATAACTCAGTAAATATAGTGCAGAGTTCACCCTGTC': -11.768194245066248,
 'AAACGGTTTCTTTTGCAGGAAAGGAATATGAGTTAAAGGTCATTGATGAAAAAACGCCTATTCTTTTTCAGTGGTTTGAACCTAATCCTGAACGATATAA': -16.768194245066248,
 'GAAAGATGAGGTTCCAATAGTTAATACTAAGCAGCATCCCTATTTAGATAATGTCACAAATGCGGCAAGGATAGAGAGTGATCGTATGATAGGTATTTTT': -14.768194245066248,
 'GTTGATGGCGATTTTTCAGTCAACCAAAAGACTGCTTTTTCAAAATTGGAACGAGATTTTGAAAATGTAATGATAATCTATCGGGAAGATGTTGACTTCA': -16.768194245066248,...}

And the (ec_dif = difference on gc content from total sequence and from subsequences):

get_chromosomal_gc_variantion(ec_dif)
1.7809591767493207

My principal concern It is to know if this function get_chromosomal_gc_variantion(ec_dif) was implemented as in the paper above (GCVAR):

journal excerpt

As suggested, a toy example:

s="TAGTTGTGAAGAAAATATGGATAAACAGGACGACGAATGCTTTCACCGATAAGGACAACTTTCCATAACTCAGTAAATATAGTGCAGAGTTCACCCTGTC"
seq_gc_non_overllap = get_gc_content_slide_window(s, 20, 20)
seq_gc_non_overllap
defaultdict(float,
 {'TAGTTGTGAAGAAAATATGG': 30.0,
 'ATAAACAGGACGACGAATGC': 45.0,
 'TTTCACCGATAAGGACAACT': 40.0,
 'TTCCATAACTCAGTAAATAT': 25.0,
 'AGTGCAGAGTTCACCCTGTC': 55.0})
seq_gc_total = gc_content(s, percent=True)
seq_gc_total
39.0
seq_dif_gctot_gc_slidewindow = difference_gc(seq_gc_total, seq_gc_non_overllap)
defaultdict(float,
 {'TAGTTGTGAAGAAAATATGG': -9.0,
 'ATAAACAGGACGACGAATGC': 6.0,
 'TTTCACCGATAAGGACAACT': 1.0,
 'TTCCATAACTCAGTAAATAT': -14.0,
 'AGTGCAGAGTTCACCCTGTC': 16.0})
chromosome_gc_variation = get_chromosomal_gc_variantion(seq_dif_gctot, gc_slidewindow)
chromosome_gc_variation
2.2192034840549946

Any tip to improve the code would be very welcome, once I only have been working with python for a couple of years and I am not very skilled, and much less with numpy.

• Reinderien: What are ec, ec_gc, ec_non, ec_dif?
• Reinderien: It's good for you to have described them, but it would be more helpful to show the code that loads them, and perhaps a reduced-size sample for your reviewers to be able to run your code.
• kubatucka: @PauloSergioSchlogl Could you also link your project once you've completed/published it?

Answer 1

Quite interesting!

You will benefit from introducing PEP484 type hints.

Your docstring for get_sequence_chunks is way out of whack; you'll need to revisit your parameter names.

There isn't a whole lot of benefit from tuple-combination assignment for your sequence and seq_len so I suggest just doing them separately.

Your overlap logic can be greatly simplified: notice that the only difference between the two code paths is the step size.

Minor spelling issues such as other wise that should be one word, overllap -> overlap, and variantion -> variation.

Comments such as

# make sequence upper case and getting the length of it

are less helpful than not having a comment at all. Note the difference between that and

# get the difference between the chromosomal mean and the chunks

which has content that can't necessarily be inferred from the code alone (though even this will be redundant with a sufficiently descriptive variable name).

This comprehension:

gc = sum(seq.count(x) for x in ["G", "C", "g", "c", "S", "s"])

will iterate through the entire seq string for every potential character. It's likely that refactoring this to use set membership checks and one single iteration through the string will be more efficient.

This code to prevent division-by-zero:

seq_len = len(seq) 
gc = sum(seq.count(x) for x in ["G", "C", "g", "c", "S", "s"])
if percent:
 try:
 return gc * 100.0 / seq_len
 except ZeroDivisionError:
 return 0.0 
else:
 return gc / seq_len 

first of all can trivially avoid logic-by-exception by checking if the sequence length is zero; and also should extend that check to the case where a non-percent figure is returned.

In difference_gc the initial defaultdict initialization is not needed, since you overwrite every single key. You can replace the whole function with one dictionary comprehension.

Try to avoid abbreviating words like count to cnt.

The len of a dict's keys is equal to the len of a dict itself, so you can drop the .keys() in len(difference_dic.keys()).

You can avoid the construction of a list here:

arr = np.array(list(difference_dic.values()))

by using fromiter instead.

Your example is kind of broken. get_chromosomal_gc_variation only accepts one parameter.

Consider reworking your example into a basic unit test.

Suggested

import math
from typing import Iterator, DefaultDict, Dict
import numpy as np
from collections import defaultdict
def get_sequence_chunks(
 sequence: str,
 window: int,
 step: int,
 as_overlap: bool = False,
) -> Iterator[str]:
 """GC Content in a DNA/RNA sub-sequence length k. In
 overlap windows of length k, other wise non overlapping.
 Inputs:
 sequence - a string representing a DNA sequence.
 as_overlap - boolean that represents if overlap is needed.
 k - a integer representing the lengths of overlapping bases.
 Default is 20.
 Outputs:
 subseq - a string representing a slice of the sequence with or without
 overlapping characters.
 """
 sequence = sequence.upper()
 if as_overlap:
 # overlap sequence length
 step = 1
 # iterates to the overlap region
 for i in range(0, len(sequence) - window + 1, step):
 # creates the substring
 yield sequence[i:i + window]
def gc_content(seq: str, percent: bool = True) -> float:
 """
 Calculate G+C content, return percentage (as float between 0 and 100).
 Copes mixed case sequences, and with the ambiguous nucleotide S (G or C)
 when counting the G and C content. The percentage is calculated against
 the full length.
 Inputs:
 sequence - a string representing a sequence (DNA/RNA/Protein)
 Outputs:
 gc - a float number or a percent value representing the gc content of the sequence.
 """
 seq_len = len(seq)
 gc = sum(1 for x in seq.upper() if x in {"G", "C", "S"})
 if seq_len == 0:
 return 0
 gc /= seq_len
 if percent:
 return gc * 100
 return gc
def get_gc_content_slide_window(sequence: str, window: int, step: int) -> DefaultDict[str, float]:
 """
 Calculate the GC (G+C) content along a sequence. Returns a dictionary-like
 object mapping the sequence window to the gc value (floats).
 Returns 0 for windows without any G/C by handling zero division errors, and
 does NOT look at any ambiguous nucleotides.
 Inputs:
 sequence - a string representing a sequence (DNA/RNA/Protein)
 window_size - a integer representing the length of sequence to be analyzed
 at each step in the full sequence.
 step - a integer representing the length of oerlapping sequence allowed.
 Outputs:
 gc- a dictionary-like object mapping the window size sequence to it gc values
 as floating numbers
 """
 gc = defaultdict(float)
 for s in get_sequence_chunks(sequence, window, step, as_overlap=False):
 gc[s] += gc_content(s)
 return gc
def difference_gc(mean_gc: float, gc_dict: Dict[str, float]) -> Dict[str, float]:
 """
 Calculates the difference between the mean GC content of window i, and the
 mean chromosomal GC content, as Di = GC i − GC
 """
 # iterate through all keys, gc calculated values for the
 # genome chunk
 # get the difference between the chromosomal mean and the chunks
 # add the difference to the appropriate chunk
 d_i = {
 chunk: gc_cnt - mean_gc
 for chunk, gc_cnt in gc_dict.items()
 }
 return d_i
def get_chromosomal_gc_variation(difference_dic: Dict[str, float]) -> float:
 """
 Calculates the chromosomal GC variation defined as the log-transformed average
 of the absolute value of the difference between the mean GC content of each
 non-overlapping sliding window i and mean chromosomal GC content.
 chromosomal_gc_variantion = log(1/N*sum(|Di|), where N is the maximum number of
 non-overlapping window size in bp in a sliding windows.
 """
 # get the number of chunks
 n = len(difference_dic)
 arr = np.fromiter(
 difference_dic.values(),
 dtype=np.float32,
 count=len(difference_dic),
 )
 var = np.log(np.sum(np.abs(arr)) / n)
 return var
def test() -> None:
 s = (
 "TAGTTGTGAAGAAAATATGGATAAACAGGACGACGAATGCTTTCACCGAT"
 "AAGGACAACTTTCCATAACTCAGTAAATATAGTGCAGAGTTCACCCTGTC"
 )
 seq_gc_non_overlap = get_gc_content_slide_window(s, 20, 20)
 expected = {
 'TAGTTGTGAAGAAAATATGG': 30,
 'ATAAACAGGACGACGAATGC': 45,
 'TTTCACCGATAAGGACAACT': 40,
 'TTCCATAACTCAGTAAATAT': 25,
 'AGTGCAGAGTTCACCCTGTC': 55,
 }
 assert expected.keys() == seq_gc_non_overlap.keys()
 for k, v in expected.items():
 assert math.isclose(seq_gc_non_overlap[k], v)
 seq_gc_total = gc_content(s, percent=True)
 assert math.isclose(seq_gc_total, 39)
 seq_dif_gctot_gc_slidewindow = difference_gc(seq_gc_total, seq_gc_non_overlap)
 expected = {
 'TAGTTGTGAAGAAAATATGG': -9,
 'ATAAACAGGACGACGAATGC': 6,
 'TTTCACCGATAAGGACAACT': 1,
 'TTCCATAACTCAGTAAATAT': -14,
 'AGTGCAGAGTTCACCCTGTC': 16,
 }
 assert expected.keys() == seq_dif_gctot_gc_slidewindow.keys()
 for k, v in expected.items():
 assert math.isclose(seq_dif_gctot_gc_slidewindow[k], v)
 chromosome_gc_variation = get_chromosomal_gc_variation(seq_dif_gctot_gc_slidewindow)
 assert math.isclose(chromosome_gc_variation, 2.2192034840549946)
if __name__ == '__main__':
 test()

• Paulo Sergio Schlogl: thank you for your suggestions! I will check it out and sorry because English it is not my native language...some times I do a lot of mistakes. But for me it is much important to check if the statistics function chromosome_gc_variation is doing what it is need to do. But, Is always excellent to learn to write a good/efficient and beautiful code. Thank you for your time and tips. They are awesome. Paulo