High throughput sequencing is most widely used as it saves a lot of time and provide good results, and produces a huge amount of data which is difficult to manage and especially the tasks and operations performed on it are also very difficult. To ease this purpose, a Python framework have been introduced by Simon Anders and team members, this framework is known as “HTSeq”.HTSeq is a Python library which easily develops the scripts required to fulfill a particular task on the HT data. Basically,HTSeq reads various formats and break it down into recognized strings of characters for further analysis. It also consists of different classes genomic coordinates, sequences, sequencing reads, alignments, gene model information, etc.
Two stand-alone applications have also been developed along with HTSeq, namely, htseq-qa for read quality assessment and htseq-count for preprocessing RNA-Seq alignments for analyzing differential expression.
HTSeq can read various formats such as FASTA, FASTQ (short reads), SAM/BAM (short-read alignments). Wherever appropriate, different parsers will yield the same type of record objects. For example, the record class SequenceWithQualities is used whenever sequencing read with base-call qualities needs to be presented, and hence yielded by the FastqParser class and also present as a field in the SAM_Alignment objects yielded by SAM_Reader or BAM_Reader parser objects (Fig. 1). There are some specific classes to represent Genomic Position and Genomic Intervals of the sequence. In order to achieve good performance, various parts of HTSeq is written in ‘Cython’ ( a tool which translates Python code augmented with C).
Fig. 1. ( a) The SAM_Alignment class as an example of an HTSeq data record: subsets of the content are bundled in object-valued fields, using classes (here SequenceWithQualities and GenomicInterval) that are also used in other data records to provide a common view on diverse data types. ( b) The cigar field in a SAM_alignment object presents the detailed structure of a read alignment as a list of CigarOperation. This allows for convenient downstream processing of complicated alignment structures, such as the one given by the cigar string on top and illustrated in the middle. Five CigarOperation objects, with slots for the columns of the table (bottom) provide the data from the cigar string, along with the inferred coordinates of the affected regions in read (‘query’) and reference.
HTSeq also consists of a class which deals with the gapped alignments, namely SAM_Alignment, with multiple alignments and with paired-end data. HTSeq provides a function, pair_SAM_alignments_with_buffer, to pair up the alignment records by keeping a buffer of reads whose end pair has not yet been found, and so facilitates processing data on the level of sequenced fragments rather than reads. HTSeq also facilitates the storage of genome-position-dependent data, which means that each base pair position on the genome can be given a particular value that can be easily stored and retrieved by simply entering the same value.
The script htseq-qa is a simple tool for initial quality assessment of sequencing runs. It produces plots that summarize the nucleotide compositions of the positions in the read and the base-call qualities. As we discussed earlier in this article that htseq-count is a tool for RNA-Seq data analysis. It counts for each gene that how many aligned reads overlap the sequence exons. Since it is designed specifically to analyse differential expression only reads mapping unambiguously to a single gene are considered and the reads aligned to multiple positions or overlapping with more than one gene are discarded. In case of paired-end data, htseq-count counts only the fragment not the reads because the two paired ends originating from the same fragment provide only evidence for one cDNA fragment and should hence be counted only once.
In this way, HTSeq offers a comprehensive solution to facilitate a wide range of programming tasks in HTS data analysis. For further reading, click here.
An exhaustive list of references for this article is available with the author and is available on personal request, for more details write to [email protected]
How to create a pie chart using Python?
How to make swarm boxplot?
With the new year, we are going to start with a very simple yet complicated topic (for beginners) in bioinformatics. In this tutorial, we provide a simple code to plot swarm boxplot using matplotlib and seaborn. (more…)
How to obtain ligand structures in PDB format from PDB ligand IDs?
How to obtain SMILES of ligands using PDB ligand IDs?
Fetching SMILE strings for a given number of SDF files of chemical compounds is not such a trivial task. We can quickly obtain them using RDKit or OpenBabel. But what if you don’t have SDF files of ligands in the first place? All you have is Ligand IDs from PDB. If they are a few then you can think of downloading SDF files manually but still, it seems time-consuming, especially when you have multiple compounds to work with. Therefore, we provide a Python script that will read all Ligand IDs and fetch their SDF files, and will finally convert them into SMILE strings. (more…)
How to get secondary structure of multiple PDB files using DSSP in Python?
vs_analysis_compound.py: Python script to search for binding affinities based on compound names.
How to download files from an FTP server using Python?
How to convert the PDB file to PSF format?
smitostr.py: Python script to convert SMILES to structures.
How to preprocess data for clustering in MATLAB?
Data preprocessing is a foremost and essential step in clustering based on machine learning methods. It removes noise and provides better results. In this article, we are going to discuss the steps involved in data preprocessing using MATLAB . (more…)
How to calculate drug-likeness using RDKit?
RDKit  allows performing multiple functions on chemical compounds. One is the quantitative estimation of drug-likeness also known as QED properties. These properties include molecular weight (MW), octanol-water partition coefficient (ALOGP), number of hydrogen bond donors (HBD), number of hydrogen bond acceptors (HBA), polar surface area (PSA), number of rotatable bonds (ROTB), number of aromatic rings (AROM), structural alerts (ALERTS). (more…)
sdftosmi.py: Convert multiple ligands/compounds in SDF format to SMILES.
tanimoto_similarities_one_vs_all.py – Python script to calculate Tanimoto Similarities of multiple compounds
We previously provided a Python script to calculate the Tanimoto similarities of multiple compounds against each other. In this article, we are providing another Python script to calculate the Tanimoto similarities of one compound with multiple compounds. (more…)
tanimoto_similarities.py: A Python script to calculate Tanimoto similarities of multiple compounds using RDKit.
RDKit  is a very nice cheminformatics software. It allows us to perform a wide range of operations on chemical compounds/ ligands. We have provided a Python script to perform fingerprinting using Tanimoto similarity on multiple compounds using RDKit. (more…)
How to commit changes to GitHub repository using vs code?
Extracting first and last residue from helix file in DSSP format.
How to extract x,y,z coordinates of atoms from PDB file?
The x, y, and z coordinates of atoms are provided in the PDB file. One way to extract them is by using the Biopython package . In this article, we will extract coordinates of C-alpha atoms for each residue from the PDB file using Biopython. (more…)
dssp_parser: A new Python package to extract helices from DSSP files.
How to calculate center of mass of a protein structure using Python script?
How to sort binding affinities based on a cutoff using vs_analysis.py script?
Previously, we have provided a Python script (vs_analysis.py) to analyze the virtual screening (VS) results of Autodock Vina. Now, we have updated this script to sort binding affinities based on user inputted cutoff value. (more…)