openomics / chrom-seek Goto Github PK

It would be awesome if we had a script to automate the process of building reference files for different/new organisms, something where we could give it a minimal set of inputs (such as a genomic fasta file, annotation, and an optional blacklist, etc.) and it would handle the rest.

At the current moment, the pipeline is using environment modules. We need to build docker images for all the software dependencies of the pipeline.

The new features to be added to the cfChIP assay include:

1. Run a DiffBind R script on all "ChIP" samples to:
   • Identify consensus peaks, produce bed file
   • Get normalized counts per sample
   • Get UMAP coordinates based upon this information
2. Annotate consensus peaks with uropa
3. Run Amos's specialized differential RMD
4. If only one group per sample and all samples have a group, create UMAP plot

This is step completes basic functional pipeline for cfChIP.

SICER2 has information pertaining to genome chromosome lengths hard-coded in the following file:
https://github.com/zanglab/SICER2/blob/master/sicer/lib/GenomeData.py

To add support for a new reference genome, a user must do the following:

I emailed the author of SICER2 asking if it would be possible to implement any of the suggestions here. It is worth also noting, if one were to make these changes, one should also update the dtype set for chr. At the current moment, it is set to U6 or a unicode string of length 6, which may not be long enough in some edge-cases (could cause a string to get truncated and cause collisions).

For the time being, we will turn off sicer2 when the reference genome is rheMac10; however, we should explore other solutions if the author cannot implement this on their side:

• Fork and add functionality described in issue 23
• Fork, update sicer/lib/GenomeData.py with rheMac10 information and build docker image

Currently using this code as a fix for ppqt rule to run the chrom-seek pipeline, but need to find a more straightforward solution.

 if paired_end:
    extensionsPPQT={"sorted": "bam", "Q5DD":"bam"}
else:
    extensionsPPQT= {"sorted": "bam", "tagAlign": "gz", "Q5DD": "tagAlign.gz"}

About

The following issue was discovered while running run_spp.R from ppqt:

The error indicates that the R package, caTools, was missing from the docker image. Initially, I could not reproduce this behavior, and from my side, it appeared that caTools was installed:

$ module purge
$ ml singularity
$ singularity exec -B $PWD /data/OpenOmics/SIFs/ppqt_v0.2.0.sif Rscript -e 'library(caTools); sessionInfo()'

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] caTools_1.18.2

loaded via a namespace (and not attached):
[1] compiler_4.1.2 bitops_1.0-7  

After more debugging, it appeared that I could reproduce the error if I module load R before running singularity:

$ module purge
$ module load R
$ ml singularity
$ singularity exec -B $PWD /data/OpenOmics/SIFs/ppqt_v0.2.0.sif Rscript -e 'library(caTools); sessionInfo()'

Error in library(caTools) : there is no package called ‘caTools’
Execution halted

It appears the following environment variable, R_LIBS_SITE, is exported when R is module loaded into a user's $PATH. If that environment variable is set, that will update what .libPaths R internally uses to find R packages. This can cause issues with how the R installation within the docker image finds/loads packages.

The following can be used to fix the problem.

$ module purge
$ module load R
$ ml singularity
$ R_LIBS_SITE='' singularity exec -B $PWD /data/OpenOmics/SIFs/ppqt_v0.2.0.sif Rscript -e 'library(caTools); sessionInfo()'
...

However, it would be better to run singularity with the --containall, -C option. This would prevent this problem in the future and other issues that can occur from sharing filesystems and env variables with the host.

Completed with 487945f and de0906b commits

Hi,

I'm trying to test the pipeline using the test files but all files on https://github.com/OpenOmics/chrom-seek/tree/main/.tests are empty except the the peakcall.tsv and the contrasts.tsv. Can you have some working test files in the github?

Best

FRIP_plot.txt

Incorporate Subrata's code into the pipeline

Sicer kills the dry-run when there are samples lacking input controls.

DiffBind with blocking currently errors out. This causes both DiffBind block and regular DiffBind to fail within the pipeline. Error seems to be that the csv file being used does not contain the blocking information.

About

Make the following changes to how the --peakcall PEAKCALL file is parsed and validated:

1. Update the names of the following columns:
   a. "ChIP" becomes "Sample"
   b. "Input" becomes "InputControl"
2. InputControl becomes an optional column.
   a. Only "Sample" and "Group" should be required
3. Documentation on the CLI and wiki should be updated to include the changes and made to encompass both ATAC and ChIP

For ATAC,

The following can be considered to probably be conditional when running the pipeline:

1. Between macs2 and Genrich- Genrich should be default, macs2 should be conditional
2. Between edgeR and deseq2 in the Diffbind rules - Deseq2 should be default, edgeR should be conditional.
3. Motif analysis (MEME/Homer? rules) could probably be conditional for ATAC.
4. Between genes, prot, protsec, protTSS - I think only prot needs to run by default. This should be a selection as the number of total jobs is ballooning in certain circumstances.

These are suggestions not demands,

Paul

Add --rerun-trigger to run sub command

This adds a mechanism to control how snakemake decides what rules to re-run. By default, this option is set to the following: ['mtime', 'params', 'input', 'software-env', 'code'] which can cause upstream rules to re-run unexpectedly.

Add something to the docs stating if this occurs, then this should be set to --rerun-trigger mtime. This allows for the user to decide how/why rules are re-run.

For ATAC part,

ppqt can fail causing a lot of downstream issues which can lead to bigwigs not being made. If this not necessary for bigwig formation, this should be removed.

Paul