Strand Direction Detection

The VirjenDB defines a field that holds information on sequence direction:

0ie field summary

We implemented a script to determine strand direction information from sequence data. The algorithm is taken from a ViralClust snippet and does the following:

flow chart1 <!– @startuml start

:fetch list of unprocessed entries from Elasticsearch → ids; :create CSV file (outfile) to collect calculated data; :write head (virjen_id , is_coding_strand) → outfile;

while (for each id in ids) is (id) :fetch sequence for id → sequence; :indentify direction of sequence → is_coding_strand; :attach id and is_coding_strand → outfile; endwhile

:close outfile; stop @enduml

–>

The identification of sequence direction works as follows:

seqeunce direction determination

The reverse complement of a strand is created as follows:

def reverseComplement(sequence):
    comp = {
        'A' : 'T',
        'C' : 'G',
        'G' : 'C',
        'T' : 'A'}
    return(''.join([comp.get(x,'N')  for x in sequence.upper()[::-1]]))

the longest coding region of a strand is calculated as follows:

longest coding region calculation

Here is the acutal python code:

CODON2PEPTIDE = {
    'ATA':'I', 'ATC':'I', 'ATT':'I', 'ATG':'M',
    'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACT':'T',
    'AAC':'N', 'AAT':'N', 'AAA':'K', 'AAG':'K',
    'AGC':'S', 'AGT':'S', 'AGA':'R', 'AGG':'R',
    'CTA':'L', 'CTC':'L', 'CTG':'L', 'CTT':'L',
    'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCT':'P',
    'CAC':'H', 'CAT':'H', 'CAA':'Q', 'CAG':'Q',
    'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGT':'R',
    'GTA':'V', 'GTC':'V', 'GTG':'V', 'GTT':'V',
    'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCT':'A',
    'GAC':'D', 'GAT':'D', 'GAA':'E', 'GAG':'E',
    'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGT':'G',
    'TCA':'S', 'TCC':'S', 'TCG':'S', 'TCT':'S',
    'TTC':'F', 'TTT':'F', 'TTA':'L', 'TTG':'L',
    'TAC':'Y', 'TAT':'Y', 'TAA':'*', 'TAG':'*',
    'TGC':'C', 'TGT':'C', 'TGA':'*', 'TGG':'W',
}

# Asterisks mark stop codons

REGEX_ORF = re.compile(r'[^*]{200,}')

def find_longest_cds(sequence):
    combined_cds_len = 0
    for frame in range(3):
        proteinSequence = ""
        for fragment in range(frame, len(sequence), 3):
            codon = sequence[fragment:fragment+3]
            if len(codon) != 3:
                continue
            try:
                proteinSequence += CODON2PEPTIDE[codon]
            except KeyError:
                proteinSequence += 'X'
        matches = REGEX_ORF.findall(proteinSequence)
        allORFs = "".join([x for x in matches if x])
        combined_cds_len += len(allORFs)
    return combined_cds_len

The whole code can be found here: github.com/VirJenDB/Curation-module

Documentation

Table of Contents