-
Beibei Chen authored5c34de27
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
Forked from
BICF / Astrocyte / chipseq_analysis
785 commits behind the upstream repository.
process.py 3.23 KiB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#!/usr/bin/python
# programmer : bbc
# usage: main function to call all the procedures for chip-seq analysis
import sys
import os
import argparse as ap
import logging
import pandas as pd
import glob
import subprocess
from multiprocessing import Pool
import runDeepTools
import runMemechip
logging.basicConfig(level=10)
def prepare_argparser():
description = "Make wig file for given bed using bam"
epilog = "For command line options of each command, type %(prog)% COMMAND -h"
argparser = ap.ArgumentParser(description=description, epilog = epilog)
argparser.add_argument("-i","--input",dest = "infile",type=str,required=True, help="input design file")
argparser.add_argument("-g","--genome",dest = "genome",type=str,required=True, help="genome", default="hg19")
argparser.add_argument("--top-peak",dest="toppeak",type=int, default=-1, help = "Only use top peaks for motif call")
#argparser.add_argument("-s","--strandtype",dest="stranded",type=str,default="none", choices=["none","reverse","yes"])
#argparser.add_argument("-n","--name",dest="trackName",type=str,default="UserTrack",help = "track name for bedgraph header")
return(argparser)
def memechip_wrapper(args):
#print args
runMemechip.run(*args)
def main():
argparser = prepare_argparser()
args = argparser.parse_args()
#dfile = pd.read_csv(args.infile)
#for testing, add testing path to all input files
test_path = "/project/BICF/BICF_Core/bchen4/chipseq_analysis/test/"
designfile = pd.read_csv(args.infile)
designfile['Peaks'] = designfile['Peaks'].apply(lambda x: test_path+x)
designfile['bamReads'] = designfile['bamReads'].apply(lambda x: test_path+x)
designfile['bamControl'] = designfile['bamControl'].apply(lambda x: test_path+x)
designfile.to_csv(args.infile+"_new",index=False)
dfile = pd.read_csv(args.infile+"_new")
#call deeptools
runDeepTools.run(args.infile+"_new", args.genome)
#call diffbind
this_script = os.path.abspath(__file__).split("/")
folder = "/".join(this_script[0:len(this_script)-1])
diffbind_command = "Rscript "+folder+"/runDiffBind.R "+args.infile+"_new"
#logging.debug(diffbind_command)
p = subprocess.Popen(diffbind_command, shell=True)
p.communicate()
#call chipseeker on original peaks and overlapping peaks
chipseeker_command = "Rscript "+folder+"/runChipseeker.R "+",".join(dfile['Peaks'].tolist())+" "+",".join(dfile['SampleID'])
#BC## logging.debug(chipseeker_command)
p = subprocess.Popen(chipseeker_command, shell=True)
p.communicate()
overlapping_peaks = glob.glob('*diffbind.xls')
overlapping_peak_names = []
for pn in overlapping_peaks:
overlapping_peak_names.append(pn.split("_diffbind")[0].replace("!","non"))
chipseeker_overlap_command = "Rscript "+folder+"/runChipseeker.R "+",".join(overlapping_peaks)+" "+",".join(overlapping_peak_names)
p = subprocess.Popen(chipseeker_overlap_command, shell=True)
p.communicate()
#MEME-chip on all peaks
meme_arglist = zip(dfile['Peaks'].tolist(),[test_path+"hg19.2bit"]*dfile.shape[0],[str(args.toppeak)]*dfile.shape[0],dfile['SampleID'].tolist())
#BC# #print meme_arglist
work_pool = Pool(min(12,dfile.shape[0]))
resultList = work_pool.map(memechip_wrapper, meme_arglist)
work_pool.close()
work_pool.join()
if __name__=="__main__":
main()