SNPector: SNP inspection tool for diagnosing gene pathogenicity and drug response in a naked sequence

Functions

The SNPector Python tool uses many packages to inspect the existence of SNPs in a given sequence. Moreover, SNPector provides users with detailed visualization figures, highlighting other SNPs with similar mutation effects on protein phosphorylation, ubiquitination, methylation, or sumoylation sites, and predicts substrates of N-acetyltransferase.

Additionally, SNPector provides the ability to visualize obtained information about the linkage disequilibrium of detected SNPs using various Python packages, such as Matplotlib¹⁷, generating a number of figures that summarize vast amounts of previously published data indicating SNPs allelic segregation, association, minor allele frequency. Figure 3 shows an example of illustrations that can be generated through SNPector.

Operation

SNPector requires at least Python 3.5, 16GB RAM, i7 cores, and 8 MB.

Workflow

SNPector was written using Pythpn3 programming language as a standalone package and can be run on different operating systems platforms supported with Python 3.x compilers. To achieve user-friendly usage, the SNPector only accepts input from FASTA sequence (Figure 1) and can be operated from a console through simple command line (Figure 2).

Figure 1. Example input naked FASTA sequence.

Figure 2. SNPector command line parameters.

Figure 3. Collective figure to show all Illustrations provided by SNPector.

(A) Circos illustration where other SNPs that have same proprieties are located. (B) Lollipop plot shows values by vertical columns (C) Counter Plot between two values creating a different coloured shade in which more contrast means higher value. (D) Numerical schematic showing the distribution between four values by plotting and scaling colour contrast according to other to values. (E) Heat map between SNP linkage disequilibrium matrix to show how two SNPs are linked. (F) Marginal plot combining column graph and plot, both showing the relationship between two values. (G) Dendrogram with heat-map showing how all SNP are linked to each other. (H) Histogram with box plot to compare visually between two values. (I) Plot illustrating the regression fit of two plotted value. (J) 3D plot of three values. (K) Annotated heat-map showing the plotted value.

SNPector uses different SNP record information collected from NCBI ClinVar (159,184 records), AWESOME (1,080,551 records), and PharmGKB¹⁶ (3,932 records). Ldlink is an online tool that can be used to assess linkage imbalance (LD) throughout ancestral populations and is a popular method of exploring population-specific genetic framework and functionally navigation disease susceptibility areas²¹. In SNPector, an Application Program Interface (API) has been programmed to download an LDhap file containing linkage disequilibrium statistics and potentially functional variants for a query variant resulting from the inputted FASTA sequence.

SNPector starts by running BLAST²² software locally to find out the genomic location of a given DNA sequence on human genome. If successful, it retrieves the SNP records located within the query genomic range using NCBI ClinVar. According to retrieved records from the database, the detected SNPs in user-provided queries are marked as wild or mutated. Additionally, more information regarding detecting SNPs records will be retrieved from different implemented databases. This information will be used to generate different illustration figures.

If the process is successfully finished, SNPector will generate four different files: (A) Text file containing the output BLAST result, where the genomic location of the user-defined sequences is predicted; (B) tab delimited file containing SNPs retrieved by NCBI ClinVar located in the same regions; (C) two files regarding specific SNPs information retrieved from AWESOME and PharmGKB databases; (D) different figures depicting SNPs with a similar mutation effect to the detected SNPs located on other genomic regions, SNP linkage disequilibrium, the relationship between SNP, drug, and phenotype (Figure 3).

To achive maximum user-friendly usage, SNPector can be run and controlled by command line. SNPector command line structure (Figure 2) is as follows: A) Python3 compiler; B) scan_dna.py: program main script that contain all functions; C) -blaston / -blastoff: in order to initiate BLAST process to provide sequence alignment against the genome to locate where the sequence is situated, if the -blastoff is chosen it will use previous BLAST results; D) -modescan: to scan the given sequence and find out whether SNP occurs or exists in sequence or not, and -modesearch: to extract all SNPs occur in this range of sequence regardless they are exist or not; E) -circoson: draws a Circos figure to illustrate where SNPs with same properties/effect are located; F) -networkon: in order to link between SNPs, diseases and drugs and produces network HTML file; G) -download: activates the API to download data for identified SNPs from LDlink database; H) -vis: produces different figures and plots; I) GivenSequence.fasta: Tte user-provided sequence in FASTA file format. Any of the previous parameters can be deactivated when replaced with -off.

Data

We use part of an EGFR gene sequence downloaded from NCBI nucleotide database in FASTA format as shown in Extended data: File 1²³. The EGFR gene FASTA sequence submitted to NCBI contains SNPs that have a clinical effect involved in disease development, such as breast cancer. SNPector uses: (i) NCBI ClinVar database that describes SNP chromosome, position, ID, reference nucleotide, alternative nucleotide, quality, filter, and information to compare and detect SNPs in EGFR sequence that has clinical complications; (ii) PharmGKB database to investigate the SNP effect on disease development and drug response; (iii) AWESOME database to explore SNP effect on phosphorylation, ubiquitination, methylation, and sumoylation sites; and (iv) Ldlink API database of SNP linkage disequilibrium to find out how detected SNPs are linked to other SNPs.

Loading SNPector libraries

SNPector uses different libraries to import, read, and read data and results. os library is used to run BLAST bash script:

import os                                                                                       

time is used to calculate the time that program.

import time                                                                                     

re refers to regular expression. This library sorts and splits input data with function re.split().

import re                                                                                       

itemgetter module is used to sort BLAST data according to identity, mismatch, and p-value.

from operator import itemgetter                                                                 

From sys library we used sys.argv[] to convert script to command line, which can be run and controlled from the terminal.

import sys                                                                                      

Then we import the Scripts package to visualize the data as follows:

from Scripts.Circos import DrawCircos                                                                              
from Scripts.Network import DrawNetwork                                                                            
from Scripts.Run_BLAST import RunBLAST                                                                             
from Scripts.Extraction import ExtractSNP                                                                          
from Scripts.APIcommands import APIcommands                                                                        
from Scripts.Visualizations import visualization                                                                   
from Scripts.DataVisualization.DownloadWithAPI.LDmatrix import LDmatrix                                            
from Scripts.DataVisualization.DownloadWithAPI.LDhap import LDhap                                                  
from Scripts.DataVisualization.DownloadWithAPI.LDproxy import LDproxy                                              
from Scripts.DataVisualization.CompleteScripts.Ready.ContourPlotWithSeaborn import CounterPlot                     
from Scripts.DataVisualization.CompleteScripts.Ready.CustomLinearRegressionFitSeaborn import LinearReg             
from Scripts.DataVisualization.CompleteScripts.Ready.CustomLollipopPlot import Lollipop                            
from Scripts.DataVisualization.CompleteScripts.Ready.DendrogramWithHeatmapAndColouredLeaves import DendoWithHeatMap
from Scripts.DataVisualization.CompleteScripts.Ready.DensityPlotWithMatplotlib import DenistyPlot                  
from Scripts.DataVisualization.CompleteScripts.Ready.HistogramWithBoxPlot import HistWithBoxPlot                   
from Scripts.DataVisualization.CompleteScripts.Ready.MarginalPlotWithSeaborn import MarginalPlot                   
from Scripts.DataVisualization.CompleteScripts.Ready.ThreeDscatterplot import ThreeDimPlot                         
from Scripts.DataVisualization.CompleteScripts.Ready.UseNormalizationOnSeabornHeatmap import SeabornHeatMap        
from Scripts.DataVisualization.CompleteScripts.Ready.AnnotatedHeatMap import AnnoHeatMap                           
from Scripts.DataVisualization.CompleteScripts.Ready.NumericalSemantics import NumSChem                            
from Scripts.DataVisualization.CompleteScripts.Ready.ThreeDscatterplot import ThreeDimPlot                         
from Scripts.DataVisualization.CompleteScripts.Ready.VolcanoLD import VolLD                                        

SNPector building blocks

RunBLAST() takes the file path of FASTA sequence and starts to align the sequence against the human genome and writes the results to BLAST_RESULT.txt (Extended data: File 2²³).

os.system('./Scripts/blastn -query GivenSequence.fasta -db ./Data/Hum_Genom38 -outfmt 6  - out ./RESULTS/BLAST_RESULT.txt')                                                                

ExtractSNP() reads BLAST_RESULT.txt and sort its with itemgetter() according to the identity, length and p-value, then stores the start and end input given sequence (the query) and subject to use later in the extraction step. It also reads the input FASTA sequence file and stores the sequence variable to use in the comparing extraction step. SNPector provides two inspection modes that can be determined from the terminal, Search and Scan. If mode was “-modesearch”, then SNPector begins to extract all SNPs within the query start and end regardless of their existence in the query. In the mode “-modescan”, SNPector will extract only SNPs that exist in the query

SNPector begins to obtain the alternative nucleotides of SNP through the input sequence and obtains the nucleotides that range from SNP position in ClinVar minus the end position of the subject, to the SNP position in ClinVar minus the end position of the subject plus alternative SNP length to ensure the capture of SNPs from the given sequence and also to detect variants with length more than one nucleotide, finally storing it in the “query_nuc_alt” variable.

query_nuc_alt = sequence[int(snp_pos) - int(subject_end):int(snp_pos) - int(subject_end) + len(snp_alt)]                                                                                   

After the process of extraction the result saved to: “FromAwesom.tsv” file (Extended data: File 3²³), in which SNPector list all other SNPs that have the same effect in different sites in proteins; “FromNCBI.tsv” (Extended data: File 4²³), which is list of the SNPs that SNPector detects in a given sequence and retrieves from NCBI ClinVar Dataset; “FromPharmGKB.tsv” (Extended data: File 5²³), which lista the effect of SNPs in disease development and drug response.

APIcommands() imports SNP IDs from “FromNCBI.tsv” and uses Ldlink API to download “LDhap.csv” file (Extended data: File 6²³), which describes the allele frequency of extracted SNPs, “LDmatrix.csv” file (Extended data: File 7²³), which shows how far detected SNPs are linked to other SNPs, and a file titles with the SNP id (e.g. rs516316.csv) (Extended data: File 8²³), which includes additional information, such as minor allele frequency, linkage disequilibrium and distance of other SNPs linked to the detected SNP.

LDmatrix('./RESULTS/FromNCBI.tsv')                                                              
												
LDhap('./RESULTS/FromNCBI.tsv')                                                                 
												
LDproxy('./RESULTS/FromNCBI.tsv')                                                                

DrawCircos() uses SNP properties from “FromAwesom.tsv” file (Extended data: File 3²³) and searches for other SNPs that have the same properties. SNPector then imports pycircos package to draw SNP location on Circos (Figure 3A).

import pycircos                                                                                 

DrawNetwork() draw network using “FromPharmGKB.tsv” (Extended data: File 5²³) to get the gene name (e.g. EGFR) and by gene name get all SNP that occur in this gene. Using SNP IDs, SNPector obtains disease names caused by these SNPs, and with the disease name SNPector can extract drugs used in treatment for this disease. Finally, with the drug name SNPector can obtain the clinical annotation of the drug. SNPector uses webweb package to draw the network and export it to .html file (Extended data: File 10).

from webweb import web                                                                          
edge_list = Network                                                                             
Web(edge_list).save("./RESULTS/%sVarPhenoDrugNetwork.html" % GeneName)                           

Visualization() uses data downloaded in the “LDmatrix.csv” file (Extended data: File 7²³), and the SNP ID file (e.g. rs516316.csv) (Extended data: File 8²¹) to draw other figures (Figures 3B–K).

Underlying data

Homo sapiens chromosome 7, GRCh38.p13 Primary Assembly, Accession number NC_000007.14: https://www.ncbi.nlm.nih.gov/nuccore/NC_000007.14?report=fasta&from=55019017&to=55211628

Extended data

Zenodo: SNPector Supplementary Data, http://doi.org/10.5281/zenodo.3569790²³.

This project contains the following extended data:

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).