3
\$\begingroup\$

I wrote this code to count the dinucletide fractions that appear in a genome (this is a sequence of two nucleotides, 'A','G','C','T', together). My code also calculates the j2 index (this is a very simple index showing the fraction of groups of dinucleotides based on an equation at a glance*).

It has currently taken 3 days to generate a dataframe of 14gb, and there is still a lot to go, I'm wondering if I can improve the performance/speed with the code at all)

This is my code:

from Bio import SeqIO
import sys
from collections import Counter

def chunks(l, n):
  for i in range(0, len(l)-(n-1)):
    yield l[i:i+n]

def species_name_function(infile):
    for record in SeqIO.parse(infile, "fasta"):
        fields = record.description.split()
        speciesname = " ".join(fields[1:3])
    return speciesname


if __name__ == '__main__':
    frequency = []
    infile = sys.argv[1]
    for fasta in SeqIO.parse(open(infile), "fasta"):
        dna = str(fasta.seq)
        freq = Counter(dna)
        freq.update(Counter(chunks(dna,2)))
        frequency.append(freq)
        species_name = species_name_function(infile)

        genomesize = freq['A'] + freq['G'] + freq['C'] + freq['T']

        FYY = (freq['TT'] + freq['CC'] + freq['TC'] + freq['CT']) / genomesize
        FRR = (freq['AA'] + freq['GG'] + freq['AG'] + freq['GA']) / genomesize
        FYR = (freq['TA'] + freq['TG'] + freq['CA'] + freq['CG']) / genomesize
        FRY = (freq['AT'] + freq['AC'] + freq['GT'] + freq['GC']) / genomesize

        J2 = FYY + FRR - FYR - FRY

        listofbases = ["A", "C", "G", "T"]

        for base in listofbases:
            for base_2 in listofbases:
                towrite = base + base_2 + '\t' + str(freq[base + base_2]/genomesize) + '\t' + species_name + '\t' + str(genomesize) + '\t' + str(J2) + '\t' + infile + '\n'
                with open("resultsdinuc.csv", "a") as myfile:
                    myfile.write(towrite)

Example input (the files are actually a lot bigger than this in reality):

>NZ_NEDJ01000100.1 Halorubrum ezzemoulense DSM 17463 NODE_100_length_8476_cov_12.335, whole genome shotgun sequence
ACCGACACCATATGAGCGACGCGCCGACGACTGCGCCCTGCGACGCCTGCGGCGAGGCCACGACGGACGCGCTCGCGCGC
ACCGTCCGGCTGAGCGTCGACCGGGCGAACATCGACACCCAGCGGCTCTGCCCCGACTGCTTCGCCGACTGGATCCAGCG
CTACCAGGACCGCCTCGGCTCCGGCGACGACGGGGGCGACGAGAGCTCCGAGATCATCGTCGACTGAGGCCGAACGCGTT
CGCGTCGGCCGGCAACGTCCGTCTCGACCGCCCGTCTTAAGCCCCGGCGGGACGGACGCCGTGGTAATGGATC
>NZ_NEDJ01000108.1 Halorubrum ezzemoulense DSM 17463 NODE_108_length_6789_cov_9.46893, whole genome shotgun sequence
TGGCGTCGAGCGGCTCGGCCCGAAATTCTATTACCCCAAGTTCCGCAAGTTCTGATAGCCTCTGGCCGAAGGCAGGACGG
TCTTCATACATACCCGTTTTTGCCGGGCCAGAGGCACTAATGCTCCTGGTTCCGCCAGTCTACTGAAGAGCGTCGTCGCT
TAACGGTCGATTCGTTCCGCTCAGCGAGCCCCCGAACGAGGTAAGAGAACGCTGTAAAGGATTTATACTGCGAGGACGAG
GCCCGAGTGTGGTCGGACTCGCACGCGGGACCGTCGAAGTCGTGCCGTATCAGGAGTCGTGGAGCGACGCGTACGACGGG
GAGGTGGCTCGGTTACGGAGCGCAGTCGGTGATCGCGTCCGTCAGTTCGAACACATCGGCAGTACCGCGGTCGAGGGGAT
GGCGGCCAAGCCGATACTCGACGTGCTCGCCGTAGTCGACGAATCGACGACCGCGAGCGACCTCGTCCCAGCGCTCGAAA
CGCACGGCTACGAACGGCGCCCCGATGAGGTGGACGGGCGGGTGTTCCTCGCGAAGGGACCGCCAGAGAATCGTACGTGC
TATCTGTCGATCGCCGAAGTCGGAAGC
>NZ_NEDJ01000109.1 Halorubrum ezzemoulense DSM 17463 NODE_109_length_6759_cov_12.5481, whole genome shotgun sequence
GGCCCGATCCCGCCCGCGAGCTGCGCCGGGACCGCCACGAACCCGTCGCCGGGAGCGAGCGTCGGCTGCATGCTCCCGGT
CTCGACGTAGCTGAGGAGGACCGGTTGGCCGAGGAGCTGTCCGACGACCAGCGAGACGACGACCAGCACCGCGGCCGCTT
GGAGCGCGACGGACAGCGTTCGTTTGAGTGACATGGTGTCGAACTCGGCTCGGAGACGGACTCGGGGCGGCGACCGCCGC
GAGGCGGTACCTGTCGCGCGGCCGTCAGGTAGTCGTCGATCGCTGAACGGCGGCGTGTCCTTATAACTTCGTGGGTGGCG
GCGAACCGGATCGGGCGGCCGCCGTCGGCCCTACTCGTCGAAGGCGCCGGCGGCGAGCAGCGCGAACGGGCCGATGAACC
CGAGGCAGAACCCGAGCAGGTGGACGTACAGGTTCACGACGCTCCCGCCGCCCGAGGGGTCGGACGGGAAGCCGACGACG
GGGTAGCCGACGGCCACCAGCCCGCCGACGACCAGCAGTTCCCCGTAGCCGGGACGCGACGCGACGGCGGCCGCGTGCGT
CCGGACCGCGGGCGGGAACCGGACCTCGCTGCTGGCCGCGTACAGCACCGCGAGGAGCGCGCCCGCGACGCCGATCGCCA
GCCCCGCGACCCCGATGCCCGTCGTCGAGACCGGCAGCGCGAGGAGCGCGATCCACCCGACGAGCGCGAAGAAGACCGCC
GGCAGCGCACGCAGCGAGGCCGCGGGCGCGAAGTGTTCGCGGGCGTAGCAGTACCACAGGATCGGGAGCAGTCCGGCGAG
CGCCATGTTGATCCCGGAGAAGCCGAAGCCGATCGCGTCCCGCGGGACGGCGAGGTTCAGCGCGGACAGCGCGAACGGGA
AGGCGCCGAGGTACGTCGCGAGCGACGTGAAGAAGAGCCGCCGTCGCCCGCCAAGGACGGCGAGCGCGTAGC

Example output:

AA  0.0141576215    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
AC  0.0624115149    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
AG  0.0366918358    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
AT  0.0165172251    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
CA  0.0284332232    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
CC  0.0970976876    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
CG  0.1910099103    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
CT  0.0486078339    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
GA  0.0777489382    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
GC  0.1178621992    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
GG  0.092732421 Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
GT  0.0658329401    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
TA  0.0093204342    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
TC  0.0878952336    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
TG  0.0337423313    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
TT  0.0198206701    Halorubrumezzemoulense  8476    -0.0503775366   ./GCF_002114285.1_ASM211428v1_genomic.fna
AA  0.0378834927    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
AC  0.0679418706    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
AG  0.0491864365    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
AT  0.0475717302    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
CA  0.0544031797    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
CC  0.0710470749    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
CG  0.1038380325    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
CT  0.06322196  Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
GA  0.0739038629    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
GC  0.0694323686    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
GG  0.0614830456    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
GT  0.0715439076    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
TA  0.0363929947    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
TC  0.0840889331    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
TG  0.0617314619    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
TT  0.0462054403    Halorubrumezzemoulense  8051    -0.0258353  ./GCF_002114285.1_ASM211428v1_genomic.fna
AA  0.018964836 Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
AC  0.0595285131    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
AG  0.0431976821    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
AT  0.0243645463    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
CA  0.0296325563    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
CC  0.0949558804    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
CG  0.1818780456    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
CT  0.0416172791    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
GA  0.0817858554    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
GC  0.1065455024    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
GG  0.0902146714    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
GT  0.0694060319    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
TA  0.0156723298    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
TC  0.0871855657    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
TG  0.0325299618    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna
TT  0.0223890425    Halorubrumezzemoulense  7593    -0.0392466746   ./GCF_002114285.1_ASM211428v1_genomic.fna

How can I best improve this code? (I'm using it for hundreds of thousands of files, and when I use this, I'm using it in a bash loop:

for FAA in $(find . -name "*.fna")
do
    python3 dinucleotidescript.py $FAA
done

*Nucleotides can be divided into purines (R) which are either A or G, and pyramidines (Y) which are either T or C, the equation for J2 index is:

J2 index = FYY + FRR - FYR - FRY

F is fraction, so FYY for example is a the fraction of the genome where a T or C is followed by another T or C.

Improve this question
\$\endgroup\$

1 Answer 1

Reset to default
4
\$\begingroup\$

Welcome back Biomage :)

I can see you have improved your Python. However there is still room left for improvement!

  • Don't put everything under the if __name__ == '__main__' part

    Doing so will make it less readable and other programs can't import anything from this script. Instead make different functions. And have only the sys.argv lines in your main.

  • Split your code into functions

    Splitting your code into multiple functions can make it easier to your test code, and now each part can be imported by another script.

  • Avoid overhead by opening a file only once

    You are opening the file when you want to calculate the frequencies, and when you want the species name. This is unnecessary since it can be done in one read.

  • What's the point of this frequency = []

    This variable is never used later on, I suggest removing it.

  • In your chunks function for i in range(0, len(l)-(n-1)) the (n-1) part is redundent

    It would yield the same chunks as for i in range(0, len(l)-n-1) or even for i in range(0, len(l)-1).

  • Your Counter can be simplified

    1. There is no need to parse the fasta.seq to a string, Counter happily reads the sequence without parsing.

    2. freq.update(Counter(chunks(dna,2))) Here the intermediate Counter is redundant. It will update anyway.

  • I think this listofbases = ["A", "C", "G", "T"] should be a global variable since it can be used multiple places.

  • Your for a in listofbases: for b in listofbases: can be done with a builtin itertools.product

  • Don't concat strings, but use the better "str".format() or f"strings"

  • At the write_output part, you keep opening and closing the file, to avoid overhead only do this once.

  • You do / genome_size 4 times, but instead you could add them all and only then divide by genome_size

Revised Code

from Bio import SeqIO
import sys
from collections import Counter
from itertools import product

LISTOFBASES = ["A", "C", "G", "T"]
RESULT_FILE = "resultsdinuc.csv"

def chunks(l, n):
    for i in range(0, len(l)-1):
        yield l[i:i+n]

def parse_file(infile):
    for fasta in SeqIO.parse(open(infile), "fasta"):
        freq = Counter(fasta.seq)
        freq.update(chunks(str(fasta.seq), 2))
        species_name = "".join(fasta.description.split()[1:3])
        genome_size = sum(freq[i] for i in LISTOFBASES)

        FYY = (freq['TT'] + freq['CC'] + freq['TC'] + freq['CT'])
        FRR = (freq['AA'] + freq['GG'] + freq['AG'] + freq['GA'])
        FYR = (freq['TA'] + freq['TG'] + freq['CA'] + freq['CG'])
        FRY = (freq['AT'] + freq['AC'] + freq['GT'] + freq['GC'])
        J2 = (FYY + FRR - FYR - FRY) / genome_size

        yield J2, species_name, freq, genome_size

def write_output(J2, species_name, frequency, genome_size, infile):
    with open(RESULT_FILE, "a") as myfile:
        for prod in product(LISTOFBASES, LISTOFBASES):
            base = ''.join(prod)
            result = f"{base}\t{frequency[base]/genome_size}\t{species_name}\t{genome_size}\t{J2}\t{infile}\n"
            # For python 3.5 and lower use .format()
            # result = "{}\t{}\t{}\t{}\t{}\t{}".format(base, frequency[base]/genome_size, species_name, genome_size, J2, infile)
            myfile.write(result)

if __name__ == '__main__':
    infile = sys.argv[1]
    for args in parse_file(infile):
        write_output(*args, infile)
Improve this answer
\$\endgroup\$
8
  • \$\begingroup\$ Nitpicking: len(l)-(n-1) is not the same as len(l)-n-1. \$\endgroup\$
    – vnp
    Aug 20, 2018 at 23:50
  • \$\begingroup\$ @vnp You are right, What I meant to say was, it will yield the same chunks. I'll clarify a bit \$\endgroup\$
    – Ludisposed
    Aug 21, 2018 at 7:08
  • \$\begingroup\$ Thanks! I'm learning! Although getting error: File "dinucleotidev2.py", line 32 result = f"{base}\t{frequency[base]/genome_size}\t{species_name}\t{genome_size}\t{J2}\t{infile}\n" ^ SyntaxError: invalid syntax File "dinucleotidev2.py", line 32 result = f"{base}\t{frequency[base]/genome_size}\t{species_name}\t{genome_size}\t{J2}\t{infile}\n" ^ SyntaxError: invalid syntax \$\endgroup\$
    – Biomage
    Aug 21, 2018 at 13:19
  • \$\begingroup\$ @Biomage What version of Phyton are you using? From 3.6 they added f"strings", you could use format instead \$\endgroup\$
    – Ludisposed
    Aug 21, 2018 at 13:31
  • 1
    \$\begingroup\$ @Biomage It's never a bad idea to think about upgrading your python version, lot's of goodies keep being added ;) But for this you can rewrite it to format \$\endgroup\$
    – Ludisposed
    Aug 21, 2018 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.