bguo068 / tskibd

Calculate IBD segments from tree sequence using tskit C API
MIT License
1 stars 0 forks source link

Suggestion for command line argument #3

Closed sdtemple closed 1 year ago

sdtemple commented 1 year ago

Here is a small suggestion as I ran into a bug in implementing your package.

I got an error in a snakemake pipeline with the following code block. The problem is that your 1.ibd file gets written to the current working directory, even if the trees file is elsewhere.

# execute the main program
rule tskibd_run:
    input:
        trees='{macro}/{micro}/{seed}/recap.trees',
    output:
        ibd='{macro}/{micro}/{seed}/1.ibd',
    params:
        datfolder='{macro}/{micro}/{seed}',
    shell:
        """
        build/tskibd \
            1 \
            1000000 \
            10000 \
            3 \
           {input.trees}     
        """

An implementation in my snakemake pipeline that worked is the following. I had to move to a new working directory. It would be nice if you could specify the output file prefix specifically in another command line argument.

# execute the main program
rule tskibd_run:
    input:
        trees='{macro}/{micro}/{seed}/recap.trees',
    output:
        ibd='{macro}/{micro}/{seed}/1.ibd',
    params:
        datfolder='{macro}/{micro}/{seed}',
    shell:
        """
        cd {params.datfolder}
        build/tskibd \
            1 \
            1000000 \
            10000 \
            3 \
            recap.trees     
        """
bguo068 commented 1 year ago

Thanks for catching this issue. This has been fixed in the commit 95511c4.

Now you should be able to either use the default output prefix which is the chromosome number, or the specified output prefix:

# 1. use chromosome number as output prefix
build/tskibd 1 15000 150 2  example_data/chr1.trees
# 2. use specified output prefix
build/tskibd 1 15000 150 2  example_data/chr1.trees chr1
# 3. use specified output prefix under a given folder (the folder will be created if it does not exist)
build/tskibd 1 15000 150 2  example_data/chr1.trees out/put/dir/chr1