美文网首页复习相关
[实战1] Polycomb associates genome

[实战1] Polycomb associates genome

作者: shine_9457 | 来源:发表于2021-02-01 08:47 被阅读0次

    写在前面:参考-# 给学徒ChIP-seq数据处理流程

    目录
    软件安装
    [公共数据的获取---sratoolkit --prefetch ]
    [得到fastq格式的测序数据----sratoolkit fastq-dump]
    fastq数据的质量控制 ----fastqc;trim_galore 多参数
    比对以及质控 ----bowtie建立索引;samtool格式转换
    讨论去除PCR重复与否
    [批处理合并多条lane的测序数据 ---samtools的merge命令]
    使用macs2找peaks
    [IGV可视化]
    [IGV-bigWig和bedgraph文件详情]
    deeptools可视化
    使用R包注释
    [使用Homer找motif]

    流程-软件

    1 需要的软件

    下载,使用conda

    source activate chipseq
    # 可以用search先进行检索
    conda install -y sra-tools  
    conda install -y trim-galore  samtools
    conda install -y deeptools homer  meme
    conda install -y macs2 bowtie bowtie2 
    

    2 数据获取
    使用srr号+prefetch
    默认下载地址:/home/yangjy/ncbi/public/sra
    数据分析地址:/home/yangjy/chipseq_test/rawdata
    移动一下:使用mv

    (chipseq) [yangjy@GSCG01 ~]$ cd /home/yangjy/ncbi/public/sra
    (chipseq) [yangjy@GSCG01 sra]$ ls
    SRR391032.sra  SRR391033.sra  SRR391034.sra  SRR391041.sra  SRR391043.sra
    (chipseq) [yangjy@GSCG01 sra]$ mv *.sra /home/yangjy/chipseq_test/rawdata
    (chipseq) [yangjy@GSCG01 sra]$ ls #无文件,已移动
    (chipseq) [yangjy@GSCG01 sra]$ cd #回到家目录
    (chipseq) [yangjy@GSCG01 ~]$ cd /home/yangjy/chipseq_test/rawdata 
    (chipseq) [yangjy@GSCG01 rawdata]$ ls
    SRR391032.sra  SRR391033.sra  SRR391034.sra  SRR391041.sra  SRR391043.sra
    (chipseq) [yangjy@GSCG01 rawdata]$ ls -lh #查看细节
    total 2.3G
    -rw-rw-r--. 1 yangjy yangjy 474M Jan 31 11:21 SRR391032.sra
    -rw-rw-r--. 1 yangjy yangjy 473M Jan 31 12:15 SRR391033.sra
    -rw-rw-r--. 1 yangjy yangjy 406M Jan 31 12:17 SRR391034.sra
    -rw-rw-r--. 1 yangjy yangjy 322M Jan 31 12:29 SRR391041.sra
    -rw-rw-r--. 1 yangjy yangjy 597M Jan 31 12:32 SRR391043.sra
    

    批量建立文件夹 mkdir{sra,raw,qc,clean,align,peacks,motif}
    3 得到fastq格式的测序数据
    制作配置文件;即把srr号与基因组的名称对应起来:

    RNAPII_S5P_1    SRR391032
    RNAPII_S5P_2    SRR391033
    RNAPII_S2P_1    SRR391034
    H2Aub1_1            SRR391041
    H3K36me3_1       SRR391043
    

    用fastq-dump把sra格式转成fastq格式(fq格式)

    ls *sra |while read id; do ~/miniconda3/envs/chipseq/bin/fastq-dump $id;done
    

    实际情况

    (chipseq) [yangjy@GSCG01 rawdata]$ which fastq-dump
    ~/miniconda3/envs/chipseq/bin/fastq-dump
    (chipseq) [yangjy@GSCG01 rawdata]$ ls *sra |while read id; do ~/miniconda3/envs/chipseq/bin/fastq-dump $id;done
    Read 13532944 spots for SRR391032.sra
    Written 13532944 spots for SRR391032.sra
    Read 13662812 spots for SRR391033.sra
    Written 13662812 spots for SRR391033.sra
    Read 26571408 spots for SRR391034.sra
    Written 26571408 spots for SRR391034.sra
    Read 9780305 spots for SRR391041.sra
    Written 9780305 spots for SRR391041.sra
    Read 12313493 spots for SRR391043.sra
    Written 12313493 spots for SRR391043.sra
    (chipseq) [yangjy@GSCG01 rawdata]$ ls
    SRR391032.fastq  SRR391033.fastq  SRR391034.fastq  SRR391041.fastq  SRR391043.fastq
    SRR391032.sra    SRR391033.sra    SRR391034.sra    SRR391041.sra    SRR391043.sra
    (chipseq) [yangjy@GSCG01 rawdata]$ ls -lh
    total 18G
    -rw-rw-r--. 1 yangjy yangjy 2.6G Jan 31 14:56 SRR391032.fastq
    -rw-rw-r--. 1 yangjy yangjy 474M Jan 31 11:21 SRR391032.sra
    -rw-rw-r--. 1 yangjy yangjy 2.6G Jan 31 14:57 SRR391033.fastq
    -rw-rw-r--. 1 yangjy yangjy 473M Jan 31 12:15 SRR391033.sra
    -rw-rw-r--. 1 yangjy yangjy 5.0G Jan 31 14:58 SRR391034.fastq
    -rw-rw-r--. 1 yangjy yangjy 406M Jan 31 12:17 SRR391034.sra
    -rw-rw-r--. 1 yangjy yangjy 1.9G Jan 31 14:58 SRR391041.fastq
    -rw-rw-r--. 1 yangjy yangjy 322M Jan 31 12:29 SRR391041.sra
    -rw-rw-r--. 1 yangjy yangjy 2.9G Jan 31 14:59 SRR391043.fastq
    -rw-rw-r--. 1 yangjy yangjy 597M Jan 31 12:32 SRR391043.sra
    

    注意:
    这里的SRRxxx.sra格式转换后为*.fastq格式,用--gzip就能输出*.fastq.gz格式, 能够节省空间的同时也不会给后续比对软件造成压力, 比对软件都支持。
    可以使用

    ls *sra |while read id; do ~/miniconda3/envs/chipseq/bin/fastq-dump --gzip --split-3 $id -O ./fastq/;done
    

    结果:

    (chipseq) [yangjy@GSCG01 rawdata]$ cd fastq/
    (chipseq) [yangjy@GSCG01 fastq]$ ls -lh
    total 3.7G
    -rw-rw-r--. 1 yangjy yangjy 736M Jan 31 16:21 SRR391032.fastq.gz
    -rw-rw-r--. 1 yangjy yangjy 741M Jan 31 16:25 SRR391033.fastq.gz
    -rw-rw-r--. 1 yangjy yangjy 855M Jan 31 16:31 SRR391034.fastq.gz
    -rw-rw-r--. 1 yangjy yangjy 508M Jan 31 16:34 SRR391041.fastq.gz
    -rw-rw-r--. 1 yangjy yangjy 879M Jan 31 16:39 SRR391043.fastq.gz
    /home/yangjy/chipseq_test/fastq     移动后*.fastq.gz所在位置
    

    4 fastq数据的质量控制
    使用trim_galore软件进行质控
    先使用fastqc看一下;[对照]

    ls *.gz | while read id ; do /home/yangjy/miniconda3/envs/chipseq/bin/fastqc $id -O ../qc/;done
    

    可以先看一下质量:

    (chipseq) [yangjy@GSCG01 fastq]$ cd ../qc/
    (chipseq) [yangjy@GSCG01 qc]$ ls -lh
    total 2.8M
    -rw-rw-r--. 1 yangjy yangjy 244K Jan 31 17:34 SRR391032_fastqc.html# 打开
    -rw-rw-r--. 1 yangjy yangjy 330K Jan 31 17:34 SRR391032_fastqc.zip
    -rw-rw-r--. 1 yangjy yangjy 243K Jan 31 17:34 SRR391033_fastqc.html# 打开
    -rw-rw-r--. 1 yangjy yangjy 331K Jan 31 17:34 SRR391033_fastqc.zip
    -rw-rw-r--. 1 yangjy yangjy 220K Jan 31 17:35 SRR391034_fastqc.html# 打开
    -rw-rw-r--. 1 yangjy yangjy 272K Jan 31 17:35 SRR391034_fastqc.zip
    -rw-rw-r--. 1 yangjy yangjy 248K Jan 31 17:36 SRR391041_fastqc.html# 打开
    -rw-rw-r--. 1 yangjy yangjy 327K Jan 31 17:36 SRR391041_fastqc.zip
    -rw-rw-r--. 1 yangjy yangjy 265K Jan 31 17:37 SRR391043_fastqc.html# 打开
    -rw-rw-r--. 1 yangjy yangjy 366K Jan 31 17:37 SRR391043_fastqc.zip
    
    fastqc报告 [!黄色为警告;× 红色为报错;]

    这个时候选择trim_galore软件进行过滤,单端测序数据的代码如下;

    ls *.gz | while read id ; do 
    nohup 
    /home/yangjy/miniconda3/envs/chipseq/bin/trim_galore $id -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired fq --gzip -o ../clean/ & 
    done 
    

    较慢,于是放在后台处理。

    结果:
    (base) [yangjy@GSCG01 clean]$ ls -lh
    total 3.6G
    -rw-rw-r--. 1 yangjy yangjy   20 Jan 31 21:16 fq_trimmed.fq.gz
    -rw-rw-r--. 1 yangjy yangjy 1.5K Jan 31 21:16 fq_trimming_report.txt
    -rw-rw-r--. 1 yangjy yangjy 2.7K Jan 31 21:03 SRR391032.fastq.gz_trimming_report.txt
    -rw-rw-r--. 1 yangjy yangjy 703M Jan 31 21:03 SRR391032_trimmed.fq.gz
    -rw-rw-r--. 1 yangjy yangjy 2.7K Jan 31 21:06 SRR391033.fastq.gz_trimming_report.txt
    -rw-rw-r--. 1 yangjy yangjy 709M Jan 31 21:06 SRR391033_trimmed.fq.gz
    -rw-rw-r--. 1 yangjy yangjy 2.4K Jan 31 21:10 SRR391034.fastq.gz_trimming_report.txt
    -rw-rw-r--. 1 yangjy yangjy 843M Jan 31 21:10 SRR391034_trimmed.fq.gz
    -rw-rw-r--. 1 yangjy yangjy 2.7K Jan 31 21:12 SRR391041.fastq.gz_trimming_report.txt
    -rw-rw-r--. 1 yangjy yangjy 491M Jan 31 21:12 SRR391041_trimmed.fq.gz
    -rw-rw-r--. 1 yangjy yangjy 3.1K Jan 31 21:16 SRR391043.fastq.gz_trimming_report.txt
    -rw-rw-r--. 1 yangjy yangjy 844M Jan 31 21:16 SRR391043_trimmed.fq.gz
    

    再次进行QC,看看处理情况;尤其时质量差的时候

    5 比对以及质控 ----bowtie建立索引;samtool格式转换
    参考基因组bowtie2_indexes官网
    参考基因组及注释下载

    下载小鼠 [fasta格式的压缩文件] (迅雷下载较快,使用winSCP上传)
     解压(tar)合并(cat)建立索引(bowtie2-build)
    得到---> Output files: "mm10.*.bt2"  # 输出的文件格式:mm10.*.bt2
      inflating: mm10.1.bt2
      inflating: mm10.2.bt2
      inflating: mm10.3.bt2
      inflating: mm10.4.bt2
      inflating: mm10.rev.1.bt2
      inflating: mm10.rev.2.bt2
      inflating: make_mm10.sh
    

    使用命令:
    cd 到align目录下(相对于目录)

    ls ../clean/*.gz |while read id;
    do
    file=$(basename $id )
    sample=${file%%.*}
    echo $file $sample
    bowtie2 -p 5 -x /home/yangjy/project/epi/align/reference/mm10 -U  $id | samtools sort  -O bam  -@ 5 -o - > ${sample}.bam
    done
    

    参数解读:
    -p 5 指5个线程;
    -x 索引文件,不是参考基因组本身,是使用bowtie built 后的文件,参考;-x 后面接着gene-index,需要指定路径(/home/yangjy/project/epi/align/reference/)及其共用文件名(mm10)
    -U 后面接着fastq格式的文件(单末端测序参数)
    注:命令把bowtie2 生成的sam文件通过管道|传递到samtools,将sam转换为bam文件,省去中间sam文件的空间占用
    双末端测序使用:
    bowtie2 -p 5 -3 5 --local -x mm10 -1 example_1.fastq -2 example_2.fastq -S example.sam
    然后再将SAM 文件转为 BAM 文件:
    samtools sort example.sam > example.bam

    处理中:
    H2Aub1_1_trimmed.fq H2Aub1_1_trimmed   # 处理输出名称
    9208673 reads; of these:
      9208673 (100.00%) were unpaired; of these:
        239111 (2.60%) aligned 0 times
        6227498 (67.63%) aligned exactly 1 time
        2742064 (29.78%) aligned >1 times
    97.40% overall alignment rate
    [bam_sort_core] merging from 0 files and 5 in-memory blocks....
    
    结果:
    (chipseq) [yangjy@GSCG01 align]$ ls -l
    total 3342040
    drwxrwxr-x. 2 yangjy yangjy       122 Feb  2 20:51 index
    -rw-rw-r--. 1 yangjy yangjy 526164473 Jan 21 14:16 H2Aub1_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 873489325 Jan 21 14:26 H3K36me3_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 765203056 Jan 21 14:51 RNAPII_S2P_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 620491240 Jan 21 14:59 RNAPII_S5P_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 636886282 Jan 21 15:02 RNAPII_S5P_2_trimmed.bam
    

    6 对bam文件进行QC----- samtools 【samtools flagstat 】

    cd ~/project/epi/align
    ls  *.bam  |xargs -i samtools index {}   # 需要构建索引
    ls  *.bam  | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
    

    nohup 命令&
    得到:后缀为.bai的文件

    -rw-rw-r--. 1 yangjy yangjy 502M Jan 21 14:16 H2Aub1_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 3.0M Feb  2 21:44 H2Aub1_1_trimmed.bam.bai
    -rw-rw-r--. 1 yangjy yangjy 834M Jan 21 14:26 H3K36me3_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 2.8M Feb  2 21:44 H3K36me3_1_trimmed.bam.bai
    drwxrwxr-x. 2 yangjy yangjy  122 Feb  2 20:51 index
    -rw-rw-r--. 1 yangjy yangjy 730M Jan 21 14:51 RNAPII_S2P_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 2.6M Feb  2 21:44 RNAPII_S2P_1_trimmed.bam.bai
    -rw-rw-r--. 1 yangjy yangjy 592M Jan 21 14:59 RNAPII_S5P_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 3.7M Feb  2 21:44 RNAPII_S5P_1_trimmed.bam.bai
    -rw-rw-r--. 1 yangjy yangjy 608M Jan 21 15:02 RNAPII_S5P_2_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 3.7M Feb  2 21:44 RNAPII_S5P_2_trimmed.bam.bai
    

    以及:后缀为.stat的文件

    total 3358004
    -rw-rw-r--. 1 yangjy yangjy 526164473 Jan 21 14:16 H2Aub1_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy   3109504 Feb  2 21:44 H2Aub1_1_trimmed.bam.bai
    -rw-rw-r--. 1 yangjy yangjy       386 Feb  2 21:48 H2Aub1_1_trimmed.stat
    -rw-rw-r--. 1 yangjy yangjy 873489325 Jan 21 14:26 H3K36me3_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy   2912232 Feb  2 21:44 H3K36me3_1_trimmed.bam.bai
    -rw-rw-r--. 1 yangjy yangjy       388 Feb  2 21:48 H3K36me3_1_trimmed.stat
    drwxrwxr-x. 2 yangjy yangjy       122 Feb  2 20:51 index
    -rw-rw-r--. 1 yangjy yangjy 765203056 Jan 21 14:51 RNAPII_S2P_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy   2697544 Feb  2 21:44 RNAPII_S2P_1_trimmed.bam.bai
    -rw-rw-r--. 1 yangjy yangjy       388 Feb  2 21:48 RNAPII_S2P_1_trimmed.stat
    -rw-rw-r--. 1 yangjy yangjy 620491240 Jan 21 14:59 RNAPII_S5P_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy   3805200 Feb  2 21:44 RNAPII_S5P_1_trimmed.bam.bai
    -rw-rw-r--. 1 yangjy yangjy       388 Feb  2 21:48 RNAPII_S5P_1_trimmed.stat
    -rw-rw-r--. 1 yangjy yangjy 636886282 Jan 21 15:02 RNAPII_S5P_2_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy   3789176 Feb  2 21:44 RNAPII_S5P_2_trimmed.bam.bai
    -rw-rw-r--. 1 yangjy yangjy       388 Feb  2 21:48 RNAPII_S5P_2_trimmed.stat
    

    7 讨论去除PCR重复与否
    参考文章:
    仔细探究picard的MarkDuplicates 是如何行使去除PCR重复reads功能的
    仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的

    在bam文件夹下:
    ls  *.bam  | while read id ;do (nohup samtools markdup -r $id $(basename $id ".bam").mergeBam  & );done
    

    得到:

    (chipseq) [yangjy@GSCG01 align]$ ls -l *.rmdup.bam
    -rw-rw-r--. 1 yangjy yangjy 474047646 Feb  2 21:59 H2Aub1_1_trimmed.rmdup.bam
    -rw-rw-r--. 1 yangjy yangjy 830716373 Feb  2 22:00 H3K36me3_1_trimmed.rmdup.bam
    -rw-rw-r--. 1 yangjy yangjy 632581765 Feb  2 22:00 RNAPII_S2P_1_trimmed.rmdup.bam
    -rw-rw-r--. 1 yangjy yangjy 368497642 Feb  2 21:59 RNAPII_S5P_1_trimmed.rmdup.bam
    -rw-rw-r--. 1 yangjy yangjy 371445769 Feb  2 21:59 RNAPII_S5P_2_trimmed.rmdup.bam
    

    通过查看*.stat 文件:

    (chipseq) [yangjy@GSCG01 align]$ cat H2Aub1_1_trimmed.stat  #使用cat查看
    9208673 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 secondary                                            0条比对到第二个地方;
    0 + 0 supplementary
    0 + 0 duplicates
    8969562 + 0 mapped (97.40% : N/A)               #比对率97.40% 
    0 + 0 paired in sequencing 
    0 + 0 read1                             因为是单端
    0 + 0 read2
    0 + 0 properly paired (N/A : N/A)   
    0 + 0 with itself and mate mapped
    0 + 0 singletons (N/A : N/A)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)
    

    来自于网页(https://www.biostars.org/p/333465/

    45084184 + 0 in total (QC-passed reads + QC-failed reads)
    4717987 + 0 secondary
    0 + 0 supplementary
    0 + 0 duplicates
    45084184 + 0 mapped (100.00%:-nan%)
    40366197 + 0 paired in sequencing
    20273012 + 0 read1
    20093185 + 0 read2
    39146254 + 0 properly paired (96.98%:-nan%)
    39644363 + 0 with itself and mate mapped
    721834 + 0 singletons (1.79%:-nan%)
    260722 + 0 with mate mapped to a different chr
    235361 + 0 with mate mapped to a different chr (mapQ>=5)
    

    解读:
    QC-passed reads的数量为9208673,未通过的0;意味着一共有9208673条reads
    下面都是0?? 可能是单端测序的原因(双末端测序下载应该有一侧有数值)
    注: 可以使用 multiqc 对 *.stat 进行可视化

    可以再对*.rmdup.bam 进行统计:【samtools flagstat 】

    ls  *.rmdup.bam  |xargs -i samtools index {}   只有建立索引,才能进行--统计 ;得到后缀为.bai 的文件
    ls  *.rmdup.bam  | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
    

    综上两个*.bam 文件:

    (chipseq) [yangjy@GSCG01 align]$ ls -lh *.bam
    -rw-rw-r--. 1 yangjy yangjy 502M Jan 21 14:16 H2Aub1_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 453M Feb  2 21:59 H2Aub1_1_trimmed.rmdup.bam #去除PCR重复的 
    -rw-rw-r--. 1 yangjy yangjy 834M Jan 21 14:26 H3K36me3_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 793M Feb  2 22:00 H3K36me3_1_trimmed.rmdup.bam
    -rw-rw-r--. 1 yangjy yangjy 730M Jan 21 14:51 RNAPII_S2P_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 604M Feb  2 22:00 RNAPII_S2P_1_trimmed.rmdup.bam
    -rw-rw-r--. 1 yangjy yangjy 592M Jan 21 14:59 RNAPII_S5P_1_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 352M Feb  2 21:59 RNAPII_S5P_1_trimmed.rmdup.bam
    -rw-rw-r--. 1 yangjy yangjy 608M Jan 21 15:02 RNAPII_S5P_2_trimmed.bam
    -rw-rw-r--. 1 yangjy yangjy 355M Feb  2 21:59 RNAPII_S5P_2_trimmed.rmdup.bam
    

    由此,可以对这些*.bam文件进行peak calling
    此前,对【一个样品分成了多个lane进行测序】情况可以把bam进行合并(merge)。
    使用:samtools merge 参数;
    参看比对率: grep 'N/A' *.stat |grep %

    (chipseq) [yangjy@GSCG01 align]$  grep 'N/A' *.stat |grep %
    H2Aub1_1_trimmed.rmdup.stat:7744878 + 0 mapped (97.01% : N/A)
    H2Aub1_1_trimmed.stat:8969562 + 0 mapped (97.40% : N/A)
    H3K36me3_1_trimmed.rmdup.stat:11007740 + 0 mapped (98.82% : N/A)
    H3K36me3_1_trimmed.stat:11737364 + 0 mapped (98.89% : N/A)
    RNAPII_S2P_1_trimmed.rmdup.stat:18461592 + 0 mapped (96.32% : N/A)
    RNAPII_S2P_1_trimmed.stat:25018815 + 0 mapped (97.26% : N/A)
    RNAPII_S5P_1_trimmed.rmdup.stat:6042684 + 0 mapped (95.30% : N/A)
    RNAPII_S5P_1_trimmed.stat:11801448 + 0 mapped (97.53% : N/A)
    RNAPII_S5P_2_trimmed.rmdup.stat:6123727 + 0 mapped (96.43% : N/A)
    RNAPII_S5P_2_trimmed.stat:12182337 + 0 mapped (98.18% : N/A)
    

    注意格式:常用的mapping软件bwa & sam格式简介


    samtools view

    (base) [yangjy@GSCG01 dum]$ samtools view H2Aub1.merge.rmdup.bam | cut -f 1 | sort | uniq -c | wc -l
    18090450
    

    8 MACS 找call peak
    有了*.bam文件后就可以进行call了

    macs2包含一系列的子命令,其中最主要的就是callpeak, 官方提供了使用实例

    macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
    

    各个参数的意义:
    -t: 实验组的输出结果
    -c: 对照组的输出结果
    -f: -t和-c提供文件的格式,可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一个。如果不提供这项,就是自动检测选择。
    -g: 基因组大小, 默认提供了hs, mm, ce, dm选项, 不在其中的话,比如说拟南芥,就需要自己提供了。
    -n: 输出文件的前缀名
    -B: 会保存更多的信息在bedGraph文件中,如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
    -q: q值,也就是最小的PDR阈值, 默认是0.05。q值是根据p值利用BH计算,也就是多重试验矫正后的结果。
    -p: 这个是p值,指定p值后MACS2就不会用q值了。
    -m: 和MFOLD有关,而MFOLD和MACS预构建模型有关,默认是5:50。

    处理:

    cd bam 文件目录
    使用命令:

    macs2 callpeak -c  Control_1_trimmed.bam -t H3K36me3_1_trimmed.bam -f BAM -B -g mm -n H3K36me3
    # macs2 callpeak -c <Control.bam> -t <test.bam> -f BAM -B -g mm -n testname
    

    注:chipseq环境下的macs一直报错;
    退出conda环境,使用命令---> 成功

    (base) [yangjy@GSCG01 align1]$ conda list macs2
    packages in environment at /home/yangjy/miniconda3:
    #########
    Name                    Version                   Build  Channel
    macs2                     2.2.7.1          py38h0213d0e_1    https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
    

    处理时的信息参考
    得到文件:

    -rw-rw-r--. 1 yangjy yangjy 565332244 Feb  3 20:40 H3K36me3_control_lambda.bdg
    -rw-rw-r--. 1 yangjy yangjy     99001 Feb  3 20:39 H3K36me3_model.r
    -rw-rw-r--. 1 yangjy yangjy    264716 Feb  3 20:40 H3K36me3_peaks.narrowPeak
    -rw-rw-r--. 1 yangjy yangjy    302627 Feb  3 20:40 H3K36me3_peaks.xls
    -rw-rw-r--. 1 yangjy yangjy    179908 Feb  3 20:40 H3K36me3_summits.bed
    -rw-rw-r--. 1 yangjy yangjy 689424712 Feb  3 20:40 H3K36me3_treat_pileup.bdg
    

    文件解读:
    1 Excel文件
    前几行为运行命令以及信息
    后面从左到右分别是染色体,peak的起始终止位置,长度,峰值位置,高度,pvalue值记忆富集程度

    (base) [yangjy@GSCG01 peaks]$ cat H2Aub1_peaks.xls
    This file is generated by MACS version 2.1.0.20150731
    Command line: callpeak -c Control.merge.bam -t H2Aub1.merge.bam -f BAM -B -g mm -n H2Aub1 --outdir ../peaks
    ARGUMENTS LIST:
    name = H2Aub1
    format = BAM
    ChIP-seq file = ['H2Aub1.merge.bam']
    control file = ['Control.merge.bam']
    effective genome size = 1.87e+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
    tag size is determined as 49 bps
    total tags in treatment: 22199430
    tags after filtering in treatment: 17515747
    maximum duplicate tags at the same position in treatment = 1
    Redundant rate in treatment: 0.21
    total tags in control: 14660284
    tags after filtering in control: 12331048
    maximum duplicate tags at the same position in control = 1
    Redundant rate in control: 0.16
    d = 83
    alternative fragment length(s) may be 83 bps
    chr     start   end     length  abs_summit      pileup  -log10(pvalue)  fold_enrichment -log10(qvalue)  name
    chr1    4492658 4492755 98      4492668 11.97   9.08008 6.48398 4.02077 H2Aub1_peak_1
    chr1    4492851 4493214 364     4493158 14.78   16.26556        10.20087        8.94022 H2Aub1_peak_2
    chr1    4493287 4493599 313     4493468 9.15    9.39303 6.56103 4.02077 H2Aub1_peak_3
    chr1    4496519 4496628 110     4496554 8.45    7.41898 5.67787 2.72721 H2Aub1_peak_4
    chr1    4497051 4497263 213     4497103 9.86    9.39303 7.01601 4.02077 H2Aub1_peak_5
    

    2 BED文件
    染色体,peak的起始终止位置,

    (base) [yangjy@GSCG01 mergeBam]$ less RNAPII_S5P_summits.bed
    chr1    3150842 3150843 RNAPII_S5P_peak_1       8.17172
    chr1    3953635 3953636 RNAPII_S5P_peak_2       4.27895
    chr1    4085551 4085552 RNAPII_S5P_peak_3       9.22983
    ............................
    chrY    90834489        90834490        RNAPII_S5P_peak_56933   7.80872
    chrY    90835530        90835531        RNAPII_S5P_peak_56934   6.31523
    chrY    90837388        90837389        RNAPII_S5P_peak_56935   6.23419
    chrY    90838538        90838539        RNAPII_S5P_peak_56936   12.47469
    chrY    90839706        90839707        RNAPII_S5P_peak_56937   5.36220
    chrY    90840114        90840115        RNAPII_S5P_peak_56938   6.03569
    

    可以使用wc -l *.bed 查看一下每个基因有多少个peaks

           0 Control_summits.bed
        1089 H2Aub1_summits.bed
       40846 H3K36me3_summits.bed
       26027 Ring1B_summits.bed
       41816 RNAPII_8WG16_summits.bed
       19973 RNAPII_S2P_summits.bed
       56942 RNAPII_S5P_summits.bed
       72262 RNAPII_S7P_summits.bed
    

    MockIP是control,所以它自己跟自己比较,是没有peaks的;

    NAMEpeaks.xls: 以表格形式存放peak信息,虽然后缀是xls,但其实能用文本编辑器打开,和bed格式类似,但是以1为基,而bed文件是以0为基.也就是说xls的坐标都要减一才是bed文件的坐标
    NAMEpeaks.narrowPeak NAMEpeaks.broadPeak 类似。后面4列表示为, integer score for display, fold-change,-log10pvalue,-log10qvalue,relative summit position to peak start。内容和NAMEpeaks.xls基本一致,适合用于导入R进行分析
    NAMEsummits.bed:记录每个peak的peak summits,也就是记录极值点的位置。MACS建议用该文件寻找结合位点的motif
    NAME_model.r,能通过NAME_model.r作图,得到是基于你提供数据的peak模型
    参考:https://www.jianshu.com/p/2b8e2ea26665


    9 将BAM文件转化为bw文件,使用deeptools

    deeptools提供bamCoveragebamCompare进行格式转换,为了能够比较不同的样本,需要对先将基因组分成等宽分箱(bin),统计每个分箱的read数,最后得到描述性统计值。对于两个样本,描述性统计值可以是两个样本的比率,或是比率的log2值,或者是差值。如果是单个样本,可以用SES方法进行标准化

    bamCoverage的基本用法
    `bamCoverage -b test.bam -o test.bw

     activate chipseq
    bamCoverage -e 170 -bs 10 -b ap2_chip_rep1_2_sorted.bam -o ap2_chip_rep1_2.bw
    # ap2_chip_rep1_2_sorted.bam是前期比对得到的BAM文件
    得到的bw文件就可以送去IGV/Jbrowse进行可视化。
    参数使用 -e/--extendReads和-bs/--binSize 即拓展原来的read长度,且设置分箱的大小
    

    如果需要以100为分箱,并且标准化到1x,且仅统计某一条染色体区域的正链,输出格式为bedgraph

    bamCoverage -e 170 -bs 100 -of bedgraph -r Chr4:12985884:12997458 --normalizeTo1x 100000000 -b 02-read-alignment/ap2_chip_rep1_1_sorted.bam -o chip.bedgraph
    bamCompare和bamCoverage类似
    需要提供两个样本,并且采用SES方法进行标准化,多了--ratio参数。
    

    首先把bam文件转为bw文件,详情:wig、bigWig和bedgraph文件详解

    cd  bam文件目录下
     activate chipseq
    ls  *.bam  |xargs -i samtools index {}  # 如果没有索引,需要先建立索引
    ls *.bam |while read id;do
    nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.bw & 
    done 
    

    参数 --normalizeUsing CPM


    以上的bed,bdg 以及bw 文件可以载入IGV进行可视化;
    head *.bed 查看
    通过:awk '{print $1":"$2"-"$3}' RNAPII_S5P_summits.bed |head可以在屏幕显示出 chr 和坐标信息;输入IGV来检查这个位置的peaks
    bw文件反应bam的测序深度


    10 查看TSS附件信号强度:使用deeptools
    需要下载参考的注释文件tss的bed格式
    下载方法
    首先对单一样本绘图:

    使用computeMatrix函数对tss上下10kb(10000)计算信号强度
    mkdir -p  ~/project/epi/tss 
    cd  ~/project/epi/tss 
    computeMatrix reference-point  --referencePoint TSS  -p 15  \
    -b 10000 -a 10000    \
    -R /home/yangjy/project/epi/peaks/mm10.tss.bed  \
    -S /home/yangjy/project/epi/align/mergeBam/H2Aub1.bw  \
    --skipZeros  -o matrix1_test_TSS.gz  \
    --outFileSortedRegions regions1_test_genes.bed  输出文件
    进行批量处理:
    ## 
    both plotHeatmap and plotProfile will use the output from   computeMatrix
    plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.png
    plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.pdf --plotFileFormat pdf  --dpi 720  
    plotProfile -m matrix1_test_TSS.gz  -out test_Profile.png
    plotProfile -m matrix1_test_TSS.gz  -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720 
    
    bed=/home/yangjy/project/epi/peaks/mm10.tss.bed
    for id in /home/yangjy/project/epi/align/mergeBam/*.bw ;
    do
    echo $id
    file=$(basename $id )
    sample=${file%%.*}
    echo $sample
    computeMatrix reference-point  --referencePoint TSS  -p 15  \
    -b 2000 -a 2000    \
    -R $bed   \
    -S $id  \
    --skipZeros  -o matrix1_${sample}_TSS_2K.gz  \
    --outFileSortedRegions regions1_${sample}_TSS_2K.bed
    plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.png
    plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.pdf --plotFileFormat pdf  --dpi 720
    plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.png
    plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720
    done
    

    输出:

    -rw-rw-r--. 1 yangjy yangjy 10108967 Jan 27 09:46 Ring1B_Heatmap_10K.pdf
    -rw-rw-r--. 1 yangjy yangjy  1093343 Jan 27 09:44 Ring1B_Heatmap_10K.png
    -rw-rw-r--. 1 yangjy yangjy  4434818 Jan 27 00:33 Ring1B_Heatmap_2K.pdf
    -rw-rw-r--. 1 yangjy yangjy   639686 Jan 27 00:32 Ring1B_Heatmap_2K.png
    -rw-rw-r--. 1 yangjy yangjy   392676 Jan 27 09:47 Ring1B_Profile_10K.pdf
    -rw-rw-r--. 1 yangjy yangjy    38805 Jan 27 09:46 Ring1B_Profile_10K.png
    -rw-rw-r--. 1 yangjy yangjy   385755 Jan 27 00:33 Ring1B_Profile_2K.pdf
    -rw-rw-r--. 1 yangjy yangjy    35869 Jan 27 00:33 Ring1B_Profile_2K.png
    

    图片:


    Control_Profile_2K.png

    为了统计全基因组范围的peak在基因特征的分布情况,也就是需要用到computeMatrix计算,用plotHeatmap以热图的方式对覆盖进行可视化,用plotProfile以折线图的方式展示覆盖情况。

    computeMatrix具有两个模式:scale-region和reference-point。前者用来信号在一个区域内分布,后者查看信号相对于某一个点的分布情况。无论是那个模式,都有有两个参数是必须的,-S是提供bigwig文件,-R是提供基因的注释信息。还有更多个性化的可视化选项。

    11 使用R包对找到的peaks文件进行注释(还不熟悉)
    12 peaks相关基因集的注释

    homer软件来 寻找motif
    homer软件处理的原理:参考

    先准备好格式化文件,再使用命令;
    1 提供包含基因组坐标的文件 比如 peak文件或BED文件
    2 分析peak文件中富集的motif:
    findMotifsGenome.pl <peak/BED file> <genome> <output directory> -size

    • 简单例子:
      findMotifsGenome.pl ERpeaks.txt hg18 ER_MotifOutput/ -size 200 -mask
      -mask 使用repeated-mask序列
      -size 设置motif长度
      -len : motif的长度设置,默认8,10,12,越大越消耗计算资源。

    • 处理后得到的文件

    homerMotifs.motifs <#> : these are the output files from the de novo motif finding, separated by motif length, and represent separate runs of the algorithm.
    homerMotifs.all.motifs : Simply the concatenated file composed of all the homerMotifs.motifs<#> files.
    motifFindingParameters.txt :         记录执行findMotifsGenome.pl的命令
    knownResults.txt :                   记录motif的统计数据,text file(open in EXCEL).
    seq.autonorm.tsv :                    autonormalization statistics for lower-order oligo normalization.
    homerResults.html :                   de novo motif finding的格式化输出
    
    • 文件解读:
      1 txt文件: 记录执行findMotifsGenome.pl的命令以及motif的统计数据,
    (chipseq) [yangjy@GSCG01 motif]$ cat H2Aub1.annLog.txt
    
            Peak file = homer_peaks.tmp
            Genome = mm10
            Organism = mouse
            Peak/BED file conversion summary:
                    BED/Header formatted lines: 0
                    peakfile formatted lines: 1089
                    Duplicated Peak IDs: 0
    
            Peak File Statistics:
                    Total Peaks: 1089
                    Redundant Peak IDs: 0
                    Peaks lacking information: 0 (need at least 5 columns per peak)
                    Peaks with misformatted coordinates: 0 (should be integer)
                    Peaks with misformatted strand: 0 (should be either +/- or 0/1)
    
            Peak file looks good!
    
            Reading Positions...
            -----------------------
            Finding Closest TSS...
            Annotating:....................
                    Annotation      Number of peaks Total size (bp) Log2 Ratio (obs/exp)    LogP enrichment (+values depleted)
                    3UTR    38.0    20746162        2.146   -29.987
                    miRNA   0.0     31126   -0.018  0.013
                    ncRNA   22.0    3439613 3.950   -42.243
                    TTS     20.0    27176332        0.830   -4.492
                    pseudo  0.0     540203  -0.291  0.224
                    Exon    205.0   34540955        3.842   -376.610
                    Intron  400.0   947608369       0.028   -1.128
                    Intergenic      163.0   1564039797      -1.990  464.277
                    Promoter        183.0   30063639        3.878   -339.028
                    5UTR    58.0    2354993 5.895   -184.382
                    snoRNA  0.0     19      -0.000  0.000
                    rRNA    0.0     5631    -0.003  0.002
            NOTE: If this part takes more than 2 minutes, there is a good chance
                    your machine ran out of memory: consider hitting ctrl+C and rerunning
                    the command with "-noann"
            To capture annotation stats in a file, use "-annStats <filename>" next time
            Annotating:....................
                    Annotation      Number of peaks Total size (bp) Log2 Ratio (obs/exp)    LogP enrichment (+values depleted)
                    3UTR    38.0    20746162        2.146   -30.001
                    Other   0.0     7154386 -1.989  2.964
                    RC?     0.0     10979   -0.007  0.005
                    RNA     0.0     113962  -0.066  0.047
                    miRNA   0.0     31126   -0.018  0.013
                    ncRNA   22.0    3439613 3.950   -42.252
                    TTS     20.0    27176332        0.831   -4.497
                    LINE    0.0     521060815       -8.076  240.272
                    LINE?   0.0     8168    -0.005  0.003
                    srpRNA  0.0     43307   -0.026  0.018
                    SINE    0.0     193760064       -6.452  83.282
                    RC      0.0     65629   -0.039  0.027
                    tRNA    0.0     266147  -0.151  0.110
                    DNA?    0.0     142070  -0.082  0.059
                    pseudo  0.0     540203  -0.291  0.224
                    DNA     1.0     28358286        -3.553  9.244
                    Exon    205.0   34540955        3.842   -376.698
                    Intron  243.0   601725627       -0.035  1.055
                    Intergenic      103.0   784206321       -1.656  135.581
                    Promoter        183.0   30063639        3.879   -339.106
                    5UTR    58.0    2354993 5.895   -184.408
                    snoRNA  0.0     19      -0.000  0.000
                    LTR?    0.0     192563  -0.111  0.080
                    scRNA   0.0     577291  -0.309  0.239
                    CpG-Island      206.0   3107258 7.324   -865.124
                    Low_complexity  1.0     18489411        -2.936  5.514
                    LTR     3.0     292929648       -5.336  115.531
                    Simple_repeat   6.0     55283218        -1.931  10.524
                    snRNA   0.0     236859  -0.135  0.098
                    Unknown 0.0     1234764 -0.596  0.511
                    SINE?   0.0     29758   -0.018  0.012
                    Satellite       0.0     3685416 -1.338  1.526
                    rRNA    0.0     165655  -0.096  0.069
            Counting Tags in Peaks from each directory...
            Organism: mouse
            Loading Gene Informaiton...
    readline() on closed filehandle IN at /home/yangjy/miniconda3/envs/chipseq/bin/annotatePeaks.pl line 3741, <IN> line 2178.
    readline() on closed filehandle IN at /home/yangjy/miniconda3/envs/chipseq/bin/annotatePeaks.pl line 3751.
            Outputing Annotation File...
            Done annotating peaks file
    

    2 xls文件:记录染色体,起始位点,链+/-, peaks落在什么位置(UTR, promoter-TSS,exon , intron),距TSS最近PromoterID的注释距离...等。

    (chipseq) [yangjy@GSCG01 motif]$ head H2Aub1.peakAnn.xls
    PeakID (cmd=annotatePeaks.pl homer_peaks.tmp mm10)      Chr     Start   End     Strand  Peak Score      Focus Ratio/Region Size Annotation      Detailed Annotation  Distance to TSS Nearest PromoterID      Entrez ID       Nearest Unigene Nearest Refseq  Nearest Ensembl Gene Name       Gene Alias      Gene Description     Gene Type
    H2Aub1.log_peak_2       chr1    4493157 4493158 +       0       NA      exon (NM_001289464, exon 3 of 4)        exon (NM_001289464, exon 3 of 4)        4197NM_001289464
    H2Aub1.log_peak_62      chr10   19849898        19849899        +       0       NA      intron (NM_029529, intron 1 of 1)       intron (NM_029529, intron 1 of 1)    1561    NM_029529
    H2Aub1.log_peak_141     chr11   96354139        96354140        +       0       NA      promoter-TSS (NR_155303)        promoter-TSS (NR_155303)        -745NR_155303
    H2Aub1.log_peak_1020    chr9    32542410        32542411        +       0       NA      promoter-TSS (NM_008026)        promoter-TSS (NM_008026)        -958NM_008026
    H2Aub1.log_peak_269     chr13   112651870       112651871       +       0       NA      intron (NM_010029, intron 1 of 20)      intron (NM_010029, intron 1 of 20)   440     NM_001145885
    H2Aub1.log_peak_129     chr11   93099621        93099622        +       0       NA      5' UTR (NM_028296, exon 1 of 9) 5' UTR (NM_028296, exon 1 of 9) 331 NM_028296
    H2Aub1.log_peak_442     chr18   64346742        64346743        +       0       NA      intron (NM_194268, intron 1 of 1)       intron (NM_194268, intron 1 of 1)    6378    NM_194268
    

    3 homer*.motifs文件:
    每个motif信息是一块,均以>开头。>所在行的信息以tab分隔。 motif首行信息解释:

    >ACCTGCTTGAAA   19-ACCTGCTTGAAA 10.757736       -456.639967     0       T:120.0(0.46%),B:1.8(0.01%),P:1e-198
    0.997   0.001   0.001   0.001
    0.001   0.995   0.001   0.003
    0.001   0.997   0.001   0.001
    0.001   0.001   0.001   0.997
    0.003   0.001   0.995   0.001
    0.001   0.997   0.001   0.001
    0.001   0.001   0.001   0.997
    0.003   0.001   0.003   0.993
    0.001   0.001   0.997   0.001
    

    一致性序列:如图上的>ACCTGCTTGAAA
    Motif名称:如图上的19-ACCTGCTTGAAA
    比值检测概率的log值:10.757736
    P-value得log值: -456.639967
    占位符:如上图得0,不具有任何信息

    逗号分隔得富集信息,如:T:120.0(0.46%),B:1.8(0.01%),P:1e-198

    其中富集信息的解读: T:120.0(0.46%),B:1.8(0.01%),P:1e-198

    • T表示带有该motif的目标序列在总的目标序列(target)中得百分比
    • B表示带有该motif的背景序列在总的背景序列(background)中得百分比
    • P表示最终富集的p-value

    逗号分隔的motif统计信息,如:Tpos:68.2,Tstd:51.5,Bpos:123.7,Bstd:43.3,StrandBias:3.0,Multiplicity:1.00

    • Tpos:motif在目标序列中的平均位置
    • Tstd:motif在目标序列中位置的标准差
    • Bpos:motif在背景序列中的平均位置
    • Bstd:motif在背景序列中位置的标准差
    • StrandBias: 链偏好性,在正义链上的motif数与反义链motif数的比值的log
    • Multiplicity:具有一个或多个结合位点的序列中每个序列的平均出现次数;(大于等于1)

    其他可视化软件:
    1 ngs.plot

    bam文件建立索引
    根据功能原件的基因组坐标索引bam文件
    计算功能原件区域富集信号丰度

    参考画图:富集轮廓图和热图

    2 将peakAnn.xls文件导入到本地,通过excel透视功能进行可视化

    ATAC参考:https://www.it610.com/article/1224820226680524800.htm
    callpeaks文件 解读: https://blog.csdn.net/sunyu_03/article/details/82633799
    热图软件 deeptools:https://blog.csdn.net/xiaomotong123/article/details/108734727

    相关文章

      网友评论

        本文标题:[实战1] Polycomb associates genome

        本文链接:https://www.haomeiwen.com/subject/wcoltltx.html