mgkit.utils.sequence module

Module containing functions related to sequence data

Note

For those functions without a docstring, look at the same with a underscore (‘_’) prepended.

class mgkit.utils.sequence.Alignment(seqs=None)

Bases: object

Simple alignment class

add_seq(name, seq)

Add a sequence to the alignment

Parameters:
  • name (str) – name of the sequence
  • seq (str) – sequence
add_seqs(seqs)

Add sequences to the alignment

Parameters:seqs (iterable) – iterable that returns (name, seq)
get_consensus(nucl=True)

Changed in version 0.1.16: added nucl parameter

The consensus sequence is constructed by checking the nucleotide that has the maximum number of counts for each position in the alignment.

Parameters:nucl (bool) – specify if the alignment is nucleotidic
Returns:consensus sequence
Return type:str
get_position(pos)

Get all characters at a position

Parameters:pos (int) – position to return (0-based)
Return str:all characters occuring at the position
get_seq_len()

Get the length of the alignment

get_snps(ref_seq=None, full_size=False)

A SNP is called for the nucleotide that has the most counts among the ones that differ in the each site of the alignment. If two nucleotides have the same maximum count, one is randomly chosen.

Parameters:
  • ref_seq (str) – a reference sequence can be provided, if None, a consensus sequence is produced for the alignment
  • full_size (bool) – if True a tuple is returned for each position in the alignment. If there is no SNP at a position the value for the SNP is None
Return list:

a list of tuples (position, SNP)

mgkit.utils.sequence._get_kmers(seq, k)

New in version 0.2.6.

Returns a generator, with every iteration yielding a kmer of size k

Parameters:
  • seq (str) – sequence
  • k (int) – kmer size
Yields:

str – a portion of seq, of size k with a step of 1

mgkit.utils.sequence._sequence_signature(seq, w_size, k_size=4, step=None)

New in version 0.2.6.

Returns the signature of a sequence, based on a kmer length, over a sliding window. Each sliding window signature is placed in order into a list, with each element being a collections.Counter instance whose keys are the kmer found in that window.

Parameters:
  • seq (str) – sequence for which to get the signature
  • w_size (int) – size of the sliding window size
  • k_size (int) – size of the kmer to use get_kmers()
  • step (int) – step to use in sliding_window()
Returns:

a list of collections.Counter instances, for each window used

Return type:

list

mgkit.utils.sequence._signatures_matrix(seqs, w_size, k_size=4, step=None)

New in version 0.2.6.

Return a matrix (pandas.DataFrame) where the columns are the kmer found in all sequences seqs and the rows are the a MultiIndex with the first level being the sequnce name and the second the index of the sliding window for which a signature was computed.

Parameters:
  • seqs (iterable) – iterable that yields a tuple, with the first element being the sequence name and the second the sequence itself
  • w_size (int) – size of the sliding window size
  • k_size (int) – size of the kmer to use get_kmers()
  • step (int) – step to use in sliding_window(), defaults to half of the window size
Returns:

a DataFrame where the columns are the kmers and the rows are the signatures of each contigs/windows.

Return type:

pandas.DataFrame

mgkit.utils.sequence._sliding_window(seq, size, step=None)

New in version 0.2.6.

Returns a generator, with every iteration yielding a subsequence of size size, with a step of step.

Parameters:
  • seq (str) – sequnece
  • size (int) – size of the sliding window
  • step (int, None) – the step to use in the sliding window. If None, half of the sequence length is used
Yields:

str – a subsequence of size size and step step

mgkit.utils.sequence.calc_n50(seq_lengths)

Calculate the N50 statistics for a numpy.array of sequence lengths.

The algorithm finds in the supplied array the element (contig length) for which the sum all contig lengths equal or greater than it is equal to half of all assembled base pairs.

Parameters:seq_lengths (array) – an instance of a numpy array containing the sequence lengths
Return int:the N50 statistics value
mgkit.utils.sequence.check_snp_in_seq(ref_seq, pos, change, start=0, trans_table=None)

Check a SNP in a reference sequence if it is a synonymous or non-synonymous change.

Parameters:
  • ref_seq (str) – reference sequence
  • pos (int) – SNP position - it is expected to be a 1 based index
  • change (str) – nucleotide change occuring at pos
  • start (int) – the starting position for the coding region - 0 based index
  • trans_table (dict) – translation table used - codon->AA
Return bool:

True if it is a synonymous change, False if non-synonymous

mgkit.utils.sequence.convert_aa_to_nuc_coord(start, end, frame=0)

Converts aa coordinates to nucleotidic ones. The coordinates must be from ‘+’ strand. For the ‘-‘ strand, use reverse_aa_coord() first.

Parameters:
  • start (int) – start of the annotation (lowest number)
  • end (int) – end of the annotation (highest number)
  • frame (int) – frame of the AA translation (0, 1 or 2)
Returns:

the first element is the converted start and the second element is the converted end

Return type:

tuple

Note

the coordinates are assumed to be 1-based indices

mgkit.utils.sequence.extrapolate_model(quals, frac=0.5, scale_adj=0.5)

New in version 0.3.3.

Extrapolate a quality model from a list of qualities. It uses internally a LOWESS as the base, which is used to estimate the noise as a normal distribution.

Parameters:
  • quals (list) – list of arrays of qualities, sorted by position in the corresponding sequence
  • frac (float) – fraction of the data used for the LOWESS fit (uses statsmodels)
  • scale_adj (float) – value by which the scale of the normal distribution will be multiplied. Defaults to halving the scale
Returns:

the first element is the qualities fit with a LOWESS, the second element is the distribution

Return type:

tuple

mgkit.utils.sequence.get_contigs_info(file_name, pp=False)

Changed in version 0.2.4: file_name can be a dict name->seq or a list of sequences

New in version 0.2.1.

Given a file name for a fasta file with sequences, a dictionary of name->seq, or a list of sequences, returns the following information in a tuple, or a string if pp is True:

  • number of sequences
  • total base pairs
  • max length
  • min length
  • average length
  • N50 statistic
Parameters:
  • file_name (str) – fasta file to open
  • pp (bool) – if True, a formatted string is returned
Returns:

the returned value depends on the value of pp, if True a formatted string is returned, otherwise the tuple with all values is.

Return type:

str, tuple

mgkit.utils.sequence.get_seq_expected_syn_count(seq, start=0, syn_matrix=None)

Calculate the expected number of synonymous and non-synonymous changes in a nucleotide sequence. Assumes that the sequence is already in the correct frame and its length is a multiple of 3.

Parameters:
  • seq (iterable) – nucleotide sequence (uppercase chars)
  • start (int) – frame of the sequence
  • syn_matrix (dict) – dictionary that contains the expected number of changes for a codon, as returned by get_syn_matrix()
Return tuple:

tuple with counts of expected counts (syn, nonsyn)

mgkit.utils.sequence.get_seq_number_of_syn(ref_seq, snps, start=0, trans_table=None)

Given a reference sequence and a list of SNPs, calculates the number of synonymous and non-synonymous SNP.

Parameters:
  • ref_seq (str) – reference sequence
  • snps (iterable) – list of tuples (position, SNP) - zero based index
  • start (int) – the frame used for the reference {0, 1, 2}
  • trans_table (dict) – translation table used - codon->AA
Return tuple:

synonymous and non-synonymous counts

mgkit.utils.sequence.get_syn_matrix(trans_table=None, nuc_list=None)

Returns a dictionary containing the expected count of synonymous and non-synonymous changes that a codon can have if one base is allowed to change at a time.

There are 9 possible changes per codon.

Parameters:
  • trans_table (dict) – a tranlation table, defaults to seq_utils.TRANS_TABLE
  • nuc_list (iterable) – a list of nucleotides in which a base can change, default to the keys of seq_utils.REV_COMP
Return dict:

returns a dictionary in which for each codon a dictionary {‘syn’: 0, ‘nonsyn’: 0} holds the number of expected changes

mgkit.utils.sequence.get_syn_matrix_all(trans_table=None)

Same as get_syn_matrix() but a codon can change in any of the ones included in trans_table.

There are 63 possible changes per codon.

mgkit.utils.sequence.get_variant_sequence(seq, *snps)

New in version 0.1.16.

Return a sequence changed in the positions requested.

Parameters:
  • seq (str) – a sequence
  • *snps (tuple) – each argument passed is a tuple with the first element as a position in the sequence (1-based index) and the second element is the character to substitute in the sequence
Returns:

string with the changed characters

Return type:

str

Example

>>> get_variant_sequence('ACTGATATATGCGCGCATCT', (1, 'C'))
'CCTGNTGTATGCGCGCATCT'

Note

It is used for nucleotide sequences, but it is valid to use any string

mgkit.utils.sequence.make_reverse_table(tbl=None)

Makes table to reverse complement a sequence by reverse_complement(). The table used is the complement for each nucleotide, defaulting to REV_COMP

mgkit.utils.sequence.put_gaps_in_nuc_seq(nuc_seq, aa_seq, trim=True)

Match the gaps in an amino-acid aligned sequence to its original nucleotide sequence. If the nucleotide sequence is not a multiple of 3, the trim option by default trim those bases from the output.

Parameters:
  • nuc_seq (str) – original nucleotide sequence
  • aa_seq (str) – aligned amino-acid sequence
  • trim (bool) – if True trim last nucleotide(s)
Return str:

gapped nucleotide sequence

mgkit.utils.sequence.qualities_model_constant(length=150, scale=1, loc=35)

New in version 0.3.3.

Model with constant quality

Parameters:
  • length (int) – length of the qualities
  • scale (float) – base level of the qualities
  • loc (float) – loc parameter of the normal distribution
Returns:

first element is sequence qualities, the second element contains the distribution used to randomise them

Return type:

tuple

mgkit.utils.sequence.qualities_model_decrease(length=150, scale=None, loc=35)

New in version 0.3.3.

The model is a decreasing one, from 35 and depends on the length of the sequence.

Parameters:
  • length (int) – length of the qualities
  • scale (float) – base level of the qualities
  • loc (float) – loc parameter of the normal distribution
Returns:

first element is sequence qualities, the second element contains the distribution used to randomise them

Return type:

tuple

mgkit.utils.sequence.random_qualities(n=1, length=150, model=None)

New in version 0.3.3.

Parameters:
  • n (int) – number of quality arrays to yield
  • length (int) – length of the quality array
  • model (tuple) – a tuple specifying the qualities and error distribution, if None qualities_model_decrease() is used
Yields:

numpy.array – numpy array of qualities, with the maximum value of 40

mgkit.utils.sequence.random_sequences(n=1, length=150, p=None)

New in version 0.3.3.

Returns an iterator of random squences, where each nucleotide probability can be customised in the order (A, C, T, G)

Parameters:
  • n (int) – number of sequences to yield
  • length (int) – length of each sequence
  • p (tuple) – tuple with the probability of a nucleotide to occur, in the order A, C, T, G
Yields:

str – string representing a nucleotidic sequence

mgkit.utils.sequence.random_sequences_codon(n=1, length=150, codons=['CTT', 'TAG', 'ACA', 'AAA', 'ATC', 'AAC', 'ATA', 'AGG', 'CCT', 'ACT', 'AGC', 'AAG', 'AGA', 'CAT', 'AAT', 'ATT', 'CTG', 'CTA', 'CTC', 'CAC', 'TGG', 'CAA', 'AGT', 'CCA', 'CCG', 'CCC', 'TAT', 'GGT', 'TGT', 'CGA', 'CAG', 'TCT', 'GAT', 'CGG', 'TTT', 'TGC', 'GGG', 'TGA', 'GGA', 'TAA', 'ACG', 'TAC', 'TTC', 'TCG', 'TTA', 'TTG', 'TCC', 'ACC', 'TCA', 'GCA', 'GTA', 'GCC', 'GTC', 'GGC', 'GCG', 'GTG', 'GAG', 'GTT', 'GCT', 'GAC', 'CGT', 'GAA', 'ATG', 'CGC'], p=None, frame=None)

New in version 0.3.3.

Returns an iterator of nucleotidic sequences, based on a defined genetic code (passed as parameter, defaults to the universal one). The sequence if first sampled with replacement from the codon list, with a number of codons that covers the length chosen plus an additional one to allow a frame shift as set by frame

Note

If the probability (for each codon) are supplied, the number of sequences required to match those probabilities within a 10% margin of error is of at least 10.000 sequences, for 5% at leas 100.000

Parameters:
  • n (int) – number of sequences to yield
  • length (int) – length of the sequences
  • codons (iterable) – codons used when generating the sequences
  • p (tuple) – probability of each codon occurence, in the same order as codons
  • frame (int or None) – used to define a specific frame shift occuring in the sequence (0 to 2) or a random one (if None)
Yields:

str – string representing a nucleotidic sequence

mgkit.utils.sequence.reverse_aa_coord(start, end, seq_len)

Used to reverse amino-acid coordinates when parsing an AA annotation on the - strand. Used when the BLAST or HMMER annotations use AA sequences.

Parameters:
  • start (int) – start of the annotation
  • end (int) – end of the annotation
  • seq_len (int) – aa sequence length
Returns:

reversed (from strand - to strand +) coordinates. The first element is the converted start and the second element is the converted end

Return type:

tuple

Note

  • start and end are 1-based indices
mgkit.utils.sequence.reverse_complement(seq, tbl='\x00\x01\x02\x03\x04\x05\x06\x07\x08\t\n\x0b\x0c\r\x0e\x0f\x10\x11\x12\x13\x14\x15\x16\x17\x18\x19\x1a\x1b\x1c\x1d\x1e\x1f !"#$%&\'()*+, -./0123456789:;<=>?@TBGDEFCHIJKLMNOPQRSAUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~\x7f\x80\x81\x82\x83\x84\x85\x86\x87\x88\x89\x8a\x8b\x8c\x8d\x8e\x8f\x90\x91\x92\x93\x94\x95\x96\x97\x98\x99\x9a\x9b\x9c\x9d\x9e\x9f\xa0\xa1\xa2\xa3\xa4\xa5\xa6\xa7\xa8\xa9\xaa\xab\xac\xad\xae\xaf\xb0\xb1\xb2\xb3\xb4\xb5\xb6\xb7\xb8\xb9\xba\xbb\xbc\xbd\xbe\xbf\xc0\xc1\xc2\xc3\xc4\xc5\xc6\xc7\xc8\xc9\xca\xcb\xcc\xcd\xce\xcf\xd0\xd1\xd2\xd3\xd4\xd5\xd6\xd7\xd8\xd9\xda\xdb\xdc\xdd\xde\xdf\xe0\xe1\xe2\xe3\xe4\xe5\xe6\xe7\xe8\xe9\xea\xeb\xec\xed\xee\xef\xf0\xf1\xf2\xf3\xf4\xf5\xf6\xf7\xf8\xf9\xfa\xfb\xfc\xfd\xfe\xff')

Returns the reverse complement of a nucleotide sequence

Parameters:
Return str:

returns the reverse complement of a nucleotide sequence

mgkit.utils.sequence.reverse_complement_old(seq, tbl=None)

Returns the reverse complement of a nucleotide sequence

Parameters:
  • seq (str) – nucleotide sequence with uppercase characters
  • tbl (dict) – dictionary of complement bases, like REV_COMP
Return str:

returns the reverse complement of a nucleotide sequence

mgkit.utils.sequence.sequence_composition(sequence, chars=('A', 'C', 'T', 'G'))

New in version 0.1.13.

Returns the number of occurences of each unique character in the sequence

Parameters:
  • sequence (str) – sequence
  • chars (iterable, None) – iterable of the chars to test, default to (A, C, T, G). if None checks all unique characters in the sequencce
Yields:

tuple – the first element is the nucleotide and the second is the number of occurences in sequence

mgkit.utils.sequence.sequence_gc_content(sequence)

Changed in version 0.3.3: in case of ZeroDivisionError returns .5

New in version 0.1.13.

Calculate GC content information for an annotation. The formula is:

(1)\[\frac {(G + C)}{(G + C + A + T)}\]
Parameters:sequence (str) – sequence
Returns:GC content
Return type:float
mgkit.utils.sequence.sequence_gc_ratio(sequence)

New in version 0.1.13.

Calculate GC ratio information for a sequence. The formula is:

(2)\[\frac {(A + T)}{(G + C)}\]
Parameters:sequence (str) – sequence
Returns:GC ratio, or numpy.nan if G = C = 0
Return type:float
mgkit.utils.sequence.translate_sequence(sequence, start=0, tbl=None, reverse=False)

Translate a nucleotide sequence in an amino acid one.

Parameters:
  • sequence (str) – sequence to translate, it’s expected to be all caps
  • start (int) – 0-based index for the translation to start
  • tbl (dict) – dictionary with the translation for each codon
  • reverse (bool) – if True, reverse_complement() will be called and the returned sequence translated
Return str:

the translated sequence