Question: criteria for choosing proper parameters for MASC2 -- chip-seq ?
0
gravatar for jimmy_zeng
17 months ago by
jimmy_zeng90
University of Macau
jimmy_zeng90 wrote:

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 !

chip-seq macs • 782 views
ADD COMMENTlink modified 17 months ago by tangming20052.2k • written 17 months ago by jimmy_zeng90

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.

ADD REPLYlink written 17 months ago by jimmy_zeng90
1
gravatar for tangming2005
17 months ago by
tangming20052.2k
Houston/MD Anderson Cancer Center
tangming20052.2k wrote:

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

ADD COMMENTlink written 17 months ago by tangming20052.2k

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!

ADD REPLYlink written 17 months ago by Wet&DryImmunology200
0
gravatar for Wet&DryImmunology
17 months ago by
Japan
Wet&DryImmunology200 wrote:

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.

ADD COMMENTlink written 17 months ago by Wet&DryImmunology200

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 ?

ADD REPLYlink written 17 months ago by jimmy_zeng90

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.

ADD REPLYlink written 17 months ago by Wet&DryImmunology200

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

ADD REPLYlink written 17 months ago by jotan1.1k
0
gravatar for Wet&DryImmunology
17 months ago by
Japan
Wet&DryImmunology200 wrote:

@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.

ADD COMMENTlink written 17 months ago by Wet&DryImmunology200

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 ?

ADD REPLYlink written 17 months ago by jimmy_zeng90

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.

ADD REPLYlink modified 17 months ago • written 17 months ago by Wet&DryImmunology200

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

ADD REPLYlink written 17 months ago by jimmy_zeng90
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1020 users visited in the last hour