Function call from GenomicFeatures causing DESeq script failure for T2T (hs1) genome

[wong-nw]

Describe the bug
Workflow stops at the differential peak calling step when using the hs1 genome for T2T. The error message is as follows:

processing file: _diff_markdown.Rmd
Error in `open.connection()`:
! cannot open the connection
Backtrace:
  1. GenomicFeatures::makeTxDbFromGFF(...)
  3. BiocIO::import(file, format = format, colnames = colnames, feature.type = GFF_FEATURE_TYPES)
  5. rtracklayer::import(FileForFormat(con, format), ...)
  6. rtracklayer (local) .local(con, format, text, ...)
  7. BiocIO:::connection(m, con, "r")
  8. BiocIO:::connectionForResource(manager, resource(x), open = open)
10. base::open.connection(con, open)

Quitting from lines 263-345 [annotate] (_diff_markdown.Rmd)
Execution halted
[Tue Sep 23 09:54:29 2025]
Error in rule DESeq:

Error seems isolated to roughly line 294-300 in scripts/_diff_markdown.Rmd. The GenomicFeatures::makeTxDbFromGFF function is now attempting to call txdbmaker::makeTxDbFromGFF, possibly by trying to reach out to BioConductor to retrieve the function. Ref: https://bioconductor.org/packages/release/bioc/manuals/GenomicFeatures/man/GenomicFeatures.pdf#Rfn.makeTxDbFromGFF.1

To Reproduce

This error is only produced when running within the container. Manual testing via R environment, command line Rscript calls, and slurm submissions were all successful. Isolated script that was obtained via adapting the Snakemake rule is as follows. Minor modifications to file paths were needed to test in singularity shell based on binded (?) directories (singularity shell -B):

#!/bin/bash

module load R

       if [[ -d "/lscratch/$SLURM_JOB_ID" ]]; then
            TMPDIR="/lscratch/$SLURM_JOB_ID/$dirname"
        else
            TMPDIR="/dev/shm/$dirname"
        fi
        TMPDIR_AUC=$TMPDIR/AUC
        TMPDIR_FRAG=$TMPDIR/FRAG
        TMPDIR_VENN=$TMPDIR/VENN

        # pull conditions
        condition1=`echo sgRACK1_MsCENPA_vs_sgNeg_MsCENPA | awk -F"_vs_" '{print $1}'`
        condition2=`echo sgRACK1_MsCENPA_vs_sgNeg_MsCENPA | awk -F"_vs_" '{print $2}'`

        ## Run AUC method
        mkdir -p ${TMPDIR_AUC}
        mkdir -p ${TMPDIR_AUC}/intermediates_dir
        mkdir -p ${TMPDIR_AUC}/knit_root_dir
        cd $TMPDIR_AUC

        Rscript /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/scripts/_diff_markdown_wrapper.R \
            --rmd /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/scripts/_diff_markdown.Rmd \
            --carlisle_functions /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/scripts/_carlisle_functions.R \
            --countsmatrix /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.countsmatrix.csv \
            --sampleinfo /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.sampleinfo.csv \
            --dupstatus dedup \
            --condition1 ${condition1} \
            --condition2 ${condition2} \
            --fdr_cutoff 0.05 \
            --log2fc_cutoff 0.59 \
            --results /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.AUCbased_diffresults.csv \
            --report /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.AUCbased_diffanalysis.html \
            --elbowlimits /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.AUCbased_diffanalysis_elbowlimits.yaml \
            --spiked SPIKEIN \
            --rawcountsprescaled \
            --scalesfbymean \
            --contrast_data /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.txt \
            --tmpdir $TMPDIR_AUC \
            --species hs1 \
            --gtf /data/CCBR_Pipeliner/db/PipeDB/Indices/hs1/genes.gtf

        # change elbow limits to provided log2fc if limit is set to .na.real
        sed -i "s/low_limit: .na.real/low_limit: -0.59/" /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.AUCbased_diffanalysis_elbowlimits.yaml
        sed -i "s/up_limit: .na.real/up_limit: 0.59/g" /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.AUCbased_diffanalysis_elbowlimits.yaml

        ## Run FRAG method
        mkdir -p ${TMPDIR_FRAG}
        mkdir -p ${TMPDIR_FRAG}/intermediates_dir
        mkdir -p ${TMPDIR_FRAG}/knit_root_dir
        cd $TMPDIR_FRAG

        # Do not use --rawcountsprescaled as these counts are not prescaled!
        Rscript /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/scripts/_diff_markdown_wrapper.R \
            --rmd /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/scripts/_diff_markdown.Rmd \
            --carlisle_functions /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/scripts/_carlisle_functions.R \
            --countsmatrix /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.fragmentscountsmatrix.csv \
            --sampleinfo /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.sampleinfo.csv \
            --dupstatus dedup \
            --condition1 ${condition1} \
            --condition2 ${condition2} \
            --fdr_cutoff 0.05 \
            --log2fc_cutoff 0.59 \
            --results /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.fragmentsbased_diffresults.csv \
            --report /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.fragmentsbased_diffanalysis.html \
            --elbowlimits /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.fragmentsbased_diffanalysis_elbowlimits.yaml \
            --spiked SPIKEIN \
            --scalesfbymean \
            --contrast_data /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.txt \
            --tmpdir $TMPDIR_FRAG \
            --species hs1 \
            --gtf /data/CCBR_Pipeliner/db/PipeDB/Indices/hs1/genes.gtf

        # change elbow limits to provided log2fc if limit is set to .na.real
        sed -i "s/low_limit: .na.real/low_limit: -0.59/" /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.fragmentsbased_diffanalysis_elbowlimits.yaml
        sed -i "s/up_limit: .na.real/up_limit: 0.59/g" /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.fragmentsbased_diffanalysis_elbowlimits.yaml

        ## Run venns
        mkdir -p ${TMPDIR_VENN}
        mkdir -p ${TMPDIR_VENN}/intermediates_dir
        mkdir -p ${TMPDIR_VENN}/knit_root_dir
        cd $TMPDIR_VENN
        Rscript /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/scripts/_plot_results_venn.R \
            --aucresults /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.AUCbased_diffresults.csv \
            --fragmentsresults /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.fragmentsbased_diffresults.csv \
            --pdf /vf/users/CCBR/projects/ccbr1431/T2T_202509/sgRACK1/results/peaks/0.05/contrasts/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup/sgRACK1_MsCENPA_vs_sgNeg_MsCENPA.dedup.gopeaks_broad.venn.pdf \
            --title "${condition1}_vs_${condition2}__dedup__gopeaks_broad"

Expected behavior
Solution may be as simple as including txdbmaker package in the container.

Ideal behavior is for Snakemake and Rmarkdown runs to completion, producing the differential peaks and downstream bed/bedgraph/bigwig files

[kelly-sovacool]

workaround: run makeTxdbFromGff() interactively, save the object as an rds file in a shared dir, add it to the config file, and modify the logic in the Rmd & wrapper script to use it when the genome is hs1.

[kelly-sovacool]

The actual problem is that carlisle was not properly setting the singularity bind paths. The GTF could not be accessed regardless of which genome was set.

Fix: b56af1c