## 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