Managing and substringing multiple strings in a matrix (Python & Numpy)

user_10514661 0
Viewed 117

I have a file titled "test.fa" that reads:

>Sequence 1
TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT
>Sequence 2
CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG
>Sequence 3
TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG
>Sequence 4
AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC

This test file only has 8 lines, but could have more. I also don't now the length of all sequences. So first I count the number of lines, and create a matrix based on the number of lines.

import numpy as np
filename = "test.fa"
with open(filename, 'r') as file:
    n1 = 0
    n2 = 0
    for line in file.readlines():
        if line[0] != '>':
            m1 = 1
            n1 = n1 + m1
        else:
            m2 = 1
            n2 = n2 + m2
    seq = np.chararray((n1,2),itemsize = 99)

Next, I add values to the matrix.

with open(filename, 'r') as file:
    n3 = -1
    n4 = -1
    for line in file.readlines():
        if line[0] != '>':
            m3 = 1
            n3 = n3 + m3
            seq[n3,1] = line
        else:
            m4 = 1
            n4 = n4 + m4
            seq[n4,0] = line

If I call 'seq' I get:

chararray([[b'>Sequence 1', b'TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT'],
       [b'>Sequence 2', b'CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG'],
       [b'>Sequence 3', b'TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG'],
       [b'>Sequence 4', b'AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC']], dtype='|S99')

seq[2,1]: b'TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG'
seq[2,1][0] : 84
seq[2,1][0:8] : b'TCGACCCT'
seq[2,8][8] : 67
seq[2,1][8:9] : C

This doesn’t feel like the best way manage this data (and why does [2,1][0] return 84?). I am not sure np.chararray with itemsize of 99 is the best approach.

My question is: What is a better way to organize/manage this data? I eventually need to count the occurrences of each nucleotide ( i.e. how many As, Ts, Cs, Gs), and pull substrings from each sequence. For context, this relates to Expectation Maximization and Gibbs Sampling.

substring numpy python
3 Answers

it returns 84 while the type is S99 because the maximum length found was 99, so the type would be the S99, how you read from the file is not that right, and a better way to manage ur data is using dict, I dont know if it will suit your next needs but here is it:

seq={}
with open("file") as f:
    for line in f:
        seq[line]=next(f)

here it will read two lines at a time, so the first line has the sequence 1, next(f) has the seq TGC....., if happened that the file is not consistently ordered use if statement

user_84116090
user_86724570

By using np.chararray, you get an array with an S dtype, bytestrings:

In [145]: x = np.chararray((1,),10)                                                            
In [146]: x[0]='ACT'                                                                           
In [147]: x                                                                                    
Out[147]: chararray([b'ACT'], dtype='|S10')
In [148]: x[0]                                                                                 
Out[148]: b'ACT'
In [149]: x[0][0]                                                                              
Out[149]: 65
In [150]: x[0][1]                                                                              
Out[150]: 67
In [151]: x[0][2]                                                                              
Out[151]: 84

That's one byte per character. ASCII code of A is 65; T is 19 characters further.

Did you see the chararray note:

The chararray class exists for backwards compatibility with Numarray, it is not recommended for new development. Starting from numpy 1.4, if one needs arrays of strings, it is recommended to use arrays of dtype object_, string_ or unicode_, and use the free functions in the numpy.char module for fast vectorized string operations.

In [159]: x = np.array([b'ACT'])                                                               
In [160]: x                                                                                    
Out[160]: array([b'ACT'], dtype='|S3')
In [161]: x[0]                                                                                 
Out[161]: b'ACT'
In [162]: x[0][0]                                                                              
Out[162]: 65

If instead we use the Py3 default string type (unicode):

In [154]: x = np.array(['ACT'])                                                                
In [155]: x                                                                                    
Out[155]: array(['ACT'], dtype='<U3')
In [156]: x[0]                                                                                 
Out[156]: 'ACT'
In [157]: x[0][0]                                                                              
Out[157]: 'A'
In [158]: x[0][1]                                                                              
Out[158]: 'C'

===

I think collecting values in a list or dictionary would be simpler:

In [165]: txt = """>Sequence 1 
     ...: TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT 
     ...: >Sequence 2 
     ...: CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG 
     ...: >Sequence 3 
     ...: TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG 
     ...: >Sequence 4 
     ...: AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC"""                                         
In [166]: alist = []                                                                           
In [167]: for row in txt.splitlines(): 
     ...:     if row.startswith('>'): 
     ...:         x = [row] 
     ...:     else: 
     ...:         x.append(row) 
     ...:         alist.append(x) 
     ...:                                                                                      
In [168]: alist                                                                                
Out[168]: 
[['>Sequence 1', 'TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT'],
 ['>Sequence 2', 'CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG'],
 ['>Sequence 3', 'TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG'],
 ['>Sequence 4', 'AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC']]
In [169]: dd = {key:val for key,val in alist}                                                  
In [170]: dd                                                                                   
Out[170]: 
{'>Sequence 1': 'TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT',
 '>Sequence 2': 'CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG',
 '>Sequence 3': 'TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG',
 '>Sequence 4': 'AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC'}

So basic processing of those strings:

In [171]: [len(val) for key,val in alist]                                                      
Out[171]: [41, 41, 44, 41]
In [172]: [val.count('T') for key,val in alist]                                                
Out[172]: [16, 7, 6, 6]
In [173]: [val.count('C') for key,val in alist]                                                
Out[173]: [12, 21, 15, 10]
user_9019250
Related