/",
urlfile = "./results/IGVurl.txt"
)
},
step_name = "bam_urls",
dependency = "bwa_alignment",
run_step = "optional"
)
```
### Variant calling with GATK
This section executes variant calling utilizing GATK, with configurations adapted
from the GATK Best Practices guidelines [here](https://software.broadinstitute.org/gatk/best-practices/).
The process includes ten
integrated substeps, with all corresponding command-line parameter files
located in the param/cwl/gatk directory. It is important to note that Base
Quality Score Recalibration (BQSR) and Variant Quality Score Recalibration
(VQSR) are not included in this default workflow configuration. These steps
rely heavily on high-quality, species-specific known variant resources (often
limited to model organisms like humans); their exclusion ensures the workflow
remains broadly applicable to a diverse range of species.
### Step1: fastq to ubam
This step converts raw fastq files into unmapped bam (ubam) format. This
process attaches essential Read Group (RG) metadata - such as sample ID and
sequencing platform - directly to the reads while preserving original base
qualities, which is critical for GATK's downstream error modeling. It is
vital to accurately specify the sequencing platform, which defaults to
`illumina`. If a different platform is utilized, users must modify the
configuration file `param/cwl/gatk/gatk_fastq2ubam.cwl` accordingly to ensure
the variant caller applies the correct parameter adjustments in later stages.
```{r fastq2ubam, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "fastq2ubam",
targets = "targetsPE_varseq.txt",
wf_file = "gatk/workflow_gatk_fastq2ubam.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
inputvars = c(
FileName1 = "_FASTQ_PATH1_",
FileName2 = "_FASTQ_PATH2_",
SampleName = "_SampleName_"
),
dependency = c("faidx_index")
)
```
### Step2: Merge bam and ubam
This step merges the aligned bam file (generated by the aligner) with the
unmapped ubam file (from Step 1) to create a unified bam file. This process
restores critical metadata and original sequence information - such as original
base quality scores and read group tags - that are often discarded or altered
by aligners like BWA. Retaining this "removed" information is essential for
accurate downstream variant statistics and error modeling. While adhering to
these data-merging steps is highly recommended for optimal results, variant
calling can technically be performed without them, albeit with potentially
reduced accuracy.
```{r merge_bam, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "merge_bam",
targets = c("bwa_alignment", "fastq2ubam"),
wf_file = "gatk/workflow_gatk_mergebams.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
inputvars = c(
bwa_men_sam = "_bwasam_",
ubam = "_ubam_",
SampleName = "_SampleName_"
),
rm_targets_col = c("preprocessReads_1", "preprocessReads_2"),
dependency = c("bwa_alignment", "fastq2ubam")
)
```
### Step3: Sort bam files by genomic coordinates
This step sorts the aligned bam files according to their genomic coordinates. Coordinate sorting is a mandatory prerequisite for many downstream analysis
tools, including `MarkDuplicates` and `HaplotypeCaller`. These tools require
reads to be ordered sequentially along the reference genome to efficiently
identify duplicates and call variants.
```{r sort, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "sort",
targets = "merge_bam",
wf_file = "gatk/workflow_gatk_sort.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
inputvars = c(merge_bam = "_mergebam_", SampleName = "_SampleName_"),
rm_targets_col = c(
"bwa_men_sam", "ubam", "SampleName_fastq2ubam",
"Factor_fastq2ubam", "SampleLong_fastq2ubam",
"Experiment_fastq2ubam", "Date_fastq2ubam"
),
dependency = c("merge_bam")
)
```
### Step4: Mark duplicates
This step identifies and marks PCR duplicates introduced during
library preparation and sequencing. The duplicate information (see
`duplicate_metrics` file) is not used in the subsequent steps. It is
generated here to assess the duplication rate and complexity of the
sequencing libraries.
```{r mark_dup, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "mark_dup",
targets = "sort",
wf_file = "gatk/workflow_gatk_markduplicates.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
inputvars = c(sort_bam = "_sort_", SampleName = "_SampleName_"),
rm_targets_col = c("merge_bam"),
dependency = c("sort")
)
```
### Step5: Fixing tags
This step calculates NM, MD, and UQ tags for the BAM files. They provide the
edit distance to the reference (NM), the string encoding of mismatched bases
(MD), and the Phred likelihood of alignment uniqueness (UQ). The tags can be
used for precise variant calling and downstream filtering. While highly
recommended to ensure the alignment file contains standard metadata required
by many analysis tools, this step is technically optional and can be skipped
if the specific variant calling pipeline does not mandate these
pre-calculated metrics.
```{r fix_tag, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "fix_tag",
targets = "mark_dup",
wf_file = "gatk/workflow_gatk_fixtag.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
inputvars = c(mark_bam = "_mark_", SampleName = "_SampleName_"),
rm_targets_col = c("sort_bam"),
dependency = c("mark_dup")
)
```
This completes the sample preprocessing phase. All necessary
analysis-ready BAM files and their corresponding .bai indix files have been
generated. The workflow now proceeds to the variant discovery stage, where
`HaplotypeCaller` is utilized for both individual and cohort-based variant
calling.
### Step6: HaplotypeCaller
In this step, the HaplotypeCaller is executed in gVCF mode . The "g" denotes
"genomic," indicating that the generated file
contains records for all sites, including non-variant blocks, rather than
just the variant sites. By retaining information on homozygous reference
sites, the gVCF format provides the essential data required for the
subsequent cohort calling step. This allows the joint genotyping tool to
accurately validate variants across the population by distinguishing between
"no variant present" and "no data available."
```{r hap_caller, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "hap_caller",
targets = "fix_tag",
wf_file = "gatk/workflow_gatk_haplotypecaller.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
inputvars = c(fixtag_bam = "_fixed_", SampleName = "_SampleName_"),
rm_targets_col = c("mark_bam"),
dependency = c("fix_tag")
)
```
### Step7: Import all gvcfs
To facilitate efficient cohort-wide variant calling in the subsequent step,
it is highly recommended to consolidate all single-sample gVCF files into a
GenomicsDB workspace (see [TileDB](https://github.com/Intel-HLS/GenomicsDB/wiki)).
This database, built upon the TileDB array framework, optimizes data storage
and access speeds for large-scale analyses.
Note: For analyses involving non-diploid data or specific scenarios where
GenomicsDB is not suitable, users should utilize the CombineGVCFs function
from GATK as an alternative . If this method is chosen, the `gvcf_db_folder`
parameter in the `param/cwl/gatk/gatk.yaml` configuration file must be updated
to point to the resulting combined gVCF file path.
Important: Before initiating this step, ensure that all sample `.g.vcf.gz`
files and their corresponding tabix index (.tbi) files are successfully
generated and located within the results directory.
```{r import, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "import",
targets = NULL, dir = FALSE,
wf_file = "gatk/workflow_gatk_genomicsDBImport.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
dependency = c("hap_caller")
)
```
### Step8: Joint genotyping (cohort calling)
This step performs joint genotyping — commonly referred to as cohort
calling - utilizing the GenotypeGVCFs tool. By evaluating genomic
information across the entire cohort simultaneously (accessed via the
GenomicsDB or combined gVCF from the previous step), the algorithm leverages
population-wide statistics to resolve ambiguous sites and improve calling
accuracy. The output is a unified, multi-sample VCF file (default name:
samples.vcf.gz) containing the raw variant calls for all samples, which
serves as the input for subsequent filtering and annotation.
```{r call_variants, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "call_variants",
targets = NULL, dir = FALSE,
wf_file = "gatk/workflow_gatk_genotypeGVCFs.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
dependency = c("import")
)
```
### Step9: Hard filter (cohort) variants
Since Variant Quality Score Recalibration (VQSR) requires large datasets and
highly curated known-variant training resources—which are often unavailable
for non-model organisms. To be applicable to most organisms, this workflow
employs _hard filtering_ as the standard alternative.
This step applies fixed thresholds to specific variant annotations (e.g.,
QualByDepth, FisherStrand, ReadPosRankSum) to identify and remove
low-confidence calls. The default filtering parameters are configured based
on GATK's [generic hard-filtering recommendations](https://gatkforums.broadinstitute.org/gatk/discussion/2806/howto-apply-hard-filters-to-a-call-set).
Users are encouraged to review these settings and adjust the corresponding parameter files to match
the specific error profiles of their dataset. For more information on the
requirements and logic behind VQSR versus hard filtering, please refer to the
official [GATK documentation](https://gatk.broadinstitute.org/hc/en-us/articles/360035531612-Variant-Quality-Score-Recalibration-VQSR-).
```{r filter, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "filter",
targets = NULL, dir = FALSE,
wf_file = "gatk/workflow_gatk_variantFiltration.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
dependency = c("call_variants")
)
```
### Step10: Extract variants per sample
Following the joint genotyping and hard filtering stages, all variant calls
for the cohort are consolidated within a single, multi-sample VCF file. This
step separates this aggregate data, extracting variants for each individual
sample and saving them into distinct files. Importantly, this process is
configured to retain only those variants that successfully passed the
filtering criteria (marked as PASS), ensuring that the final per-sample files
contain only high-confidence calls ready for downstream interpretation.
```{r create_vcf, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "create_vcf",
targets = "hap_caller",
wf_file = "gatk/workflow_gatk_select_variant.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
inputvars = c(SampleName = "_SampleName_"),
dependency = c("hap_caller", "filter")
)
```
This completes the variant calling phase.
## Inspect VCF files
The following section provides example code for inspecting and manipulating
VCF files within R. Please note that this step is distinct
from the automated pipeline and serves as a guide for downstream exploratory
analysis.
### Import VCFs into R
To import VCF data into R, the `readVcf` function is utilized. This loads
the data into `VCF` or `VRanges` class objects, which provide highly
structured and efficient formats for handling variant data, facilitating
tasks such as custom SNP quality filtering and genomic annotation.
```{r inspect_vcf, eval=FALSE}
library(VariantAnnotation)
vcf_raw <- getColumn(sal, "create_vcf")
vcf <- readVcf(vcf_raw[1], "Homo sapiens")
vcf
vr <- as(vcf, "VRanges")
vr
```
### Filter variants in R
The `filter_vars` function facilitates flexible, post-hoc filtering of VCF files
based on user-defined quality parameters. This utility sequentially imports
each VCF file into the R environment, converts the data into an internal
`VRanges` object, applies the specified filtering logic, and exports the
remaining variants to a new, subsetted VCF file.
Filtering criteria are passed to the function as character strings. These
strings are evaluated dynamically to subset the VRanges object using standard R
indexing syntax (e.g., vr[filter, ]).
Example Scenario: The code below demonstrates filtering for variants that meet the following criteria:
(1) Read Support: minimum total read depth >= x; (2) Allele Fraction: At least 80%
of the reads support the called variant; (3) Upstream Status: the variant must have
passed a specific subset of "soft filters" (e.g. PASS tags) inherited from the
upstream GATK analysis.
Note: While the GATK pipeline (specifically Steps 9 and 10) typically performs
the primary hard filtering, this R-based step is provided as an optional tool
for exploratory data analysis or secondary fine-tuning without the need to
re-execute the command-line workflow.
```{r filter_vcf, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
source("VARseq_helper.R")
vcf_raw <- getColumn(sal, "create_vcf")
library(VariantAnnotation)
filter <- "totalDepth(vr) >= 20 & (altDepth(vr) / totalDepth(vr) >= 0.8)"
vcf_filter <- suppressWarnings(filter_vars(vcf_raw, filter, organism = "Homo sapiens", out_dir = "results/vcf_filter"))
},
step_name = "filter_vcf",
dependency = "create_vcf",
run_step = "optional"
)
```
This validation code allows users to compare the VCF files before and after the
filtering process. It serves as a quality control measure to verify the impact
of the applied filters. Please ensure that the optional `filter_vcf` step has
been successfully executed before running this comparison.
```{r check_filter, eval=FALSE}
copyEnvir(sal, "vcf_raw", globalenv())
copyEnvir(sal, "vcf_filter", globalenv())
length(as(readVcf(vcf_raw[1], genome = "Homo sapiens"), "VRanges")[, 1])
length(as(readVcf(vcf_filter[1], genome = "Homo sapiens"), "VRanges")[, 1])
```
### Filter statistics
While optional, this step is included in the default workflow configuration to
provide insight into the quality control process. It facilitates a comparative
analysis between the raw variant calls and the filtered set produced in Step 9.
By generating summary statistics and visualizations, this step allows users to
quantitatively assess the impact of the applied filtering criteria.
```{r summary_filter, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
# read in the cohort VCF file
vcf_all <- suppressWarnings(VariantAnnotation::readVcf("./results/samples_filter.vcf.gz", "Homo sapiens"))
filter_values <- VariantAnnotation::filt(vcf_all)
filter_values[is.na(filter_values)] <- "" # ensure character comparisons work
overall_counts <- table(filter_values) |>
dplyr::as_tibble() |>
dplyr::arrange(dplyr::desc(n))
vcf_all_ft <- VariantAnnotation::geno(vcf_all)$FT
passes_per_sample <- apply(vcf_all_ft, 2, function(x) sum(x == "PASS", na.rm = TRUE))
fails_per_sample <- apply(vcf_all_ft, 2, function(x) sum(x != "PASS", na.rm = TRUE))
sample_filter_summary <- dplyr::tibble(
sample = names(passes_per_sample),
passed = passes_per_sample,
filtered = fails_per_sample
)
write.table(overall_counts, file = "results/summary_filter_overall.tsv", sep = "\t", quote = FALSE, row.names = FALSE)
write.table(sample_filter_summary, file = "results/summary_filter_per_sample.tsv", sep = "\t", quote = FALSE, row.names = FALSE)
library(ggplot2)
source("VARseq_helper.R")
png("results/summary_filter_plot.png", width = 800, height = 600)
p_filter <- plot_summary_filter_plot(sample_filter_summary)
print(p_filter)
dev.off()
# clean up RAM
try(rm(vcf_all, vcf_all_ft, filter_values, passes_per_sample, fails_per_sample), silent = TRUE)
invisible(gc())
},
step_name = "summary_filter",
dependency = "filter"
)
```
Figure 2: Variant Filtering Summary per Sample
## Annotate variants
### Annotate variants with SnpEff
In this required step, the workflow utilizes [SnpEff](https://pcingola.github.io/SnpEff/)
to annotate the filtered variants. SnpEff is a robust tool designed to predict the functional effects
of genetic variants on genes and proteins, providing essential biological
context (e.g. amino acid changes, impact severity) for the identified
mutations.
```{r annotate_vcf, eval=FALSE, spr=TRUE}
appendStep(sal) <- SYSargsList(
step_name = "annotate_vcf",
targets = "create_vcf", dir = TRUE,
wf_file = "gatk/snpeff.cwl",
input_file = "gatk/gatk.yaml",
dir_path = "param/cwl",
inputvars = c(SampleName = "_SampleName_", vcf_raw = "_vcf_raw_"),
dependency = c("create_vcf")
)
```
### Combine annotation results among samples
Upon completion of the annotation process for each VCF file, the subsequent
step involves consolidating the results. All annotated VCFs can be imported and
organized into a single list containing large `VRanges` objects, facilitating
efficient downstream analysis and manipulation within R.
```{r combine_var, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
vcf_anno <- getColumn(sal, "annotate_vcf", position = "outfiles", column = "ann_vcf")
clean_vcf_file <- function(path) {
lines <- readLines(path, warn = FALSE)
if (!length(lines)) return(path)
header_start <- which(grepl("^##fileformat=", lines, perl = TRUE))[1]
if (is.na(header_start) || header_start <= 1) return(path)
writeLines(lines[header_start:length(lines)], path)
path
}
vcf_anno <- vapply(vcf_anno, clean_vcf_file, FUN.VALUE = character(1))
vr_from_vcf <- function(path) {
message("Importing annotated VCF: ", path)
vcf <- suppressWarnings(VariantAnnotation::readVcf(path, "Homo sapiens"))
ft <- VariantAnnotation::geno(vcf)$FT
if (!is.null(ft) && !isTRUE(is.character(ft))) {
ft_vec <- as.character(ft)
ft_clean <- matrix(ft_vec, ncol = 1)
if (!is.null(dim(ft))) {
ft_clean <- matrix(ft_vec, nrow = nrow(ft), dimnames = dimnames(ft))
}
VariantAnnotation::geno(vcf)$FT <- ft_clean
}
suppressWarnings(as(vcf, "VRanges"))
}
vcf_vranges <- lapply(vcf_anno, vr_from_vcf)
},
step_name = "combine_var",
dependency = "annotate_vcf"
)
```
### Summary table of variants
To transform the variant objects into a human-readable summary, key
annotation fields and essential variant metadata are extracted and reorganized
into a long-format table. In this structured output, each row corresponds to a
single variant annotation for a specific sample, facilitating easier inspection
and downstream analysis.
```{r summary_var, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
source("VARseq_helper.R")
summary_var <- extract_var_table(vcf_vranges)
utils::write.table(summary_var, file = "results/variant_summary_long.tsv", sep = "\t", quote = FALSE, row.names = FALSE)
},
step_name = "summary_var",
dependency = "combine_var"
)
```
### Plot raw variants table by consequence
This visualization step generates a summary plot displaying the frequency of
variants categorized by their predicted functional consequence for each sample.
This allows for a rapid assessment of the mutational profile and potential biological
impact across the cohort.
```{r plot_var_consequence, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
library(ggplot2)
effect_counts <- summary_var |>
dplyr::count(sample, effect, name = "n")
png("./results/var_consequence_log.png", width = 1000, height = 600)
p_var_consequence_log <- ggplot(effect_counts, aes(sample, n, fill = effect)) +
geom_col(position = "dodge", alpha = 0.8) +
geom_text(aes(label = scales::comma(n)), position = position_dodge(width = 0.9), vjust = -0.2, size = 3) +
scale_y_continuous(trans = "log10", labels = scales::comma_format()) +
labs(y = "Variant count (log10 scale)", title = "Variant consequences (log scale)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p_var_consequence_log)
dev.off()
},
step_name = "plot_var_consequence",
dependency = "summary_var"
)
```
Figure 3: variant consequence summary (log scale)
### Plot variants stats
This step provides a focused visualization of biologically significant
variants. Leveraging the above summary table, it quantifies and plots the
distribution of nonsynonymous changes and variants classified as HIGH impact
across all samples. This enables researchers to identify samples with a high
burden of potentially functional or pathogenic mutations.
```{r plot_var_stats, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
library(ggplot2)
plot_summary_data <- summary_var |>
dplyr::filter(
consequence %in% c(
"nonsynonymous_variant", "stop_gained", "frameshift_variant",
"splice_acceptor_variant", "splice_donor_variant",
"start_lost", "stop_lost"
),
effect == "HIGH"
) |>
dplyr::mutate(
sample = factor(sample, levels = getColumn(sal, step = "annotate_vcf", position = "targetsWF", column = "SampleName")),
sex = getColumn(sal, step = "annotate_vcf", position = "targetsWF", column = "Factor")[sample]
)
png("./results/var_summary.png")
p_var_summary <- ggplot(plot_summary_data) +
geom_bar(aes(x = sample, fill = sex), alpha = 0.75) +
scale_fill_brewer(palette = "Set2") +
theme_minimal()
print(p_var_summary)
dev.off()
},
step_name = "plot_var_stats",
dependency = "summary_var"
)
```
Figure 4: high impact variant summary
### Compare high-impact variants by group
To investigate potential disparities in the mutational burden between sample
groups (e.g., Male vs. Female), this step utilizes the `plot_summary_data` object
to generate a comparative boxplot . The visualization displays the distribution
of high-impact variants for each group and overlays a Wilcoxon rank-sum test
p-value, providing an immediate statistical assessment of the difference in
variant loads.
```{r plot_var_boxplot, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
source("VARseq_helper.R")
# change the sample and group columns as needed
p_summary_boxplot <- plot_summary_boxplot(plot_summary_data, sample_col = "sample", group_col = "sex")
library(ggplot2)
png("./results/var_summary_boxplot.png")
print(p_summary_boxplot)
dev.off()
},
step_name = "plot_var_boxplot",
dependency = "plot_var_stats"
)
```
Figure 5: distribution of high-impact variants per sex.
## Venn diagram of variants
The venn diagram utilities defined by the `systemPipeR` package can be used to
identify common and unique variants reported for different samples
and/or variant callers. The below generates a 3-way venn diagram
comparing 3 samples.
```{r venn_diagram, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
unique_samples <- summary_var |> dplyr::distinct(sample) |> dplyr::pull(sample)
if (!length(unique_samples)) {
stop("No samples available in `summary_var`; cannot draw Venn diagram.")
}
top_n <- min(3, length(unique_samples))
selected_samples <- unique_samples[seq_len(top_n)]
variant_df <- summary_var |>
dplyr::filter(sample %in% selected_samples) |>
dplyr::distinct(sample, seqnames, start, ref, alt) |>
dplyr::mutate(variant_id = paste0(seqnames, ":", start, "_", ref, "/", alt))
variant_sets <- split(variant_df$variant_id, variant_df$sample)
vennset <- overLapper(variant_sets, type = "vennsets")
png("./results/vennplot_var.png")
vennPlot(vennset, mymain = "Venn Plot of First 3 Samples", mysub = "", colmode = 2, ccol = c("red", "blue"))
dev.off()
},
step_name = "venn_diagram",
dependency = "summary_var"
)
```
Figure 6: Venn Diagram for 3 samples
## Plot variants programmatically
This section demonstrates how to visualize specific genomic variants in a
genome-browser style layout using the `ggbio` package. The resulting example
plot organizes the data into three synchronized panes: the variant location, the
corresponding read alignments, and the coverage depth. The input requirements for
this step are sorted and indexed BAM files generated in GATK Step 5, as well as
variant information from the `create_vcf` step.
```{r plot_variant, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
source("VARseq_helper.R")
first_high <- summary_var |>
dplyr::filter(effect == "HIGH") |>
head(n = 1)
library(ggbio)
library(VariantAnnotation)
mychr <- as.character(first_high$seqnames)
mystart <- as.numeric(first_high$start) - 500
myend <- as.numeric(first_high$end) + 500
bam_path <- getColumn(sal, "fix_tag")[first_high$sample]
vcf_path <- getColumn(sal, step = "create_vcf")[first_high$sample]
vcf <- suppressWarnings(readVcf(vcf_path, "Homo sapiens"))
ga <- readGAlignments(bam_path, use.names = TRUE, param = ScanBamParam(which = GRanges(mychr, IRanges(mystart, myend))))
vcf_chr <- normalize_ft(simplify_info(vcf[seqnames(vcf) == mychr]))
vr <- suppressWarnings(as(vcf_chr, "VRanges"))
vr_region <- vr[start(vr) >= mystart & end(vr) <= myend]
if (!length(vr_region)) {
vr_region <- vr
}
p1 <- autoplot(ga, geom = "rect")
p2 <- autoplot(ga, geom = "line", stat = "coverage")
p3 <- autoplot(vr_region, type = "fixed") +
xlim(mystart, myend) +
theme(legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank())
p1_3 <- tracks(
place_holder = ggplot2::ggplot(),
Reads = p1,
Coverage = p2,
Variant = p3,
heights = c(0, 0.3, 0.2, 0.1)
) + ylab("")
ggbio::ggsave(p1_3, filename = "./results/plot_variant.png", units = "in")
},
step_name = "plot_variant",
dependency = "summary_var"
)
```
Figure 7: Variant plot for selected region
## Non-Synonymous Variants
This step filters the `variant_summary_long.tsv` table to isolate genes affected
specifically by non-synonymous variants. By focusing on these protein-altering
changes, the workflow identifies functionally relevant candidates. The
corresponding Entrez IDs for these genes are extracted and written to a file
(`common_nonsyn_entrez`), serving as the input for downstream pathway analysis.
```{r non_syn_vars, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
vardf <- read.delim("results/variant_summary_long.tsv")
source("VARseq_helper.R")
common_nonsyn_entrez <- filterNonSyn(df=vardf)
writeLines(common_nonsyn_entrez, "results/common_nonsyn_entrez")
},
step_name = "non_syn_vars",
dependency = "plot_var_consequence"
)
```
## Pathway Enrichment
Reactome pathway enrichment analysis is performed on the gene set identified in
the previous step (genes harboring non-synonymous changes). This step utilizes
the `fora` function (overrepresentation analysis) from the `fgsea` package to
identify statistically significant biological pathways. The results, including
enrichment scores and p-values, are exported to a tabular file (`foraRes.xls`).
```{r pathenrich, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
common_nonsyn_entrez <- readLines("results/common_nonsyn_entrez")
source("VARseq_helper.R")
reacdb <- load_reacList(org="R-HSA")
library(fgsea); library(data.table); library(ggplot2)
foraRes <- fora(genes=common_nonsyn_entrez, universe=unique(unlist(reacdb)), pathways=reacdb)
if (!dir.exists("results/fea")) dir.create("results/fea", recursive = TRUE)
foraRes$overlapGenes <- vapply(foraRes$overlapGenes, toString, FUN.VALUE = character(1))
write.table(foraRes, file = "results/fea/foraRes.xls", row.names = FALSE, sep = "\t", quote = FALSE)
foraRes$pathway <- gsub("\\(.*\\) ", "", foraRes$pathway)
foraRes$pathway <- factor(foraRes$pathway, levels = rev(foraRes$pathway))
png("./results/fea/pathenrich.png", width = 680)
ggplot(head(foraRes, 15), aes(pathway, overlap, fill = padj)) +
geom_bar(position="dodge", stat="identity") + coord_flip() +
scale_fill_distiller(palette = "RdBu", direction=-1, limits = range(head(foraRes$padj, 15))) +
theme(axis.text=element_text(angle=0, hjust=1, size=12), axis.title = element_text(size = 14))
dev.off()
},
step_name = "pathenrich",
dependency = "non_syn_vars"
)
```
Figure 8: Enrichment result of top 15 pathways.
## Drug-Target Analysis
This step utilizes the `drugTargetInteractions` package to identify drugs that
target the proteins affected by non-synonymous changes (identified above).
While the tool supports high-throughput batch analysis of large protein sets,
this specific analysis focuses on a biologically relevant subset: proteins
mapping to the "Diseases of metabolism" pathway. This pathway was selected
based on its high enrichment ranking in the previous Reactome analysis and the
high frequency with which its proteins serve as established therapeutic targets
for genetic diseases.
```{r drug_target_analysis, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
## Configure paths for drugTargetInteractions. Under chembldb provide path to chembl_xx.db on your system
# chembldb <- system.file("extdata", "chembl_sample.db", package="drugTargetInteractions")
resultsPath <- "results/drug_target/"
config <- drugTargetInteractions::genConfig(chemblDbPath=chembldb, resultsPath=resultsPath)
downloadUniChem(config=config)
cmpIdMapping(config=config)
foraRes <- read.delim("results/fea/foraRes.xls")
entrez_ids <- unlist(strsplit(foraRes[13,7], ", ")) # select here pathway of interest
source("VARseq_helper.R")
drugMap <- runGeneTargetDrug(entrez=entrez_ids)[[1]]
drugMap <- drugMap[!grepl("Query_", drugMap$GeneName), c("GeneName", "UniProt_ID", "Target_Desc", "Drug_Name", "CHEMBL_CMP_ID", "MOA", "Mesh_Indication")]
drugMap <- drugMap[!is.na(drugMap$CHEMBL_CMP_ID),]
if (!dir.exists("results/drug_target")) dir.create("results/drug_target", recursive = TRUE)
write.table(drugMap, file="results/drug_target/drug_target.xls", row.names=FALSE, sep="\t", quote=FALSE)
},
step_name = "drug_target",
dependency = "pathenrich"
)
```
The resulting table below lists the identifiers, names, and mechanisms of
action for drugs known to target these specific metabolic proteins.
```{r read_drug_table, eval=TRUE}
drugMap <- read.delim("results/drug_target.xls")
DT::datatable(drugMap)
```
## Version information
```{r sessionInfo, eval=FALSE, spr=TRUE}
appendStep(sal) <- LineWise(
code = {
sessionInfo()
},
step_name = "sessionInfo",
dependency = "drug_target")
```
# Additional details
## Running workflows
The `runWF` function is the primary tool for executing workflows. It runs the
code of the workflow steps after loading them into a `SAL` workflow container.
The workflow steps can be loaded interactively one by one or in batch mode with
the `importWF` function. The batch mode is more convenient and is the intended
method for loading workflows. It is part of the standard routine for running
workflows introduced in the [Quick start](#quick-start) section.
```{r runWF, eval=FALSE}
sal <- runWF(sal)
```
## Parallelization on clusters
The processing time of computationally expensive steps can be greatly
accelerated by processing many input files in parallel using several CPUs
and/or computer nodes of an HPC or cloud system, where a scheduling system is
used for load balancing. To simplify for users the configuration and execution
of workflow steps in serial or parallel mode, `systemPipeR` users for both the
same `runWF` function. Parallelization simply requires appending of the
parallelization parameters to the settings of the corresponding workflow steps
each requesting the computing resources specified by the user, such as the
number of CPU cores, RAM and run time. These resource settings are stored in
the corresponding workflow step of the `SAL` workflow container. After adding the
parallelization parameters, `runWF` will execute the chosen steps in parallel
mode as instructed.
The following example applies to an alignment step of a VAR-Seq workflow.
In the chosen alignment example, the parallelization
parameters are added to the alignment step (here `bwa_alignment`) of `SAL` via
a `resources` list. The given parameter settings will run 8 processes (`Njobs`) in
parallel using for each 4 CPU cores (`ncpus`), thus utilizing a total of 32 CPU
cores. The `runWF` function can be used with most queueing systems as it is based on
utilities defined by the `batchtools` package, which supports the use of
template files (_`*.tmpl`_) for defining the run parameters of different
schedulers. In the given example below, a `conffile` (see
_`.batchtools.conf.R`_ samples [here](https://mllg.github.io/batchtools/)) and
a `template` file (see _`*.tmpl`_ samples
[here](https://github.com/mllg/batchtools/tree/master/inst/templates)) need to be present
on the highest level of a user's workflow project. The following example uses the sample
`conffile` and `template` files for the Slurm scheduler that are both provided by this
package.
The `resources` list can be added to analysis steps when a workflow is loaded into `SAL`.
Alternatively, one can add the resource settings with the `addResources` function
to any step of a pre-populated `SAL` container afterwards. For workflow steps with the same resource
requirements, one can add them to several steps at once with a single call to `addResources` by
specifying multiple step names under the `step` argument.
```{r runWF_cluster, eval=FALSE}
# wall time in mins, memory in MB
resources <- list(conffile=".batchtools.conf.R",
template="batchtools.slurm.tmpl",
Njobs=8,
walltime=120,
ntasks=1,
ncpus=4,
memory=1024,
partition = "short"
)
sal <- addResources(sal, c("bwa_alignment"), resources = resources)
sal <- runWF(sal)
```
## Visualize workflow
_`systemPipeR`_ workflows instances can be visualized with the `plotWF` function.
```{r plotWF, eval=FALSE}
plotWF(sal, rstudio = TRUE)
```
## Checking workflow status
The following commands can be used to check the workflow status. The first one gives a summary
and the second one a more detailed outline.
```{r statusWF, eval=FALSE}
sal
statusWF(sal)
```
## Accessing logs report
`systemPipeR` consolidates all workflow execution logs into a central location
within the `.SPRproject` directory. This organization simplifies the inspection
of standard output (stdout) and error (stderr) streams generated by both
command-line and R steps. Additionally, this log data serves as the
foundation for generating `systemPipeR's` technical reports with the `renderLogs`
function.
```{r logsWF, eval=FALSE}
sal <- renderLogs(sal)
```
## Tools used
To inspect the command-line tools utilized in this workflow, the
`listCmdTools` function can be used. Additionally, when using a module
environment is used one can verify the loaded modules with `listCmdModules`.
The following code block will generate a report listing the tools required for
running this workflow. If the workflow is used for the first time or the user wishes
to browse the workflow's requirements, the code below will display the default
toolset configuration.
```{r list_tools, eval=TRUE}
if(file.exists(file.path(".SPRproject", "SYSargsList.yml"))) {
local({
sal <- systemPipeR::SPRproject(resume = TRUE)
systemPipeR::listCmdTools(sal)
systemPipeR::listCmdModules(sal)
})
} else {
cat(crayon::blue$bold("Tools and modules required by this workflow are:\n"))
cat(c("trimmomatic/0.39", "samtools/1.14", "gatk/4.2.0.0", "bcftools/1.15",
"bwa/0.7.17", "snpEff/5.3"), sep = "\n")
}
```
## Report Session Info
This is the session information that will be included when rendering this report.
```{r report_session_info, eval=TRUE}
sessionInfo()
```
# Funding
This project is funded by awards from the National Science Foundation
([ABI-1661152](https://www.nsf.gov/awardsearch/showAward?AWD_ID=1661152)], and
the National Institute on Aging of the National Institutes of Health
([U19AG023122](https://reporter.nih.gov/project-details/9632486)).
# References