criteria for choosing proper parameters for MASC2 -- chip-seq ?

Question

[jimmy_zeng]

Hi, All:

I'm now tring to using MACS2 call peaks for CHIP-seq data. But most of the information for this software are base on MACS1.4 or lower.

So, I just do a small test for it, and it runs very fast, I have enough time and storage to try it. So, the question is : Is there any criteria for choosing proper parameters for MASC2,

ARGUMENTS LIST:

name = Xu_MUT_rep1

format = SAM

ChIP-seq file = ['SRR1042593.sam']

control file = ['SRR1042594.sam']

effective genome size = 2.70e+09

band width = 300

model fold = [5, 50]

qvalue cutoff = 5.00e-02

Larger dataset will be scaled towards smaller dataset.

Range for calculating regional lambda is: 1000 bps and 10000 bps

Broad region calling is off

Paired-End mode is off

For the code below, I got less than 100 peaks, that I don't know what's wrong with the parameters.

nohup time ~/.local/bin/macs2 callpeak -c SRR1042594.sorted.bam -t SRR1042593.sorted.bam -f BAM -B -g hs -n Xu_MUT_rep1 2>Xu_MUT_rep1.masc2.log &
nohup time ~/.local/bin/macs2 callpeak -c SRR1042596.sorted.bam -t SRR1042595.sorted.bam -f BAM -B -g hs -n Xu_MUT_rep2 2>Xu_MUT_rep2.masc2.log &
nohup time ~/.local/bin/macs2 callpeak -c SRR1042598.sorted.bam -t SRR1042597.sorted.bam -f BAM -B -g hs -n Xu_WT_rep1  2>Xu_WT_rep1.masc2.log &
nohup time ~/.local/bin/macs2 callpeak -c SRR1042600.sorted.bam -t SRR1042599.sorted.bam -f BAM -B -g hs -n Xu_WT_rep2  2>Xu_WT_rep2.masc2.log &

should I do some filter for the alignment files ?

## cat >run_bowtie2.sh
ls *.fastq | while read id ; 
do  
echo $id
~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -p 8 -x ~/biosoft/bowtie/hg19_index/hg19 -U $id   -S ${id%%.*}.sam  2>${id%%.*}.align.log; 
samtools view -bhS -q 30  ${id%%.*}.sam > ${id%%.*}.bam  
samtools sort   ${id%%.*}.bam ${id%%.*}.sort  ## prefix for the output      
samtools index ${id%%.*}.sorted.bam 
done

I don't know whether I should filter out the alignment quality less than 30 , I just do it.

And I don't know whether I should remove PCR duplicated reads or not !

[jimmy_zeng]

Then I change the p value criteria just like :

nohup time ~/.local/bin/macs2 callpeak -c SRR1042594.sam -t SRR1042593.sam -f SAM -p 0.01  -g hs -n Xu_MUT_rep1 2>Xu_MUT_rep1.masc2.log &
nohup time ~/.local/bin/macs2 callpeak -c SRR1042596.sam -t SRR1042595.sam -f SAM -p 0.01  -g hs -n Xu_MUT_rep2 2>Xu_MUT_rep2.masc2.log &
nohup time ~/.local/bin/macs2 callpeak -c SRR1042598.sam -t SRR1042597.sam -f SAM -p 0.01  -g hs -n Xu_WT_rep1  2>Xu_WT_rep1.masc2.log &
nohup time ~/.local/bin/macs2 callpeak -c SRR1042600.sam -t SRR1042599.sam -f SAM -p 0.01  -g hs -n Xu_WT_rep2  2>Xu_WT_rep2.masc2.log &

Then I can get more than 10,000 peaks , I know the criteria is too broad. But I just don't know how to deal it.

Answer 1

[Ming Tommy Tang]

I wrote something might be helpful for you https://github.com/crazyhottommy/ChIP-seq-analysis/blob/master/part1.3_MACS2_peak_calling_details.md

[Wet&DryImmunology]

Hi Tangming, I wondered whether you could have a quick look the question I asked. Visualization of ChIP-seq data (mapped reads) Basically, the original paper had unusually sharp peaks, and I could not recapitulate the results. why made me wonder whether they use some special tricks or I had something wrong. Tsk a lot!

Answer 2

[Wet&DryImmunology]

I guess my previous question was sort of relevant, you might find some useful information here enter link description here My suggestion is "try and see", see enter link description here, "Finding Enriched Regions of Variable Length", the instruction is based on Homer, but I think the general idea is the same for other peak callers: basically you play with different combination of parameters see which one work best for you.

[jimmy_zeng]

yeh, I didn't choose HOMER, just becaues it's a little complicate to install it, I will try it as soon as possible.

We have the same problem about how to deal with the alignment files , yes or right ? fiter the bad quality? remove PCR duplicated reads ?

[Wet&DryImmunology]

the peak caller does not have to be Homer, what I am trying to say is try different parameters, check the peaks by your self to see the peak calling is proper not by loading the bed file of called peaks in some genome browser such as IGV. start with the default parameter.

For filter, enter link description herehere people recommend to use -q 5 or 10. For "remove PCR duplicated reads ?" Sorry I have no idea.

[jotan]

I would remove the PCR duplicates. A lot of these will map to repetitive elements in the genome for ChIP seq.

Answer 3

[Wet&DryImmunology]

@Jotan 1982, could we assign the parameter to remove the PCR duplicates, if they will map to receptive elements, does that mean they map to multiple loci? (not uniquely mapped), How to extract unique mapped results from Bowtie2 bam results? here people seems to suggest that by setting -q 5or 10 "samtools view -bhS -q ", reads which not map uniquely will be eliminated.

[jimmy_zeng]

I've use the command : samtools view -bhS -q 30 ${id%%.}.sam > ${id%%.}.bam in my code

but I don't know why I should use it ?

[Wet&DryImmunology]

you type the command: "samtools view " then you will get something like: Usage: samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]

Options: -b output BAM ~ -q INT only include reads with mapping quality >= INT [0] ~

" -q INT only include reads with mapping quality >= INT [0]" if you use "samtools view -bhS -q 30" which means you filter out reads with mapping quality lower than 30, some people suggest that more than 5 or 10 will be good enough.

[jimmy_zeng]

That's exactly what confuse me .

Some guys in ENCODE project suggest me to only filter out reads with mapping quality lower than 30. While others suggest that more than 5 or 10 should be keep .

maybe it doest matter, I will have a try for both