Simulations guide ChIPsim

ChIPsim Guide

• download the ChIPsim code

• All the chromosome files must be in the same directory, FASTA format and have the chr?.fa extension, in order to be accepted. The chr? part will be used to create the chr?.fastq and control_chr?.fastq files in the same directory.

• It is required to specify the path to the directory with the chromosomes, and a number for the seed while it is optional to specify a number for the cores that will be used (default 1).

.\chipsim.py --dir /path-to-your-dir/ -s INT -n INT

usage: ChIPsim [-h] -d DIR -s S [-n N]

Simulated ChIPseq data for treatment and control using R package "ChIPsim"

optional arguments:
  -h, --help         show this help message and exit
  -n N               number of cores

required arguments:
  -d DIR, --dir DIR  directory with chromosome .fa files
  -s S               seed number for the simulation

• The python code does not perform any additional check if the files provided are FASTA format, if there are nucleotides other than ATGCN etc, all these are handled by the ChIPsim package. Also we should note that although the user selects only one seed numbed, in reality (#chromosomes * 5 scripts) = # seed will be generated, and this is the reason that a dictionary assigning a number to each chromosome file is created before calling the ChIPsim package.

#!/usr/bin/env python3

import sys
import os
import subprocess
from multiprocessing import Pool
import argparse
from pathlib import Path

parser = argparse.ArgumentParser(description="Simulated ChIPseq data for treatment and control using R package \"ChIPsim\"", prog="ChIPsim")

parser.add_argument("-d", "--dir", help = "directory with chromosome .fa files", type = Path, required = True)
parser.add_argument("-s", help = "seed number for the simulation", type = int, default = 1)
parser.add_argument("-n", help ="number of cores", type = int, default = 1)

args = parser.parse_args()

try:
    directory=args.dir
    print(directory)
    os.chdir(directory)
    fafiles = []
    for file in os.listdir(directory):
        if file.endswith(".fa"):
            fafiles.append(file)
    num_chromosomes = len(fafiles)
    if num_chromosomes !=0:
        print(f"There are {num_chromosomes} chromosomes.")
    else:
        print("No chromosome files in this directory. Files must have .fa  extension.")
        sys.exit()
    try:
        path = f'{directory}/FeatFiles'
        os.mkdir(path)
    except:
        print("FeatFiles directory exists.")

except:
    print('Please pass directory.')
    sys.exit()

num_processes = args.n
seed = args.s

chr_dir = {}
for chrom in enumerate(fafiles,1):
    chr_dir[chrom[1]] = chrom[0]

def calculate(func, args):
    result = func(*args)
    return result

def chipsim1(fafile, chr, seed):
    seed = chr_dir[f'{fafile}'] * seed
    return subprocess.run(['Rscript', f'{directory}/rscripts/chipsim1.r', f'{fafile}', f'{chr}', f'{seed}'])

def chipsim2(fafile, chr, seed):
    seed = chr_dir[f'{fafile}'] * seed
    return subprocess.run(['Rscript', f'{directory}/rscripts/chipsim2.r', f'{fafile}', f'{chr}', f'{seed}'])

def chipsim3(fafile, chr, seed):
    seed = chr_dir[f'{fafile}'] * seed
    return subprocess.run(['Rscript', f'{directory}/rscripts/chipsim3.r', f'{fafile}', f'{chr}', f'{seed}'])

def chipsim4(fafile, chr, seed):
    seed = chr_dir[f'{fafile}'] * seed
    return subprocess.run(['Rscript', f'{directory}/rscripts/chipsim4.r', f'{fafile}', f'{chr}', f'{seed}'])
    
def chipsim5(fafile, chr, seed):
    seed = chr_dir[f'{fafile}'] * seed
    return subprocess.run(['Rscript', f'{directory}/rscripts/chipsim5.r', f'{fafile}', f'{chr}', f'{seed}'])


def test():
    with Pool(processes=num_processes) as pool:
        TASKS = [(chipsim1, (file, file.split(".")[0], seed + 1)) for file in fafiles]+ \
        [(chipsim2, (file, file.split(".")[0], seed + 2)) for file in fafiles] + \
        [(chipsim3, (file, file.split(".")[0], seed + 3)) for file in fafiles] + \
        [(chipsim4, (file, file.split(".")[0], seed + 4)) for file in fafiles] + \
        [(chipsim5, (file, file.split(".")[0], seed + 5)) for file in fafiles]
        results = [pool.apply_async(calculate, t) for t in TASKS]
        for r in results:
            r.get()
    


if __name__ == '__main__':
    test()