Genomic sequences

The genomic sequences refer to the sequence of the chromosome or fragment where your SNVs (or SNPs) are located. They can be complete chromosome sequences or large fragments extracted from chromosomes. We will predict G4 structures from the sequences you provide and analyze whether these SNVs can affect their formation.

Please note that the coordinates of your SNVs must be relative to the genomic sequence you provide, and they should be in 1-based coordinates.

The genomic sequences need to be processed into a DNAStringSet object. We provide a built-in function, loadSequence, to facilitate this processing. Of course, you can also load your customized sequence file and convert it into a DNAStringSet object without using loadSequence function.

loadSequence mainly accepts three types of input:

• A two-column data.frame, with the first column containing the sequence identifiers and the second column containing their corresponding sequences. This should be specified using the genome_seq parameter.

• The path to a stored FASTA file. The FASTA file must have a .fa, .fna, or .fasta extension. This should be specified using the seq_path parameter.

• A text file (.txt) that stores the sequence identifiers and their corresponding sequences. The first column should contain the sequence identifiers, and the second column should contain the sequences. This should also be specified using the seq_path parameter.
  Note: Please do not include column names in the file!

Here are some examples for you to load sequences into a DNAStringSet object using loadSequence function:

Load from a data.frame object

seq_df <- data.frame(chr = c("seq1", "seq2"),
                     sequence = c(paste0(rep("G", 100), collapse = ""), 
                                  paste0(rep("A", 100), collapse = "")))
seq <- loadSequence(genome_seq = seq_df)

Load from a fasta file

# File path to the sequences in fasta format
fa_path <- system.file("extdata", "seq.fa", package = "G4SNVHunter")
seq <- loadSequence(seq_path = fa_path)

Load from a txt file

# File path to the sequences in txt format
txt_path <- system.file("extdata", "seq.txt", package = "G4SNVHunter")
seq <- loadSequence(seq_path = txt_path)

We can also obtain genome sequences for some common species from Bioconductor Annotation Packages. While this is convenient, it requires you to install some related packages in advance. For example:

# Load sequence for chromosome 21 (hg19)
hg19 <- BSgenome.Hsapiens.UCSC.hg19
chr21_seq <- DNAStringSet(hg19$chr21)
# Chromosome names are needed for analysis
names(chr21_seq) <- "chr21"

SNV data

The SNV (or SNP) data needs to be processed into a GRanges object. Since the form of SNV data is too flexible, we do not provide a function to load SNVs here. However, you can easily load and convert SNV data into a GRanges object yourself. For example,

# Path to SNPs
snp_path <- system.file("extdata", "snp.txt", package = "G4SNVHunter")
# Load SNPs into memory
snp <- read.table(snp_path, sep = "\t", header = FALSE)
# Convert snp to GRanges
snp <- GRanges(seqnames = snp$V1,
               ranges = IRanges(start = snp$V2, width = 1),
               rsid = snp$V3,
               ref = snp$V4,
               alt = snp$V5)
print(snp)

## GRanges object with 15 ranges and 3 metadata columns:
##        seqnames    ranges strand |         rsid         ref         alt
##           <Rle> <IRanges>  <Rle> |  <character> <character> <character>
##    [1]    chr21  15397335      * | rs1244782663           G           A
##    [2]    chr21  39931215      * |  rs910271214           G           A
##    [3]    chr21  22915361      * | rs1210166884           A           G
##    [4]    chr21  42644592      * |  rs554035836           A           G
##    [5]    chr21  10743898      * |  rs955290928           A           G
##    ...      ...       ...    ... .          ...         ...         ...
##   [11]    chr21  30836340      * | rs1212424980           T           C
##   [12]    chr21  36426864      * | rs1478442431           G           T
##   [13]    chr21  27584238      * | rs1191559536           T           C
##   [14]    chr21  36106373      * |  rs561521189           C           T
##   [15]    chr21  21840835      * | rs1301205088           T           C
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

SNV data validation

We provide a built-in checkSNV function for you to verify whether your SNV (or SNP) data meets the requirements of our program.

However, it is not necessary to run this function explicitly, as our program will automatically perform the checks at the appropriate stages. Nevertheless, performing a proactive check is always a good practice.

gr <- GRanges(
  seqnames = Rle("seq1"),
  ranges = IRanges(c(100, 200, 300), width=1),
    ref = c("A", "C", "G"),
    alt = c("T", "T", "A")
  )
# Check width ('w'), ref ('r'), and alt ('a')
# Returns TRUE if it passed the validation
checkSNV(gr, mode = "wra", ref_col = "ref", alt_col = "alt")

## [1] TRUE

Leverage G4Hunter for G4 detection

You can use the G4HunterDetect function to directly predict G4s from a DNAStringSet object containing sequences based on G4Hunter algorithm.
For example,

# Sequence file in fasta file format 
fa_path <- system.file("extdata", "seq.fa", package = "G4SNVHunter")
# Load sequences
seq <- loadSequence(seq_path = fa_path)

# Predict G4s
G4_detected <- G4HunterDetect(seq)

## Start processing...

## Merging predicted G4s...

## Done!

Then, we can examine the prediction results, which are stored as a GRanges object.

You can use functions like print to directly view the results. However, in this instance, we’ll leverage the datatable function to view the predicted G4s, as this method provides a more user-friendly display.

Let’s take a look at the G4s,

datatable(as.data.frame(G4_detected), options = list(scrollX = TRUE))

This function will return a GRanges object containing all potential G4s predicted by the G4Hunter algorithms. It mainly includes:

• The sequence identifier where the G4 is located (seqnames),

• The position of the G4 relative to the sequence you provided (ranges, 1-based coordinate),

• The strand on which the G4 is located (strand), with + indicating the G4 is on the positive strand, and - indicating it is on the negative strand,

• The G4Hunter score for the G4 (score),

• The maximum G4Hunter score of the windows covering that G4 during the execution of the G4Hunter algorithm (max_score, used to determine if the G4 surpasses the threshold),

• The predicted G4 sequence on the positive strand (sequence, with C-rich means G4s on the reverse strand).

Please note that score reflects the precise score of the G4 sequence, while max_score represents the highest score obtained by sliding a window across that G4.

See the illustration below.

Some parameters for prediction can be customized. For example,

# Predict G4 by customizing other parameters
G4_detected <- G4HunterDetect(seq, threshold = 1.5, window_size = 20)

## Start processing...

## Merging predicted G4s...

## Done!

• threshold: G4Hunter will search for G4s in windows above this threshold (absolute value). Default is 1.5. Unless there are special needs, we do not recommend setting the threshold below 1.2.

• window_size: The window size (bp) for G4Hunter prediction. Default is 25. Another commonly used window size is 20. However, 25 is generally preferred unless otherwise required.

• include_sequences: Whether to include the predicted G4 sequences in the output. Default is TRUE.

• strands: Indicate which strand requires prediction, with b for both strands and p for positive strand only. Please note that if your genome is single-stranded, this parameter should be set to p to prevent the program from predicting G4s on a non-existent negative strand.

We generally do not recommend modifying certain parameters, such as window_size and threshold, as their default settings are already optimal.

Export prediction results

Since we store the prediction results as a GRanges object, exporting is very straightforward.

For example, you can export it into a CSV file,

# export as csv format
write.csv(as.data.frame(G4_detected), "/path/to/your/file.csv")

or export it into other formats as well.

# export as bed format
export(G4_detected, "/path/to/your/file.bed", format = "bed")
# export as bigwig format
export(G4_detected, "/path/to/your/file.bw", format = "bigWig")

Visualize G4Hunter prediction results

You can use plotG4Info function for visualizing the statistical information of G4s predicted by the G4HunterDetect function.

plotG4Info(G4_detected)

The left panels (A, C, E) show the overall distribution of the max_score (absolute value), score (absolute value), and length of G4s, respectively. The right panels (B, D, F) illustrate the distribution of max_score, score, and length for G4 structures on both the positive and negative strands. The scores for G4s on the positive strand are positive, while those on the negative strand are negative, resulting in a bimodal distribution of max_score and score. Please note that the positive or negative score reflects only the type of strand the G4 is on, whereas the capability of the G4 structure to form is determined by the absolute value of the score.

Variant-centric assessment

The default mode of SNVImpactG4 is variant-centric (mode = 's'), which requires only the G4 data in GRanges format (as returned by the G4HunterDetect function), the SNV data in GRanges format, and the column name for the alternate bases.

Then, the assessment can be done easily by,

snp_eff <- SNVImpactG4(chr21_G4, 
                       snp_gr, 
                       alt_col = "alt")

Let’s view the first three records,

datatable(as.data.frame(snp_eff[1:3]), options = list(scrollX = TRUE))

SNVImpactG4 function will return a GRanges object that includes detailed information about the SNVs (core columns and SNV.info.* columns), their associated G4 regions (G4.info.* columns), and the changes in G4 sequences (mut.G4.seq and mut.G4.anno.seq columns) and G4Hunter scores (mut.score and score.diff column).

Please note that SNVs located outside of G4 regions will be automatically ignored.

Sample-centric assessment

You can set the mode parameter in SNVImpactG4 function to 'm' to enable sample-centric mode. This mode is particularly useful for specific scenarios, such as analyzing SNVs derived from cancer patients.

When using this mode, you must also specify the column names for sample IDs (sampleid_col) and SNV IDs (snvid_col).

# Column names of the Sample ID and SNV ID must be specified
snv_eff <- SNVImpactG4(chr21_G4, 
                       snv_gr, 
                       alt_col = "alt", 
                       mode = "m", 
                       sampleid_col = "sampleid", 
                       snvid_col = "snv_id")

Let’s view the first three records,

datatable(as.data.frame(snv_eff[1:3]), options = list(scrollX = TRUE))

It will return a GRanges object containing detailed information about the G4s (core columns and G4.info.* columns). The object will also include the IDs of the SNVs (snv.ids column) located in the G4 regions, the sample IDs (sample.ids column), the mutated G4 sequences (mut.G4.seq and mut.G4.anno.seq columns), and the changes in G4Hunter scores (mut.score and score.diff columns).

Please note that in sample-centric mode, when multiple SNVs from the same sample occur within a single G4 region, we will consider the combined impact of these SNVs. For example, the following G4 overlaps with two SNVs (id_3608 and id_49857) from sample GQ3Y, so the final G4Hunter score is based on the combined effect of these two SNVs.

datatable(as.data.frame(snv_eff[528]), options = list(scrollX = TRUE))

Filtering potentially affected G4 sequences

Given that some SNVs may have minimal effects on G4 formation, we need to filter out those SNVs that are worth investigating further for experimental validation or additional analysis. In other words, identifying those variants that have the potential to disrupt the G4 structures.

We provide a convenient function, filterSNVImpact, to help with this filtering process.

There are three threshold parameters for you to adjust: raw_score_threshold, mut_score_threshold, and score_diff_threshold.

If raw_score_threshold (a positive number, but no more than 4) is specified, filterSNVImpact will filter out records where the absolute value of the original G4Hunter score is below this threshold.

If mut_score_threshold (a positive number, but no more than 4) is specified, filterSNVImpact will retain records where the G4Hunter score of the mutated G4 sequences does not exceed this threshold.

If score_diff_threshold (a negative number, but no less than -4) is specified, filterSNVImpact will retain records where the decrease in the G4Hunter score after mutation exceeds this threshold. For example, if score_diff_threshold = -0.2, those records will be retained: |G4HunterScore_{mut_seq}| − |G4HunterScore_{raw_seq}| ≤ −0.2

Our recommendation is that raw_score_threshold should be greater than 1.5, and the mut_score_threshold should be less than 1.2. This is an empirically based guideline, and you are certainly free to adjust them to be more stringent or to use these parameters for flexible customization as needed.

Please note that you must specify at least one threshold parameter, but you do not need to specify all of them.

For example, using the results returned in mode m

filtered_snv_eff <- filterSNVImpact(snv_eff,
                                    mut_score_threshold = 1.2,
                                    score_diff_threshold = -0.35)

We can examine the filtered records to identify SNVs that have a substantial impact on G4 formation.

datatable(as.data.frame(filtered_snv_eff), options = list(scrollX = TRUE))

The final results can be saved to a file similarly to how G4 prediction results are saved, as described in Section 5.2.

Visualize changes in G4Hunter scores

You can use the plotSNVImpact function we provide to visualize changes in G4Hunter scores. This function is suitable for both variant-centric and sample-centric output data, as well as their filtered outputs.

We can observe that mutated G4 structures generally exhibit a trend of decreased formation potential.

plotSNVImpact(snv_eff)

Let’s examine how the G4Hunter scores have changed in the filtered_snv_eff object.

We observe that the scores for the filtered G4s have decreased significantly.

plotSNVImpact(filtered_snv_eff)

Plot the seqlogo for mutated G4 sequences

Finally, we can use plotImpactSeq to visualize the G4 sequences affected by SNVs.

For example, we can plot the sequence logos for the top five G4s with the most significant changes in their G4Hunter scores.

top5_snv_eff <- filtered_snv_eff[order(filtered_snv_eff$score.diff)[1:5]]
plotImpactSeq(top5_snv_eff, ncol = 2)

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the ggseqlogo package.
##   Please report the issue at <https://github.com/omarwagih/ggseqlogo/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.