Content Comparison

...

...

...

...

...

...

****NOTE: The questions for Session V  homework are more comprehensive to help you get a better feel for the workflow and will involve working with full datasets. Therefore, for this homework we suggest that you do not work in your home directory. Follow the instructions below to work on /n/scratch2.****

ChIP-seq Practice Exercises: Answer Key

# ChIP-Seq Analysis Workflow

1.  Create a directory on /n/scratch2 using your eCommons user ID as the directory name. Within that directory create a directory called HCFC1_chipseq.

$ mkdir -p /n/scratch2/your_eCommons_ID/HCFC1_chipseq

2. You have strong evidence that HCFC1 is the transcription co-factor that associates with your protein of interest. To confirm this hypothesis you need to find binding regions for HCFC1 and see if they overlap with your current list of regions. The ENCODE project has ChIP-seq data for HCFC1 available for download. The dataset is using a human liver cancer cell line HepG2 and the data contains 32 bp single end reads.  We have downloaded this data and made it available for you on Orchestra. 

a.  Setup a project directory structure within the HCFC1_chipseq directory as shown below and copy over the raw FASTQ files from /n/groups/hbctraining/ngs-data-analysis2016/chipseq/HCFC1/ into the appropriate directory:

...

$ mkdir fastQC bowtie2 unique_align macs2 sambamba

$ cp cp /n/groups/hbctraining/ngs-data-analysis2016/chipseq/HCFC1/* /n/scratch2/your_eCommons_ID/HCFC1_chipseq/data/raw_fastq

...

1. 1. 1. Run FASTQC on all FASTQ files
      2. Align reads with bowtie2 using the parameters we used in class. You will need to copy over the hg19 index files from /n/groups/shared_databases/igenome/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/
      3. Change alignment file format from SAM to BAM (can be done using samtools or sambamba)
      4. Sort the BAM file by read coordinate locations (can be done using sambamba or with samtools)
      5. Filter to keep only uniquely mapping reads (this will also remove any unmapped reads) using sambamba
      6. Index the final BAM file. This will be useful for visualization and QC (deepTools)

...

module load seq/fastqc/0.11.3

module load gcc/6.2.0 bowtie2/2.2.9 

module load samtools/1.3.1

export PATH=/optn/app/bcbio/centostools/bin:$PATH

# Set paths

baseDir=/n/scratch2/your_eCommons_ID/HCFC1_chipseq

# FastQC on raw file

...

# Align trimmed file 

echo "****Running Bowties2 Bowtie2 to align raw reads****"

bowtie2 -p 6 -q -x /n/groups/shared_databases/igenome/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome -U $fq -S ${baseDir}/results/bowtie2/${NAME}_aln_unsorted.sam

...

sambamba view -h --nthreads 1 -f bam -F "[XS] == null and not unmapped and not duplicate" ${baseDir}/results/sambamba/${NAME}_aln_sorted.bam > ${baseDir}/results/unique_align/${NAME}_aln.bam

...

c. Create a separate job submission script to run the shell script you created in b) on all four .fastq files. You have the option of running this in serial or parallel. Take a look at the lessons from the RNA-Seq session the automation lesson to help with setting up the job submission script. 

## IN PARALLEL

#!/bin/bash

for fq in /n/scratch2/your_eCommons_ID/HCFC1_chipseq/data/raw_fastq/*.fastq

dobsub -q priority

    sbatch -p short -t 0-12:00 -n 6 -

...

-job-name chipseq-analysis -o %j.out -e

...

%j.err

...

--wrap="sh chipseq_align_on_input_file.sh $fq"

sleep 1

done

d. Once the job submission script has finished running, use MACS2 to call peaks for each replicate using the input and treat BAM files and setting a q-value threshold of 0.05. How many peaks are being called for each replicate?

$ cd /n/scratch2/your_eCommons_ID/HCFC1_chipseq

$ macs2 callpeak -t results/unique_align/HCFC1_treat_rep1_aln.bam -c results/unique_align/HCFC1_input_rep1_aln.bam -f BAM -g hs -q 0.05 --bdg --outdir results/macs2 -n HCFC1-rep1

$ macs2 callpeak -t results/unique_align/HCFC1_treat_rep2_aln.bam -c results/unique_align/HCFC1_input_rep2_aln.bam -f BAM -g hs -q 0.05 --bdg --outdir results/macs2 -n HCFC1-rep2

...

Replicate 2 has 10899 peaks

e.  Sort each of the narrowPeak files using:

sort -k8,8nr HCFC1-rep1_peaks.narrowPeak > HCFC1-rep1_peaks_sorted.narrowPeak

sort -k8,8nr HCFC1-rep2_peaks.narrowPeak > HCFC1-rep2_peaks_sorted.narrowPeak

f.  Use phantompeakqualtools to check the quality of the ChIP signal for HCFC1-rep1 and HCFC1-rep2. Report the NSC, RSC, and QualityTag values. What can you conclude about the quality of the data based on these quality measures?

...

Values of NSC > 1.1 and RSC > 1 indicate ChIP experiments with good signal to noise for the peak calls. This is reflected in the QualityTag value of 1, which indicates 'High' quality ChIP data.

g. Use deepTools to explore the correlation between samples using a heatmap. How well do the read coverages for the replicate samples correlate? Upload the heatmap to your homework folder.

# Created a 'deeptools' directory within 'unique_align' folder and changed directories into it

# Running mutliBamSummary to create coverage estimates for each sample

multiBamSummary bins --ignoreDuplicates -p 6 --bamfiles ../*aln.bam -out deeptools_multiBAM.out.npz --outRawCounts readCounts.tab

# Creating heatmap

plotCorrelation --corData deeptools_multiBAM.out.npz --plotFile deepTools_heatmap.png --corMethod pearson --whatToPlot heatmap --labels Input_Rep1 Input_Rep2 Treat_Rep1 Treat_Rep2 --plotNumbers

The coverage profiles for the replicate samples correlate very well, and we see stark differences in coverage profiles between the input and ChIP groups.

Image Removed

h 

f. Sort each of the narrowPeak files using:

sort -k8,8nr HCFC1-rep1_peaks.narrowPeak > HCFC1-rep1_peaks_sorted.narrowPeak

sort -k8,8nr HCFC1-rep2_peaks.narrowPeak > HCFC1-rep2_peaks_sorted.narrowPeak

g. Use IDR to assess reproducibility between replicates using the sorted peaks in 2e. **NOTE: you would normally set up IDR on peaks called using a more liberal threshold, but for this question we will work with what we have.** 2f. Use the awk command from class to determine how many peaks meet and IDR of 0.05?

$ mkdir /n/scratch2/your_eCommons_ID/HCFC1_chipseq/results/idr

$ cd /n/scratch2/your_eCommons_ID/HCFC1_chipseq/results/

$ idr --samples macs2/HCFC1-rep1_peaks_sorted.narrowPeak macs2/HCFC1-rep2_peaks_sorted.narrowPeak \

...

awk '{if($12 > 1.3) print $0}' idr/HCFC1-idr | wc -l

ih.  SUGGESTION: Use the resulting These high confidence peaks from IDR can be used to explore downstream tools for functional enrichment (GREAT) and motif discovery (MEME suite). Explain how these tools might help in drawing conclusions about your original hypothesis.

 ## OPTIONAL: 

3. Use the full dataset available at /groups/hbctraining/ngs-data-analysis2016/rnaseq/full_dataset/ to perform the MOV10 analysis as listed below. Please make sure your working directory is NOT your home directory, and that it is in the group directory for your lab. 

a. Using the workflow and submission scripts we generated in class, parallelize the RNA-Seq analysis of the whole files. Start with QC and trimming, perform the alignment with STAR, and obtain gene counts. You will have to modify a few things, like the directories, the number of cores, the location of the index, and any other pertinent things. Use the STAR index available at /groups/shared_databases/igenome/Homo_sapiens/UCSC/hg19/Sequence/starIndex/.

Once you have completed all the steps on Orchestra, transfer the files over to your personal computer and perform the differential expression analysis on the counts data using DESeq2.  Make sure you submit your workflow on orchestra asjobs as scripts

Report the DE genes identified for the KD and OE samples, in comparison to the Irrel control, at an FDR threshold of 0.05, by uploading a text file with only significant genes and all the other columns output by DESeq2.

b. Use the raw data to run through Sailfish, followed by tximport and DESeq2. Report the DE genes identified for the KD and OE samples, in comparison to the Irrel control, at an FDR threshold of 0.05, by uploading a text file with only significant genes and all the other columns output by DESeq2.

Note: You can create a Sailfish index for yourself using the transcript fasta file here: /groups/bcbio/bcbio/genomes/Hsapiens/GRCh37/rnaseq/ref-transcripts.fa

c. Get the overlap between the 2 methods after separating the lists into UP-regulated and DOWN-regulated genes, and report the numbers and lists of genes for the comparisons below:

UP-regulated genes in KD identified at FDR 0.05 from Sailfish (2b) and the standard approach (2a)

Down-regulated genes in KD identified at FDR 0.05 from Sailfish (2b) and the standard approach (2a)

UP-regulated genes in OE identified at FDR 0.05 from Sailfish (2b) and the standard approach (2a)

...

Use GREAT to annotate the peaks and write those annotations to file. Take a look at binding site locations relative to the TSS, what does this tell us?

Evaluate the GO enrichment analysis from GREAT, what terms do you see over-represented?

i. **OPTIONAL:** Try motif analysis using the MEME suite and comment on the results you find.