Xueyd's Blog

0%

MeRIP-seq数据分析之peak calling

发表于 分类于
本文字数： 1.1k 阅读时长 ≈ 1 分钟

m6A-seq上游的数据处理流程和Ribo-seq差不多，包括去接头、数据质控、比对。比对之前也有文献会对m6A-seq数据去rRNA、tRNA，根据需求处理即可，这里不再总结，直接从得到比对的结果之后开始。

1. 提取高质量比对结果

得到比对的结果out.sam之后，要提取比对质量高的reads进行peak calling。
-f 2:提取完全mapping上的序列；
-F 256：过滤掉比对到多个位置上的序列；
-q 30：提取比对的质量值在30以上的序列。

1
2
samtools view -@ 10 -bS out.sam -f 2 -F 256 -q 30 -o out.highQ.bam
samtools sort -@ 10 -l 9 out.highQ.bam -o out.highQ.sorted.bam

2. MACS2进行peak calling

MACS2是为ChIP-seq设计的，已经被广泛使用。主要参数：
-t：处理组样本bam文件，即IP样本；
-c：对照组样本bam文件，即input样本；
-B：save extended fragment pileup, and local lambda tracks (two files) at every bp into a bedGraph file；
–SPMR：this option is recommended for displaying normalized pileup tracks across many datasets,可以在样本间进行比较peaks，需要-B参数；
–keep-dup：保留重复的tag；
-n：输出的文件前缀；
–outdir：输出文件路径；
-f：文件格式，这里为BAM；
-g：基因组大小，可以是数字或物种，如：hs、mm、ce、dm；
–nomodel：不建立shifting model；
–extsize：与nomodel一起使用，从5’->3’延伸fragment到指定长度；
-q：qvalue;
–fe-cutoff：fold enrichment。

1
macs2 callpeak -t treatment.highQ.sorted.bam -c control.highQ.sorted.bam -B --SPMR --keep-dup all -n outname --outdir out_peak -f BAM -g mm --nomodel --extsize 150 -q 0.05 --fe-cutoff 2

3. peak calling结果文件

  1. __control_lambda.bdg/__treat_pileup.bdg
    Input和IP样本的bedgraph格式数据，可以载入IGV可视化。

  2. __peak.narrowPeak
    每一列信息如下：
    chrom start end peak_name peak_length +/-strand Foldchange -log10pvalue -log10qvalue peak_start位置距离summit的距离

    1
    2
    3
    chr1    3682490 3682640 MII-m6A-1_peak_50       212     .       6.90033 24.1656 21.2778 127
    chr1    3683230 3684604 MII-m6A-1_peak_51       9634    .       19.4958 967.498 963.478 454
    chr1    3741597 3741747 MII-m6A-1_peak_52       122     .       7.52285 15.0429 12.2639 48

  3. __peak.xls
    xls格式的peak数据。abs_summit为顶峰的位置，pileup为顶峰的值。

    1
    2
    3
    4
    chr     start   end     length  abs_summit      pileup  -log10(pvalue)  fold_enrichment -log10(qvalue)  name
    chr1	3682491 3682640 150     3682618 43.89   24.1656 6.90033 21.2778 MII-m6A-1_peak_50
    chr1    3683231 3684604 1374    3683685 1090.77 967.498 19.4958 963.478 MII-m6A-1_peak_51
    chr1    3741598 3741747 150     3741646 21.94   15.0429 7.52285 12.2639 MII-m6A-1_peak_52

  4. __summits.bed
    每个peak顶峰的位置，最后一列为-log10qvalue。

    1
    2
    3
    chr1    3682617 3682618 MII-m6A-1_peak_50       21.2778
    chr1    3683684 3683685 MII-m6A-1_peak_51       963.478
    chr1    3741645 3741646 MII-m6A-1_peak_52       12.2639

4. 生物学重复replicates样本的peaks取交集

取交集用bedtools intersect,取至少有50%区域overlsap的peak。

1
bedtools intersect -a rep1_peaks.narrowPeak -b rep2_peaks.narrowPeak -f 0.5 -F 0.5 -e > overlap.peak