.. _sequences: ========= Sequences ========= `EncodedArray` ============== Sequence data in BioNumPy are are represented by `EncodedArray` objects. These are basically `numpy` arrays of integers, that have an encoding that specifies which character each integer represent. This representataion allows us to do fast `numpy` operation on the sequences, while still allowing for human readable representation of them. The easiest way to create an `EncodedArray` is to use the `as_encoded_array` function. >>> import bionumpy as bnp >>> encoded_array = bnp.as_encoded_array("actggtcc") >>> encoded_array encoded_array('actggtcc') We see that the encoded array represents the text "actggtcc", but under the hood these are just integers with an ASCII encoding. You seldom have to think about this internal representation, but it is what allow us to write for instance: >>> print(encoded_array == "g") [False False False True True False False False] And get numpy-fast performance for the query. We can also use numpy-like indexing on encoded arrays, so that we can for instance trim the first and last two characters from the sequence: >>> encoded_array[2:-2] encoded_array('tggt') `EncodedRaggedArray` ==================== When working with multiple sequences we usually have to use `EncodedRaggedArray` objects. These are much like `EncodedArray` objects, but instead uses `npstructures.RaggedArray` to store the integers. This allows us to store seqeunces of differing lengths. The easiest way to create an `EncodedRaggedArray` is to use the `as_encoded_array` function on a list of strings: >>> encoded_ragged_array = bnp.as_encoded_array(["ctt", "actg", "ag"]) >>> encoded_ragged_array encoded_ragged_array(['ctt', 'actg', 'ag']) These objects also behave very much like numpy arrays, in indexing and broadcasting. For instance, to get the 2nd character of the first and third seqeunce: >>> encoded_ragged_array[[0, 2], 1] encoded_array('tg') Absolute basics of `EncodedRaggedArray` ======================================= As shown above, sequences are represented as `EncodedRaggedArray` objects. Below are some absolute basics of `EncodedRaggedArray` objects: For this example, instead of creating a `EncodedRaggedArray` object directly, we will simulate some sequences using the BioNumpy's `simulate_sequences` function: >>> import bionumpy as bnp >>> import numpy as np >>> from bionumpy.simulate import simulate_sequences >>> rng = np.random.default_rng(seed=1) >>> np.random.seed(1) >>> aa_alphabet = bnp.encodings.alphabet_encoding.AminoAcidEncoding.get_labels() >>> named_seqs = simulate_sequences(aa_alphabet, {f's{i}': np.random.randint(5,20) for i in range(10)}, rng) >>> my_seqs = named_seqs.sequence >>> named_seqs # print the sequences SequenceEntry with 10 entries name sequence s0 LMSYAEVYGH s1 WKGVGKQNCAWSVNVH s2 LTDHDL*DKKWFMGASC s3 GMMD*S*CSHNYG s4 SEH*KMHDKQLTIP s5 TYKASNWLICLQTVFP s6 TGIVPMRM*S s7 CENVC s8 RSTWF s9 NTIFMC Indexing and slicing of `EncodedRaggedArray` objects ----------------------------------------------------- We can index and slice `EncodedRaggedArray` objects in a similar way to numpy arrays. Below are some examples: >>> my_seqs[0:2] # first 2 sequences encoded_ragged_array(['LMSYAEVYGH', 'WKGVGKQNCAWSVNVH'], AlphabetEncoding('ACDEFGHIKLMNPQRSTVWY*')) >>> my_seqs[-4:] # last 4 sequences encoded_ragged_array(['TGIVPMRM*S', 'CENVC', 'RSTWF', 'NTIFMC'], AlphabetEncoding('ACDEFGHIKLMNPQRSTVWY*')) Some basic properties of `EncodedRaggedArray` objects ------------------------------------------------------ >>> my_seqs.shape # number of sequences and length of each sequence (10, array([10, 16, 17, 13, 14, 16, 10, 5, 5, 6])) >>> my_seqs.lengths # lengths of each sequence array([10, 16, 17, 13, 14, 16, 10, 5, 5, 6]) >>> my_seqs.size # total number of elements (amino acid residues across all sequences in the encoded ragged array) 112 >>> my_seqs.encoding # the encoding used for the sequences AlphabetEncoding('ACDEFGHIKLMNPQRSTVWY*') Concatenation of `EncodedRaggedArray` objects ------------------------------------------------ >>> np.concatenate([my_seqs, my_seqs[-2:]]).shape # concatenate two encoded ragged arrays and get the shape (12, array([10, 16, 17, 13, 14, 16, 10, 5, 5, 6, 5, 6])) Getting unique elements and counting occurrences ------------------------------------------------- >>> bnp.count_encoded(my_seqs.get_column_values(0)) # count the number of occurrences of each amino acid at the first position (similar to numpy.unique) EncodedCounts(alphabet=['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y', '*'], counts=array([0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 1, 0, 0, 1, 1, 2, 0, 1, 0, 0]), row_names=None) Counting the number of occurrences of a specific element in each sequence -------------------------------------------------------------------------- >>> np.sum(my_seqs == "F", axis=-1) # count the number of occurrences of the amino acid "F" in each sequence array([0, 0, 1, 0, 0, 1, 0, 0, 1, 1]) Filtering `EncodedRaggedArray` objects based on a mask ------------------------------------------------------ >>> mask = my_seqs.lengths < 8 >>> short_seqs = my_seqs[mask] >>> short_seqs encoded_ragged_array(['CENVC', 'RSTWF', 'NTIFMC'], AlphabetEncoding('ACDEFGHIKLMNPQRSTVWY*')) Broadcasting and one-hot encoding ---------------------------------- >>> short_seqs[1][..., np.newaxis] == "ACDEFGHIKLMNPQRSTVWY" # one-hot encoding of the second sequence array([[False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, False, False, False], [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, False, False], [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, False], [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False], [False, False, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False]]) Reading sequences from file =========================== Usually we get sequences directly from file. BioNumPy supports a range of file formats containing sequence data including fasta, fastq, indexed fasta and bam files. Reading sequence entries (.fa, .fq, .bam) ----------------------------------------- To read in a set of sequence entries, we can just use the `bnp.open` method (here we read a fastq file, but this works the same for fasta and bam files): >>> entries = bnp.open("example_data/reads.fq").read() >>> entries SequenceEntryWithQuality with 2 entries name sequence quality headerishere CTTGTTGA [2 2 2 2 2 2 2 2] anotherheader CGG [93 93 93] We see we get all the entries in the file, with the corresponding fields. The `sequence` field here is an `EncodedRaggedArray` and thus supports numpy-like indexing etc: >>> (entries.sequence=="T").sum(axis=-1) # Count the number of T's array([4, 0]) Reading indexed files (.fa.fai) ------------------------------- When reading a reference genome, we often can't read in the whole file (using `.read()`) and it doesn't make sense to read in the chromosome-sequences as entries in chunks (using `.read_chunks()`) but we rather want to read specific parts of the genome. In these cases we can use an index to be able to read specific regions. We then use `bnp.open_indexed` function: >>> reference_sequence = bnp.open_indexed("example_data/small_genome.fa") >>> reference_sequence["2"][10:20] encoded_array('atattagcca') Functions workin on sequences ============================= A set of functions working on sequences are gathered in the `bnp.seqeunce` module.