### Optimize I/O Performance by Copying to Local Scratch Source: https://context7.com/pacificbiosciences/sawfish/llms.txt This example illustrates a common performance optimization strategy: copying large input files (like BAM files) from network storage to a faster local scratch directory before processing. This reduces latency associated with network I/O, potentially speeding up analysis steps that are I/O-bound. ```bash # Optimize for network storage # Copy to local scratch for better I//O cp /network/storage/sample.bam /scratch/sample.bam sawfish discover --bam /scratch/sample.bam ... ``` -------------------------------- ### Generate PDF Methods Document with Docker Pandoc Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/methods/README.md This command generates a PDF methods document using a Dockerized pandoc image, which is useful if pandoc and LaTeX are not installed locally. It mounts the current directory to the container and executes the pandoc command. ```shell docker run --rm -v "$(pwd)":/data -w /data -u $(id -u):$(id -g) pandoc/latex methods.md -f markdown -t pdf -o methods.pdf ``` -------------------------------- ### Generate PDF Methods Document with Pandoc Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/methods/README.md This command generates a PDF methods document from markdown files using pandoc and LaTeX. Ensure pandoc and a LaTeX distribution are installed locally. The command takes the input markdown file and outputs a PDF file. ```shell pandoc methods.md -f markdown -t pdf -o methods.pdf ``` -------------------------------- ### Setup Conda Environment for Excluded Region Tools Source: https://github.com/pacificbiosciences/sawfish/blob/main/scripts/cnv_excluded_regions/README.md Sets up the necessary Conda environment for all excluded region generation scripts. This ensures all required tools are available for processing. ```bash conda env create --file excluded_region_tools.yaml ``` -------------------------------- ### CNV Genotype Formatting Examples Source: https://context7.com/pacificbiosciences/sawfish/llms.txt Demonstrates CNV genotype conventions for different copy number states (0, 1, gain, 4) and haploid regions. It highlights the relationship between GT (Genotype), CN (Copy Number), and CNQ (Copy Number Quality). ```bash # CNV-only genotype conventions (not merged with breakpoint SV) # Copy number 0 (homozygous deletion) in diploid region # GT=0/0 with CN=0 # chr1 5000000 sawfish:0:234:0:0 N 567 PASS SVTYPE=DEL;END=5100000;SVLEN=100000;SVCLAIM=D GT:CNQ:CN 0/0:567:0 # Copy number 1 (heterozygous deletion) in diploid region # GT=0/1 with CN=1 # chr1 6000000 sawfish:0:345:0:0 N 432 PASS SVTYPE=DEL;END=6050000;SVLEN=50000;SVCLAIM=D GT:CNQ:CN 0/1:432:1 # Copy number gain (CN≥3) in diploid region # GT=./1 (one allele duplicated, other unknown) # chr1 7000000 sawfish:0:456:0:0 N 389 PASS SVTYPE=DUP;END=7200000;SVLEN=200000;SVCLAIM=D GT:CNQ:CN ./1:389:3 # Copy number 4 (both alleles may be duplicated or more complex) # GT=./1 (no allele-specific CN estimation) # chr1 8000000 sawfish:0:567:0:0 N 456 PASS SVTYPE=DUP;END=8300000;SVLEN=300000;SVCLAIM=D GT:CNQ:CN ./1:456:4 # Haploid region (e.g., chrY) with CN=0 (expected for XX sample) # chr1 9000000 sawfish:0:678:0:0 N 234 PASS SVTYPE=DEL;END=9100000;SVLEN=100000;SVCLAIM=D GT:CNQ:CN 0:234:0 # Haploid region with CN=1 (expected for XY sample) # GT=0 (reference state) # chr1 10000000 sawfish:0:789:0:0 N 123 MinQUAL SVTYPE=DEL;END=10050000;SVLEN=50000;SVCLAIM=D GT:CNQ:CN 0:123:1 ``` -------------------------------- ### Extract GC-Corrected Depth using bigWigToBedGraph Source: https://context7.com/pacificbiosciences/sawfish/llms.txt This command extracts GC-corrected depth information from a BigWig file. The first example shows how to pipe the output to 'head' to view the initial lines, useful for quick inspection. The second example demonstrates extracting depth for a specific genomic region defined by chromosome, start, and end coordinates. ```bash # View GC-corrected depth in IGV or extract specific region bigWigToBedGraph samples/sample0_HG002/gc_bias_corrected_depth.bw \ /dev/stdout | head # Extract depth for specific region bigWigToBedGraph samples/sample0_HG002/gc_bias_corrected_depth.bw \ gc_depth_region.bedgraph -chrom=chr1 -start=1000000 -end=2000000 ``` -------------------------------- ### Run Trio Joint-Call Step with Sawfish (Multiple --sample) Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/user_guide.md Executes the 'joint-call' step for a trio after all 'discover' steps are complete. This method uses multiple '--sample' arguments to specify the output directories from the discover steps. It reuses the reference and BAM paths provided during the discover phase. ```bash sawfish joint-call \ --threads 16 \ --sample HG004_discover_dir \ --sample HG003_discover_dir \ --sample HG002_discover_dir \ --output-dir HG002_trio_joint_call_dir ``` -------------------------------- ### Example Expected Copy Number BED File (hg38, XY) Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/user_guide.md This example demonstrates the format for an expected copy number BED file. It specifies regions on chromosomes X and Y with their corresponding expected copy numbers, useful for analyzing sex chromosomes. The first three columns define the genomic regions, column 5 provides the expected copy number, and column 4 can be used for labels. ```plaintext chrX 0 2781479 chrX_PAR_1 2 chrX 2781479 155701382 chrX_uniq_1 1 chrX 155701382 156040895 chrX_PAR_2 2 chrY 0 2781479 chrY_PAR_1 0 chrY 2781479 56887902 chrY_uniq_1 1 chrY 56887902 57227415 chrY_PAR_2 0 ``` -------------------------------- ### Run Trio Discover Step with Sawfish Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/user_guide.md Executes the 'discover' step for each of the three samples in a trio. This step analyzes individual BAM files, requiring specified threads, reference genome, and output directory. The commands are independent and can be run in parallel. ```bash sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam HG004.GRCh38.bam \ --expected-cn ${DISTRO_ROOT_DIR}/data/expected_cn/expected_cn.hg38.XX.bed \ --cnv-excluded-regions ${DISTRO_ROOT_DIR}/data/cnv_excluded_regions/annotation_and_common_cnv.hg38.bed.gz \ --output-dir HG004_discover_dir sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam HG003.GRCh38.bam \ --expected-cn ${DISTRO_ROOT_DIR}/data/expected_cn/expected_cn.hg38.XY.bed \ --cnv-excluded-regions ${DISTRO_ROOT_DIR}/data/cnv_excluded_regions/annotation_and_common_cnv.hg38.bed.gz \ --output-dir HG003_discover_dir sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam HG002.GRCh38.bam \ --expected-cn ${DISTRO_ROOT_DIR}/data/expected_cn/expected_cn.hg38.XY.bed \ --cnv-excluded-regions ${DISTRO_ROOT_DIR}/data/cnv_excluded_regions/annotation_and_common_cnv.hg38.bed.gz \ --output-dir HG002_discover_dir ``` -------------------------------- ### Run Sawfish Discover Step (Shell) Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/user_guide.md Executes the 'discover' step of Sawfish on a mapped sample BAM file. This step identifies candidate structural variant regions, assembles local SV haplotypes, builds sequencing coverage tracks, and performs initial copy number segmentation. Optional inputs for expected copy number and excluded regions can improve CNV output quality. ```shell sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam HG002.GRCh38.bam \ --expected-cn ${DISTRO_ROOT_DIR}/data/expected_cn/expected_cn.hg38.XY.bed \ --cnv-excluded-regions ${DISTRO_ROOT_DIR}/data/cnv_excluded_regions/annotation_and_common_cnv.hg38.bed.gz \ --output-dir HG002_discover_dir ``` -------------------------------- ### Restrict T2T Benchmark BED to Autosomes Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/accuracy.md Modify the T2T benchmark BED file to restrict confident regions to autosomes only. This is achieved by removing lines starting with 'chrX' or 'chrY' using grep. ```bash grep -v chr[XY] GRCh38_HG2-T2TQ100-V1.1_stvar.benchmark.bed > GRCh38_HG2-T2TQ100-V1.1_stvar.benchmark.only_autosomes.bed ``` -------------------------------- ### Run Trio Joint-Call Step with Sawfish (--sample-csv) Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/user_guide.md Executes the 'joint-call' step for a trio using a CSV file to specify sample inputs. This method consolidates sample information, including discover directory and BAM paths, into a single file. It requires specifying the reference genome directly. ```bash cat << END > all_samples.csv HG002_discover_dir, HG002.GRCh38.bam HG003_discover_dir, HG003.GRCh38.bam HG004_discover_dir, HG004.GRCh38.bam END sawfish joint-call \ --threads 16 \ --ref GRCh38.fa \ --sample-csv all_samples.csv \ --output-dir HG002_trio_joint_call_dir ``` -------------------------------- ### Call SVs with Sawfish Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/accuracy.md Commands to discover and joint-call structural variants (SVs) using Sawfish. Requires a conda environment with Sawfish installed and a reference genome FASTA file and a BAM file with mapped reads. ```bash conda activate sawfish_accuracy_test_v2 sawfish discover --threads 16 --ref human_GRCh38_no_alt_analysis_set/human_GRCh38_no_alt_analysis_set.fasta --bam m84011_220902_175841_s1.pbmm2-1.13.1.GRCh38.bam sawfish joint-call --threads 16 --sample sawfish_discover_output ``` -------------------------------- ### Parallel Discover Steps for Multiple Samples Source: https://context7.com/pacificbiosciences/sawfish/llms.txt This bash script demonstrates how to run the 'sawfish discover' command in parallel for multiple samples. It defines an array of sample names and then loops through them, launching each 'discover' process in the background (`&`). The `wait` command ensures that the script doesn't exit until all background jobs are completed, effectively utilizing multiple CPU cores for sample processing. ```bash # Parallel discover steps for multiple samples samples=(HG002 HG003 HG004) for sample in "${samples[@]}"; do sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam "${sample}.bam" \ --output-dir "${sample}_discover_dir" & done wait ``` -------------------------------- ### View Copy Number Segments in BEDGraph Format Source: https://context7.com/pacificbiosciences/sawfish/llms.txt Demonstrates how to view copy number segments for a sample by examining the `.bedgraph` file. The output format shows chromosome, start position, end position, and the copy number value for each segment. ```bash # View copy number segments for a sample cat samples/sample0_HG002/copynum.bedgraph # chr1 0 500000 2 # chr1 500000 750000 1 # chr1 750000 1000000 2 # chr1 1000000 1500000 3 ``` -------------------------------- ### Run Sawfish Joint-Call Step (Direct Sample Directories) (Shell) Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/user_guide.md Executes the 'joint-call' step of Sawfish using direct paths to 'discover' output directories. This step merges duplicate SV haplotypes, associates them with samples, evaluates read support, genotypes variants, identifies breakends for CNV boundary points, re-segments depth, and unifies SV/CNV calls, outputting results to a VCF file. It reuses reference and BAM paths stored during the 'discover' step. ```shell sawfish joint-call \ --threads 16 \ --sample HG002_discover_dir \ --output-dir HG002_joint_call_dir ``` -------------------------------- ### View Copy Number Summary JSON Source: https://context7.com/pacificbiosciences/sawfish/llms.txt This command displays the content of a JSON file that summarizes copy number information, likely aggregated across chromosomes. The example JSON shows counts of different copy number states (cn_0, cn_1, cn_2, cn_3) for a specific chromosome. ```bash # View copy number summary JSON cat samples/sample0_HG002/copynum.summary.json ``` -------------------------------- ### Run Sawfish Joint-Call Step (Sample CSV) (Shell) Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/user_guide.md Executes the 'joint-call' step of Sawfish using a CSV file to specify sample inputs (discover directory, BAM path). This method is useful for specifying larger numbers of samples and allows direct specification of the reference genome. The primary output is a VCF file containing genotyped SVs and CNVs. ```shell cat << END > all_samples.csv HG002_discover_dir, HG002.GRCh38.bam END sawfish joint-call \ --threads 16 \ --ref GRCh38.fa \ --sample-csv all_samples.csv \ --output-dir HG002_joint_call_dir ``` -------------------------------- ### Increase Threads for Faster Runtime Source: https://context7.com/pacificbiosciences/sawfish/llms.txt This command shows how increasing the number of threads allocated to 'sawfish discover' can directly improve runtime performance. By assigning more processing cores, the analysis can be completed approximately twice as fast, demonstrating linear scaling with respect to thread count for this operation. ```bash # Linear thread scaling - increase threads for faster runtime sawfish discover --threads 32 --bam sample.bam ... # ~2x faster ``` -------------------------------- ### Calculate Total Bases by Copy Number using awk Source: https://context7.com/pacificbiosciences/sawfish/llms.txt This command calculates the total number of base pairs at each copy number state from a BEDGraph file. It iterates through the file, summing the length of segments (end - start) for each copy number (the fourth column), and then prints the results in the END block. This is useful for understanding the genomic distribution of copy number states. ```bash # Calculate total bases at each copy number awk '{sum[$4]+=$3-$2} END {for (cn in sum) print "CN=" cn "\t" sum[cn] " bp"}' \ samples/sample0_HG002/copynum.bedgraph ``` -------------------------------- ### Sawfish Discover Step: Advanced Parameters Source: https://context7.com/pacificbiosciences/sawfish/llms.txt Configures advanced parameters for the `sawfish discover` step, including minimum indel size, disabling path canonicalization, enabling output directory clobbering, and activating debug logging. It also shows how to specify a custom output directory. ```bash # Minimum indel size for variant reporting sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam sample.bam \ --min-indel-size 50 \ --output-dir discover_dir # Disable path canonicalization for relative paths sawfish discover \ --threads 16 \ --ref ../ref/GRCh38.fa \ --bam ../bam/sample.bam \ --disable-path-canonicalization \ --output-dir discover_dir # Allow output directory clobber (overwrite existing) sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam sample.bam \ --output-dir existing_dir \ --clobber # Enable debug logging sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam sample.bam \ --debug \ --output-dir discover_dir # Custom output directory (default: sawfish_discover_output) sawfish discover \ --threads 16 \ --ref GRCh38.fa \ --bam sample.bam \ --output-dir custom_discover_results ``` -------------------------------- ### Interpret VCF Quality Scores and Filters (Bash) Source: https://context7.com/pacificbiosciences/sawfish/llms.txt Explains quality scores (QUAL, GQ, CNQ) and filter descriptions used in Sawfish's VCF output. QUAL indicates the confidence of an alternate allele or non-expected CN, while GQ and CNQ provide per-sample quality for breakpoints and copy numbers, respectively. Filters like PASS, MinQUAL, MaxScoringDepth, InvBreakpoint, and ConflictingBreakpointGT help in variant assessment. Examples show how to extract specific variant types based on these filters. ```bash # Example VCF records with quality annotations bcftools view genotyped.sv.vcf.gz | grep -v "^##" | head -5 # chr1 1651421 sawfish:0:92:0:0 GTAGCCC...TCTA G 504 PASS SVTYPE=DEL;END=1651474;SVLEN=53;HOMLEN=13;HOMSEQ=TAGCCCCTCTGAA;SVCLAIM=J GT:GQ:PL:AD 0/1:387:537,0,387:9,12 # chr2 5043821 sawfish:1:234:0:0 T 812 PASS SVTYPE=INS;END=5043821;SVLEN=1543;SVCLAIM=J GT:GQ:PL:AD 1/1:812:812,120,0:0,24 # Quality score meanings: # QUAL: Phred-scaled confidence of alternate allele in sample set # - For SV: probability of non-reference genotype # - For CNV: probability CN != expected CN # - For merged SV/CNV: maximum of component QUAL values # GQ: Per-sample genotype quality (breakpoint-based) # CNQ: Per-sample copy number quality (depth-based) # All phred scores capped at 999 # Filter descriptions: # PASS: Variant passes all filters # MinQUAL: QUAL < 10 # MaxScoringDepth: Read depth exceeds 12x GC-corrected haploid depth (max 1000) # InvBreakpoint: Part of inversion record (see EVENT ID) # ConflictingBreakpointGT: Genotypes conflict in multi-breakpoint event # Extract high-quality deletions bcftools view -i 'FILTER="PASS" && INFO/SVTYPE="DEL" && QUAL>=100' \ genotyped.sv.vcf.gz > high_quality_deletions.vcf # Extract CNV-supported variants bcftools view -i 'INFO/SVCLAIM~"D"' genotyped.sv.vcf.gz > cnv_supported.vcf # Extract breakpoint-only SVs bcftools view -i 'INFO/SVCLAIM="J"' genotyped.sv.vcf.gz > breakpoint_only.vcf ``` -------------------------------- ### Extract and Analyze Supporting Reads using jq Source: https://context7.com/pacificbiosciences/sawfish/llms.txt This series of commands demonstrates how to process a gzipped JSON file containing supporting read information using the 'jq' tool. It covers viewing the entire JSON, counting supporting reads for a specific variant, filtering variants based on read count, extracting read names, and retrieving the actual supporting reads from a BAM file. ```bash # Extract supporting reads for all variants zcat supporting_reads.json.gz | jq '.' # Count supporting reads for specific variant zcat supporting_reads.json.gz | \ jq '.["sawfish:0:1041:0:0"]["HG002"] | length' # Extract all variant IDs with >10 supporting reads in HG002 zcat supporting_reads.json.gz | \ jq -r 'to_entries[] | select(.value.HG002 | length > 10) | .key' # Get supporting read names for specific variant zcat supporting_reads.json.gz | \ jq -r '.["sawfish:0:1041:0:0"]["HG002"][]' # Extract reads from BAM using supporting read names zcat supporting_reads.json.gz | \ jq -r '.["sawfish:0:1041:0:0"]["HG002"][]' | \ while read qname; samtools view HG002.bam | grep "^${qname}\\s" done > variant_supporting_reads.sam ``` -------------------------------- ### Run Truvari T2T Assessment Bench with Singularity Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/accuracy.md Executes the `truvari bench` command within a singularity container to assess variant calls against a reference VCF and BED file. It requires the reference FASTA, inclusion BED, base VCF, comparison VCF, and an output directory. Options like `--passonly`, `--pick ac`, and `--dup-to-ins` are used for filtering and variant processing. ```bash # Activate conda environment to use singularity: # conda activate sawfish_accuracy_test_v2 sif=library://ctsa/truvari/truvari:4.2.2 ref=human_GRCh38_no_alt_analysis_set/human_GRCh38_no_alt_analysis_set.fasta bench_vcf=GRCh38_HG2-T2TQ100-V1.1_stvar.svtype.vcf.gz bench_bed=GRCh38_HG2-T2TQ100-V1.1_stvar.benchmark.only_autosomes.bed query_vcf=sawfish_joint-call_output/genotyped.sv.vcf.gz # Truvari bench requires that the output directory does not already exist: output_dir=truvari_t2t_result rm -rf $output_dir # Copy this truvari post-processing script from the sawfish repo to the local directory to simplify usage inside of the singularity container: cp -f ${SAWFISH_REPO_ROOT_DIR}/scripts/truvari_utils/process_truvari_ga4gh_vcfs.py . working_dir=$(pwd) # MAFFT (used in the refine step) can be sensitive to the space available in the default tmp directory, so create one on the regular filesystem: tmp_dir=$(mktemp -d) trap 'rm -rf $tmp_dir' EXIT # Step 1: Run bench singularity exec --bind $working_dir:/data --cleanenv --pwd /data $sif \ truvari bench \ --reference $ref \ --includebed $bench_bed \ --base $bench_vcf \ --comp $query_vcf \ --output $output_dir \ --passonly \ --pick ac \ --dup-to-ins ``` -------------------------------- ### Download and Unpack Reference Genome Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/accuracy.md This code snippet downloads the GRCh38 reference genome in a compressed tarball and then unpacks it. The reference genome is essential for downstream read mapping and SV calling. Ensure sufficient disk space for the unpacked files. ```bash wget https://downloads.pacbcloud.com/public/reference-genomes/human_GRCh38_no_alt_analysis_set.tar.2023-12-04.gz tar -xzf human_GRCh38_no_alt_analysis_set.tar.2023-12-04.gz ``` -------------------------------- ### Run Truvari T2T Assessment Refine with Singularity Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/accuracy.md Executes the `truvari refine` command within a singularity container to refine variant calls. This step uses MAFFT for alignment and requires the reference FASTA, the output BED from the bench step, and specifies threads and alignment tool. It also uses a custom temporary directory for MAFFT. ```bash # Step 2: Run refine singularity exec --bind $working_dir:/data --bind /tmp:$tmp_dir --cleanenv --pwd /data $sif \ truvari refine \ --reference $ref \ --regions $output_dir/candidate.refine.bed \ --recount \ --use-region-coords \ --use-original-vcfs \ --threads 16 \ --align mafft \ $output_dir ``` -------------------------------- ### Create Conda Environment for Sawfish Analysis Source: https://github.com/pacificbiosciences/sawfish/blob/main/docs/accuracy.md This command creates a new conda environment named 'sawfish_accuracy_test_v2' with specified versions of Python, pbmm2, sawfish, and singularity. This ensures reproducible analysis by managing software dependencies. Activating this environment is necessary before running subsequent analysis steps. ```bash conda create -y -n sawfish_accuracy_test_v2 python=3.9 bioconda::pbmm2=1.13.1 bioconda::sawfish=0.12.9 conda-forge::singularity ```