ChipSeq/chipseq.sh at master · BradyLab/ChipSeq · GitHub

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
data_path= #change to some path where your data is stored
           ##need full data path

WT= #header name of WT files
KO= #header name of KO files
index= # name of index file

##################################################
############## Convert sra to fastq ##############
##################################################

fastq-dump ${data_path}/${WT}.sra
fastq-dump ${data_path}/${KO}.sra

WT_path=${data_path}/${WT}
KO_path=${data_path}/${KO}


#################################################
############### Quality Control #################
#################################################

#Before
#TODO find out why can only do from this dir
./FastQC/fastqc ${WT_path}.fastq
./FastQC/fastqc ${KO_path}.fastq

#Trim and QC
trimReads -q 20 -m 20 -f ${data_path}/adapters.fasta ${WT_path}.fastq
trimReads -q 20 -m 20 -f ${data_path}/adapters.fasta ${KO_path}.fastq

#After
./FastQC/fastqc ${WT_path}.trimmed.fastq
./FastQC/fastqc ${KO_path}.trimmed.fastq


#################################################
################### BWA #########################
#################################################

python aln_bwa.py -n ${WT_path}.trimmed -i ${data_path}/${index}
python aln_bwa.py -n ${KO_path}.trimmed -i ${data_path}/${index}


#TODO add code for removing dups
#################################################
################### SICER #######################
#################################################

mkdir ${data_path}/SICER_res/
python determine_gapsize.py --ko ${data_path}/sorted-${KO}.sam.bed --wt ${data_path}/sorted-${WT}.sam.bed -w 200 -o ${data_path}/SICER_res

#run graph_gap_size.R
#TODO make a wrapper for this

#run with best parmas
best_gap=800
sh ~/Downloads/SICER_V1.1/SICER/SICER-df-rb.sh ${data_path}/sorted-${KO}.sam.bed ${data_path}/sorted-${WT}_head.sam.bed 200 ${best_gap} 0.1 0.001 > $data_path/SICER_res/


#################################################
############## Gene enrichment ##################
#################################################

python chipseq_peak.py \
  -e ${data_path}/sorted-${KO}.sam-W200-G${best_gap}-E0.001.scoreisland
  -c ${data_path}/sorted-${WT}_head.sam-W200-G${best_gap}-E0.001.scoreisland
  --pad 1000
  --genes ## file with start stop postion for genes/ bed file
  -o ## name of output file