Hi all, hope you are all having a nice day. I am w...
# modules
j
Hi all, hope you are all having a nice day. I am working on contributing modules leafcutter/clusterregtools (clusters splice junctions from RNA-seq) and leafcutter/differentialsplicing (differential splicing analysis based on output from leafcutter/clusterregtools). I will need to add some minimal RNA-seq test data, and I wanted to check if my approach is okay. I have used ENCODE4 knockdown data for RBFOX2 in K562, with two KD and two control replicates, and subsetted each BAM to chr22. I obtained the necessary junction files for clusterregtools input by running
regtools junctions extract -a 8 -m 50 -M 500000 -s 0
for each BAM. To get the counts file for differentialsplicing input, I ran leafcutter_cluster_regtools.py on those 4 junction files with params
-m 50 -l 500000
. These are the files that I would like to add as test data: For leafcutter/clusterregtools: • 4 junction files (
.junc
) from the ENCODE4 2 controls vs. 2 RBFOX KD samples in K562, amounting to 1000 junctions each (~80KB per file) For leafcutter/differentialsplicing: • 1 groups file (
.txt
) mapping samples to conditions (<1KB) • 1 counts file (
_perind_numers.c_ounts.gz
) from LeafCutter clustering on the above 4 junction files (<1KB) I have confirmed that leafcutter runs successfully with this minimal test data (tested with
-i 2 -g 2
params for a simple 2v2 comparison). Please let me know if this seems appropriate, in which case I can go ahead and open a PR to add this test data, or whether you would recommend any adjustments. If any additional info would be helpful just let me know. Thank you!
e
You can use the
regtools_junctionsextract
module as a setup step in the nf-test to produce the required data, except if for some reason it is very slow and you want to bypass that. Also, for very niche text files that would not be used elsewhere, it is good practice to just print them on the go (see
clipkit
module for example), if possible clipkit test example writing the file on the go:
Copy code
test("PF14720 - aln - auxiliary file (CST)") {
        when {
            params {
                module_args = "--mode cst"
            }
            process {
                """
                input[0] = GUNZIP.out.gunzip
                input[1] = []
                input[2] = Channel.of(
                    "1\tkeep", "2\ttrim", "3\tkeep", "4\ttrim", "5\tkeep"
                ).collectFile(name: 'cst_auxiliary.tsv', newLine: true)
                """
            }
        }
j
Thank you very much for the helpful feedback! I can definitely look into simplifying to use a setup step like you mention for leafcutter/clusterregtools to bypass uploading any junction files, since that runs very quickly typically. I'll omit the group file as well, good to know. Would it still be OK for me to open a PR for the counts file though? (edit: actually, I guess I could see about writing a mock file on the go for that potentially as well)
Although, I'm realizing that even if I changed things to use regtools_junctionsextract to get the junction files in a setup step, I still didn't see any suitable RNA-seq BAMs within the existing test data that I could run on, unless you had any suggestions