pqsfinder: User Guide

In pqsfinder: Identification of potential quadruplex forming sequences

Introduction

The main functionality of the r Biocpkg("pqsfinder") package is to detect DNA and RNA sequence patterns that are likely to fold into an intramolecular G-quadruplex (G4). G4 is a nucleic acid structure that can form as an alternative to the canonical B-DNA. G4s are believed to be involved in regulation of diverse biological processes, such as telomere maintenance, DNA replication, chromatin formation, transcription, recombination or mutation [@maizels_2014; @kejnovsky_etal_2015]. The main idea of our algorithmic approach is based on the fact that G4 structures arise from compact sequence motifs composed of four consecutive and possibly imperfect guanine runs (G-run) interrupted by loops of semi-arbitrary lengths. The algorithm first identifies four consecutive G-run sequences. Subsequently, it examines the potential of such G-runs to form a stable G4 and assigns a corresponding quantitative score to each. Non-overlapping potential quadruplex-forming sequences (PQS) with positive score are then reported.

It is important to note that unlike many other approaches, our algorithm is able to detect sequences responsible for G4s folded from imperfect G-runs containing bulges or mismatches and as such is more sensitive than competing algorithms [^1]. We also believe the presented solution is the most scalable, since it can be easily and quickly customized (see chapter Customizing detection algorithm for details). The program can be made to detect novel or experimental G4 types that might be discovered or studied in future.

[^1]: We have tested r Biocpkg('pqsfinder') on experimentally verified G4 sequences. The results of that work are reflected in default settings of searches. Details of these tests will be presented elsewhere.

For those interested in non-B DNA, we have previously authored a similar package that can be used to search for triplexes, another type of non-B DNA structure. For details, please see r Biocpkg('triplex') package landing page.

G-quadruplex detection

As usual, before first package use, it is necessary to load the r Biocpkg("pqsfinder") package using the following command:

library(pqsfinder)

Identification of potential quadruplex-forming sequences (PQS) in DNA is performed using the pqsfinder function. This function has one required parameter representing the studied DNA sequence in the form of a DNAString object and several modifying options with predefined values. For complete description, please see pqsfinder function man page.

Basic quadruplex detection

As a simple example, let's find all PQS in a short DNA sequence.

seq <- DNAString("TTTTGGGCGGGAGGAGTGGAGTTTTTAACCCCAAAAATTTGGGAGGGTGGGTGGGAGAA")
pqs <- pqsfinder(seq, min_score = 20)
pqs

Detected PQS are returned in the form of a PQSViews class, which represents the basic container for storing a set of views on the same input sequence based on XStringViews object from r Biocpkg("Biostrings") package. Each PQS in the view is defined by (i) start location, (ii) width, (iii) score, (iv) strand, (v) number of G-tetrads nt, (vi) number of bulges nb and (vii) number of mismatches nm. The first four values can be accessed by standard functions start(x), width(x) and score(x) and strand(x). To get other PQS features, please use elementMetadata(x) function. It additionaly provides run and loop lengths of the detected PQS (rl1, rl2, rl3, ll1, ll2, ll3).

elementMetadata(pqs)

By default, pqsfinder function reports only the locally best non-overlapping PQS, ignoring any other that would overlap it. However, it's possible to change the default behavior by setting the overlapping option to TRUE.

pqsfinder(seq, overlapping = TRUE, min_score = 30)

Alternatively, it's possible to get numbers of all overlapping PQS at each position of the input sequence. To achieve that, set deep option to TRUE and then call density(x) function on the PQSViews object:[^2]

[^2]: Clusters of overlapping PQS usually have steep edges when the number of neighboring G-runs is low, but could be more spread out in other situations.

pqs <- pqsfinder(seq, deep = TRUE, min_score = 20)
density(pqs)

The following example shows, how such density vector could be simply visualized along the input sequence using r Biocpkg('Gviz') from Bioconductor.

library(Gviz)
ss <- DNAStringSet(seq)
names(ss) <- "chr1"
dtrack <- DataTrack(
  start = 1:length(density(pqs)), width = 1, data = density(pqs),
  chromosome = "chr1", genome = "", name = "density")
strack <- SequenceTrack(ss, chromosome = "chr1", name = "sequence")
suppressWarnings(plotTracks(c(dtrack, strack), type = "h"))

Modifying basic algorithm options

Depending on the particular type of PQS you want to detect, the algorithm options can be tuned to find the PQS effectively and exclusively. The table bellow gives an overview of all basic algorithm options and their descriptions.

Option name | Description -----------------|---------------------------------------- strand | Strand specification (+, - or *). overlapping | If true, than overlapping PQS will be reported. max_len | Maximal total length of PQS. min_score | Minimal score of PQS to be reported. The default value 52 shows the best balanced accuracy on human G4 sequencing data [@chambers_2015]. run_min_len | Minimal length of each PQS run (G-run). run_max_len | Maximal length of each PQS run. loop_min_len | Minimal length of each PQS inner loop. loop_max_len | Maximal length of each PQS inner loop. max_bulges | Maximal number of runs containing a bulge. max_mismatches | Maximal number of runs containing a mismatch. max_defects | Maximum number of defects in total (#bulges + #mismatches).

The more you narrow these options in terms of shorter PQS length, narrower run or loop length ranges and lower number of defects, the faster the detection process will be, with a possible loss of sensitivity.

Important note: In each G-run, the algorithm allows at most one type of defect and at least one G-run must be perfect, that means without any defect. Therefore the values of max_bulges, max_mismatches and max_defects must fall into the range from 0 to 3.

Example 1: If you are insterested solely in G-quadruplexes with perfect G-runs, just restrict max_defects to zero:

pqsfinder(seq, max_defects = 0, min_score = 20)

Example 2: In case you don't mind defects in G-runs, but you want to report only high-quality PQS, increase min_score value:

pqsfinder(seq, min_score = 70)

Exporting results

As mentioned above, the results of detection are stored in the PQSViews object. Because the PQSViews class is only an extension of the XStringViews class, all operations applied to the XStringViews object can also be applied to the PQSViews object as well.

Additionaly, PQSViews class supports a conversion mechanism to create GRanges objects. Thus, all detected PQS can be easily transformed into elements of a GRanges object and saved as a GFF3 file, for example.

GRanges conversion and export to GFF3

In this example, the output of the pqsfinder function will be stored in a GRanges object and subsequently exported as a GFF3 file. At first, let's do the conversion using the following command:

gr <- as(pqs, "GRanges")
gr

Please note that the chromosome name is arbitrarily set to chr1, but it can be freely changed to any other value afterwards. In the next step the resulting GRanges object is exported as a GFF3 file.

library(rtracklayer)
export(gr, "test.gff", version = "3")

Please note, that it is necessary to load the rtracklayer library before running the export command. The contents of the resulting GFF3 file are:

text <- readLines("test.gff")
cat(strwrap(text, width = 80, exdent = 3), sep = "\n")

Another possibility of utilizing the results of detection is to transform the PQSViews object into a DNAStringSet object, another commonly used class of the Biostrings package. PQS stored inside DNAStringSet can be exported into a FASTA file, for example.

DNAStringSet conversion and export to FASTA

In this example, the output of the pqsfinder function will be stored in a DNAStringSet object and subsequently exported as a FASTA file. At first, let's do the conversion using the following command:

dss <- as(pqs, "DNAStringSet")
dss

In the next step, the DNAStringSet object is exported as a FASTA file.

writeXStringSet(dss, file = "test.fa", format = "fasta")

The contents of the resulting FASTA file are:

text <- readLines("test.fa")
cat(text, sep = "\n")

Please, note that all attributes of detection such as start position, end position and score value are stored as a name parameter (inside the DNAStringSet), and thus, they are also shown in the header line of the FASTA format (the line with the initial > symbol).

A real world example

In the following example, we load the human genome from the r Biocpkg('BSgenome') package and identify all potential G4 (PQS) in the region of AHNAK gene on chromose 11. We then export the identified positions into a genome annotation track (via a GFF3 file) and an additional FASTA file. Finally, we plot some graphs showing the PQS score distribution and the distribution of PQS along the studied genomic sequence.

1. Load necessary libraries and genomes.

   r library(pqsfinder) library(BSgenome.Hsapiens.UCSC.hg38) library(rtracklayer) library(ggplot2) library(Gviz)

2. Retrive AHNAK gene annotation.

   ```r gnm <- "hg38" gene <- "AHNAK"

   Load cached AHNAK region track:

   load(system.file("extdata", "gtrack_ahnak.RData", package="pqsfinder"))

   Alternatively, query biomaRt API with the following commands:

   library(biomaRt)

   gtrack <- BiomartGeneRegionTrack(genome = gnm, symbol = gene, name = gene)

   ```

3. Get AHNAK sequence from r Biocpkg('BSgenome') package extended by 1000 nucleotides on both sides.

   r extend <- 1000 seq_start <- min(start(gtrack)) - extend seq_end <- max(end(gtrack)) + extend chr <- chromosome(gtrack) seq <- Hsapiens[[chr]][seq_start:seq_end]

4. Search for PQS on both strands.

   r pqs <- pqsfinder(seq, deep = TRUE)

5. Display the results.

   r pqs

6. Sort the results by score to see the best one.

   r pqs_s <- pqs[order(score(pqs), decreasing = TRUE)] pqs_s

7. Export all PQS into a GFF3-formatted file.

   r export(as(pqs, "GRanges"), "test.gff", version = "3")

   The contents of the GFF3 file are as follows (the first three records only):

   r text <- readLines("test.gff", n = 5) cat(strwrap(text, width = 80, exdent = 3), sep= "\n")

8. Export all PQS into a FASTA format file.

   r writeXStringSet(as(pqs, "DNAStringSet"), file = "test.fa", format = "fasta")

   The contents of the FASTA file are as follows (the first three records only):

   r text <- readLines("test.fa", n = 6) cat(text, sep = "\n")

9. Show histogram for score distribution of detected PQS.

   r sf <- data.frame(score = score(pqs)) ggplot(sf) + geom_histogram(mapping = aes(x = score), binwidth = 5)

10. Show PQS score and density distribution along AHNAK gene annotation using r Biocpkg('Gviz') package.

    r strack <- DataTrack( start = start(pqs)+seq_start, end = end(pqs)+seq_start, data = score(pqs), chromosome = chr, genome = gnm, name = "score") dtrack <- DataTrack( start = (seq_start):(seq_start+length(density(pqs))-1), width = 1, data = density(pqs), chromosome = chr, genome = gnm, name = "density") atrack <- GenomeAxisTrack() suppressWarnings(plotTracks(c(gtrack, strack, dtrack, atrack), type = "h"))

    The stacked plot of the score and density distribution might help to assess the singularity of PQS. Higher density values indicates low-complexity regions (full of guanines), so it is expected to contain high-scoring PQS. On the other hand, a high-scoring PQS in low-density region might be an interesting target.

Customizing the detection algorithm

The underlying detection algorithm is almost fully customizable, it can even be set up to find fundamentally different types of G-quadruplexes. The very first option how to change the detection behavior is to tune scoring bonuses, penalizations and factors. Supported options are summarized in the table bellow:

Option name | Description ---------------------|-------------------------------------- tetrad_bonus | G-tetrad bonus, regardless the tetrade contains mismatches or not. mismatch_penalty | Penalization for a mismatch in tetrad. bulge_penalty | Penalization for a bulge. bulge_len_factor | Penalization factor of a bulge length. bulge_len_exponent | Exponent of a bulge length. loop_mean_factor | Penalization factor of a loop length mean. loop_mean_exponent | Exponent of a loop length mean.

Customizing the scoring function

A more complicated way to influence the algorithm output is to implement a custom scoring function and pass it throught the custom_scoring_fn options. Before you start experimenting with this feature, please consider the fact that custom scoring function can influence the overall algorithm performance very negatively, particularly on long sequences. The best use case of this feature is rapid prototyping of novel scoring techniques, which can be later implemented efficiently, for example in the next version of this package. Thus, if you have any suggestions how to further improve the default scoring system (DSS), please let us know, we would highly appreciate that.

Basically, the custom scoring function should take the following 10 arguments:

• subject - input DNAString object,
• score - positive PQS score assigned by DSS, if enabled,
• start - PQS start position,
• width - PQS width,
• loop_1 - loop #1 start position,
• run_2 - run #2 start position,
• loop_2 - loop #2 start position,
• run_3 - run #3 start position,
• loop_3 - loop #3 start position,
• run_4 - run #4 start position.

The function will return a new score as a single integer value. Please note that if use_default_scoring is enabled, the custom scoring function is evaluated after the DSS but only if the DSS resulted in positive score (for performance reasons). On the other hand, when use_default_scoring is disabled, custom scoring function is evaluated on every PQS.

Example: Imagine you would like to assign a particular type of quadruplex a more favourable score. For example, you might want to reflect that G-quadruplexes with all loops containing just a single cytosine tend to be more stable than similar ones with different nucleotide at the same place. This can be easily implemented by the following custom scoring function:

c_loop_bonus <- function(subject, score, start, width, loop_1,
                         run_2, loop_2, run_3, loop_3, run_4) {
  l1 <- run_2 - loop_1
  l2 <- run_3 - loop_2
  l3 <- run_4 - loop_3
  if (l1 == l2 && l1 == l3 && subject[loop_1] == DNAString("C") &&
      subject[loop_1] == subject[loop_2] &&
      subject[loop_1] == subject[loop_3]) {
    score <- score + 20
  }
  return(score)
}

Without the custom scoring function, the two PQS found in the example sequence will have the same score.

seq <- DNAString("GGGCGGGCGGGCGGGAAAAAAAAAAAAAGGGAGGGAGGGAGGG")
pqsfinder(seq)

However, if the custom scoring function presented above is applied, the two PQS are clearly distinguishable by score:

pqsfinder(seq, custom_scoring_fn = c_loop_bonus)

Complete replacement of the default scoring system

There might be use cases when it is undesirable to have the default scoring system (DSS) enabled. In this example we show how to change the detection algorithm behavior to find quite a different type of sequence motif - an interstrand G-quadruplex (isG4) [@kudlicki_2016]. Unlike standard intramolecular G-quadruplex, isG4 can be defined by interleaving runs of guanines and cytosines respectively. Its canonical form can be described by a regular expression G_nN_aC_nN_bG_nN_cC_n.

To detect isG4s by the pqsfinder function, it is essential to change three options. At first, disable the DSS by setting use_default_scoring to FALSE. Second, specify a custom regular expression defining one run of the quadruplex by setting run_re to G{3,6}|C{3,6}. The last step is to define a custom scoring function validating each PQS:

isG4 <- function(subject, score, start, width, loop_1,
                 run_2, loop_2, run_3, loop_3, run_4) {
  r1 <- loop_1 - start
  r2 <- loop_2 - run_2
  r3 <- loop_3 - run_3
  r4 <- start + width - run_4

  if (!(r1 == r2 && r1 == r3 && r1 == r4))
    return(0)

  run_1_s <- subject[start:start+r1-1]
  run_2_s <- subject[run_2:run_2+r2-1]
  run_3_s <- subject[run_3:run_3+r3-1]
  run_4_s <- subject[run_4:run_4+r4-1]

  if (length(grep("^G+$", run_1_s)) && length(grep("^C+$", run_2_s)) &&
      length(grep("^G+$", run_3_s)) && length(grep("^C+$", run_4_s)))
    return(r1 * 20)
  else
    return(0)
}

Let's see how it all works together:

pqsfinder(DNAString("AAAAGGGATCCCTAAGGGGTCCC"), strand = "+",
          use_default_scoring = FALSE, run_re = "G{3,6}|C{3,6}",
          custom_scoring_fn = isG4)

Session info

Here is the output of sessionInfo() on the system on which this document was compiled:

sessionInfo()

References

pqsfinder documentation built on Nov. 8, 2020, 5:52 p.m.