Occurrence Frequency of the Python DNA sequence neutron Sequence

There are many things we need to pay attention to when using Python DNA sequences. In fact, there are many problems in constant learning. Next we will look at in detail how to conduct related technical schools. Ms is my brother's rotation project: Given A bunch of Python DNA sequences, that is, A string consisting of characters A, C, G, T, count the occurrence frequency of all subsequences whose lengths are n.

For example, ACGTACGT, the sub-sequence length is 2, so AC = 2, CG = 2, GT = 2, TA = 1, and the sub-sequence frequency of the remaining length is 0.

The first thing that comes to mind is to create a dictionary. The key is all possible subsequences, and the value is the frequency at which this subsequence appears. Python DNA sequences, but when the sub-sequences are relatively long, such as n = 8, A 65536 (4 to the power of 8) key-value pair dictionary is required, the length of each key is 8 characters. In this case, ms is a waste of memory ..

As a result, all subsequences whose lengths are n are ordered and continuous, so they can be mapped to a list with 4 n ing power. Make A = 0, C = 1, G = 2, T = 3, then, the sub-sequence ACGT is converted to 0*4 ^ 3 + 1*4 ^ 2 + 2*4 + 3 = 27 and mapped to the 27th-bit list. In this way, the index of the list corresponds to the subsequence, and the index position of the list stores the occurrence frequency of the subsequence.

So we need to create two dictionaries to indicate the one-to-one correspondence between ACGT and 0123:

 
 

  
  
i2mD = {0:'A', 1:'C', 2:'G', 3:'T'}  
  
  
m2iD = dict(A=0,C=1,G=2,T=3)  
  
  
# This is just another way to initialize a 
dictionary 
 
 

And the following subsequences are mapped to integer functions:

 
 

  
  
def motif2int(motif):  
  
  
'''convert a sub-sequence/motif to a non-negative 
integer'''  
  
  
total = 0 
  
  
for i, letter in enumerate(motif):  
  
  
total += m2iD[letter]*4**(len(motif)-i-1)  
  
  
return total  
  
  
Test:  
  
  
>>> motif2int('ACGT') 
 
 

Although we use subsequences as positive integers to store them internally, actually, this integer is not stored in memory, but represented by its index in the list.) to make it easier for scientists to see, it is better to convert the data back to the subsequence during output.

Therefore, the following integer is mapped to a sub-Sequence function, and another function baseN () is called. The source is here. Thanks to the author ~

 
 

  
  
def baseN(n,b):  
  
  
'''convert non-negative decimal integer n to  
  
  
equivalent in another base b (2-36)'''  
  
  
return ((n == 0) and '0' ) or ( baseN(n // b, b).lstrip('0') + \  
  
  
"0123456789abcdefghijklmnopqrstuvwxyz"[n % b])  
  
  
def int2motif(n, motifLen):  
  
  
'''convert non-negative integer n to a sub-sequence/motif with length motifLen'''  
  
  
intBase4 = baseN(n,4)  
  
  
return ''.join(map(lambda x: i2mD[int(x)],'0'*(motifLen-len(intBase4))+intBase4))  
  
  
Test:  
  
  
>>> int2motif(27,4)  
  
  
'ACGT'  
 
 

The following code reads a fasta file containing a DNA sequence from the command line and the length of the sub-sequence, and outputs the sub-sequence and frequency. Note that the following code requires the Biopython module.

 
 

  
  
if __name__ == '__main__':  
  
  
import sys  
  
  
from Bio import SeqIO  
  
  
# read in the fasta file name and motif length  
  
  
# from command line parameters  
  
  
fastafile = sys.argv[1]  
  
  
motifLen = int(sys.argv[2])  
  
  
# list to store subsequence frequency  
  
  
frequencyL = [0]*4**motifLen  
  
  
# go over each DNA sequence in the fasta file  
  
  
# and count the frequency of subsequences  
  
  
it = SeqIO.parse(open(fastafile),'fasta')  
  
  
for rec in it:  
  
  
chrom = rec.seq.tostring()  
  
  
for i in range(len(chrom)-motifLen+1):  
  
  
motif = chrom[i:i+motifLen]  
  
  
frequencyL[motif2int(motif)] += 1  
  
  
# print frequency result to screen  
  
  
for i, frequency in enumerate(frequencyL):  
  
  
print int2motif(i, motifLen), frequency  
 
 

The above describes the Python DNA sequence.