Newer
Older
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
## extract chromosome from bam file
samtools view -bo test.bam /home/brandl/mnt/chip-seq_study/ChIPSeq_February_2014/alignments_trimmed_nomulti_pooled/H2Az.bam 1
########################################################################################################################
## macs2 playground
macs2 -n --broad --gsize
#macs2 -n --gsize
########################################################################################################################
### deeptools
## https://github.com/fidelram/deepTools/wiki/All-command-line-options#bamCorrelate
#bamCorrelate bins --bamfiles $bamFiles --region 10:1:100000 --plotFile="bam_correlation.png" --numberOfProcessors=4 --corMethod spearman
## see how well bam files correlate using untrimmed data
bamCorrelate bins --bamfiles $bamFiles --plotFile="bam_correlation_untrimmed.png" --numberOfProcessors=4 --corMethod spearman
########################################################################################################################
### Peak calling with SPP
Rscript /home/brandl/bin/phantompeakqualtools/run_spp.R -c="../alignments_untrimmed/H2Az_Rep1_Lane2_Lib4454.0001.bam" -savp -out="test.txt"
########################################################################################################################
### other qc
## http://www.nature.com/nmeth/journal/v11/n1/full/nmeth.2786.html
# Teytelman et al. discovered that highly expressed loci were always enriched in ChIP peaks, regardless of which protein was pulled dow
#########################################################################################################################
#### CHANCE
mcdir $baseDir/qc_chance
## convert example to sam format
bamToBed -i $baseDir/alignments_trimmed/H2Az_Rep1_Lane1_Lib4454_ca.bam > H2Az_Rep1_Lane1_Lib4454_ca.bed
samtools view -h -o H2Az_Rep1_Lane1_Lib4454_ca.sam $baseDir/alignments_trimmed/H2Az_Rep1_Lane1_Lib4454_ca.bam
/local/home/henry/bin/CHANCE/run_chance_linux.sh /usr/local/MATLAB/MATLAB_Compiler_Runtime/v717/
########################################################################################################################
## Use phantompeakqualtools to calculate quality tag (strand cross-correlation of peak)
# did the antibody-treatment enrich sufficiently so that the ChIP signal can be separated from the background signal?
## --> use https://code.google.com/p/phantompeakqualtools/
mcdir $baseDir/qc_phantompeaks
#Rscript /local/home/brandl/bin/phantompeakqualtools/run_spp.R
#Rscript run_spp.R <bamfile> -savp -out=<outfile>
ctrlBam=$baseDir/alignments_untrimmed/H3-3_Rep1_Lane1_Lib4453.bam
#DEBUG ctrlBam=$baseDir/alignments_untrimmed/H3K4me3_Rep2_Lane2_Lib4455.bam
for bamFile in $baseDir/alignments_untrimmed/*bam ; do
echo "processing $bamFile..."
# Rscript /home/brandl/bin/phantompeakqualtools/run_spp.R -c=$bamFile -savp -ou t=$(basename $bamFile).pt.log
## one output file is enough because output is tabular including inpuyt file name
# Rscript /home/brandl/bin/phantompeakqualtools/run_spp.R -c=$bamFile -i=$ctrlBam -savp -savd -odir=$(pwd) -out=phantom_qc.txt
## with control
# ( Rscript /home/brandl/bin/phantompeakqualtools/run_spp.R -c=$bamFile -savp -savd -odir=$(pwd) -out=phantom_qc_noctrl.txt &> $(basename $bamFile).pt.log ) &
## without control
# ( Rscript /home/brandl/bin/phantompeakqualtools/run_spp.R -c=$bamFile -savp -savd -odir=$(pwd) -out=phantom_qc_noctrl.txt &> $(basename $bamFile).pt.log ) &
done
wait
## todo visualize results
#phantom output format
#format:Filename<tab>numReads<tab>estFragLen<tab>corr_estFragLen<tab>PhantomPeak<tab>corr_phantomPeak<tab>argmin_corr<tab>min_corr<tab>Normalized SCC (NSC)<tab>Relative SCC (RSC)<tab>QualityTag)
mailme "phantom qc done"
########################################################################################################################
### motifs
## run meme-chip
########################################################################################################################
### profiles
## try deepTools
## differential binding diffbind