##############################
#                            #
# MA570 - Chip-Seq practical #
#                            #
##############################

# In this practical we will me using data from the ENCODE project to analyze 2 different histone modifications.
# These are H3k4me3 which is a well known marker for euchromatin,
# And Hek9me3 which is a marker for heterochromatin.

# The data for this prac consists of 3 fastQ files (One for H3k4me3, H3k9me3 and control).
# These are found in /home/nextgen/data/chip_seq

# If any of the tools are running slow or taking too long, feel free to use the files in /home/nextgen/peter/chip_seq


############################
#                          #
# Step 1: Mapping the data #
#                          #
############################

# The first step is to map the raw reads to the human genome using Bowtie

bowtie -l 32 -p 10 --chunkmbs 500 -a --best --strata -m 1 --sam /home/nextgen/data/indexes/chr18 /home/nextgen/data/chip_seq/Helas3.H3k4me3.chr18.fastq Helas3.H3k4me3.sam
bowtie -l 32 -p 10 --chunkmbs 500 -a --best --strata -m 1 --sam /home/nextgen/data/indexes/chr18 /home/nextgen/data/chip_seq/Helas3.H3k9me3.chr18.fastq Helas3.H3k9me3.sam
bowtie -l 32 -p 10 --chunkmbs 500 -a --best --strata -m 1 --sam /home/nextgen/data/indexes/chr18 /home/nextgen/data/chip_seq/Helas3.Control.chr18.fastq Helas3.Control.sam

# Options:
# -l seed length = 32; i.e., the number of bases on the high-quality end of the read is chosen to be 32bp.
# -p number of processors to use.
# -a --best --strata: report all possible alignments for any read and sort/cluster them by mapping quality.
# -m options tells bowtie that we only want uniquely mapped reads. Reads which map > 1 times are discarded.


######################################
#                                    #
# Step 2: Processing the mapped data #
#                                    #
######################################

# Next, convert from SAM to BAM format, sort the alignment and remove potential PCR duplicates using samtools
# The coversio/sorting steps are combined here using the pipe, to speed things up.

samtools view -Sb Helas3.H3k4me3.sam | samtools sort - Helas3.H3k4me3
samtools view -Sb Helas3.H3k9me3.sam | samtools sort - Helas3.H3k9me3
samtools view -Sb Helas3.Control.sam | samtools sort - Helas3.Control

# Now, remove duplicates

samtools rmdup Helas3.H3k4me3.bam Helas3.H3k4me3.rmdup.bam
samtools rmdup Helas3.H3k9me3.bam Helas3.H3k9me3.rmdup.bam
samtools rmdup Helas3.Control.bam Helas3.Control.rmdup.bam


###################################
#                                 #
# Step 3: Peak calling using MACS #
#                                 #
###################################

macs14 -t Helas3.H3k4me3.rmdup.bam -c Helas3.Control.rmdup.bam -g 70000000 -f BAM -n Helas3.H3k4me3
macs14 -t Helas3.H3k9me3.rmdup.bam -c Helas3.Control.rmdup.bam -g 70000000 -f BAM -n Helas3.H3k9me3

# The options used by MACS here are as follows:
# -t The treatment file, here the bam files for the histone modifications
# -c The control file
# -g The effective genome size. Here we set 70000000 which corresponds to 90% of the length of human chr18
#	 This value depends on the length of read and type of species, generally we choose between 80-90% of the genome length.
# -f The format of the input sequences. G=Here it is BAM format. 
# -n Name of the MACS run.


######################################
#                                    #
# Step 4: Visualise results with IGV #
#                                    #
######################################

# MACS produces a number of output files. 
# Here we are interested in the peaks.bed files which can be used by IGV and show us the peaks called
# Open IGV by typing igv and load these files.
# Dont forget fo select chr18.