search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
110 stars 30 forks source link

Phantompeakqualtools — 质量评估之cross correlation分析 #5114

Closed ixxmu closed 3 months ago

ixxmu commented 3 months ago

https://mp.weixin.qq.com/s/VaP1apm3ERu6oZrLjP2NVw

ixxmu commented 3 months ago

Phantompeakqualtools — 质量评估之cross correlation分析 by 生信菜鸟团

工欲善其事必先利其器


链交叉相关性分析是ChIP-seq质量评估的一个重要指标，基于高质量ChIP-seq实验会在蛋白质结合位点周围产生显著的富集DNA序列标签富集这一事实。这些“真实信号”的序列标签位于结合位点中心的一定距离处，这个距离取决于片段大小分布（Kharchenko等人，2008年）。这使得开发出一种基于基因组范围内链密度相关性的片段聚集（IP富集）量度成为可能。这种量度是通过计算Crick链和Watson链之间的皮尔逊线性相关性来实现的，计算时会将Watson链移动k个碱基对。这通常在交叉相关图中产生两个峰：一个对应于主要片段长度的富集峰，另一个对应于读长的峰（“幽灵”峰）。

1phantompeakqualtools

phantompeakqualtools是由 Anshul Kundaje开发的一款用于评估ChIP-seq和相关高通量测序数据的质量的工具，其具有以下功能

  • 信号与噪声评估：Phantompeakqualtools可以计算样本中的信号与噪声比例，帮助研究人员评估数据的质量。这包括计算如NSC（Normalized Strand Coefficient）和RSC（Relative Strand Coefficient）等关键指标，这些指标能够反映样本中信号的清晰度和显著性。
  • 幽灵峰的识别：该工具能够识别数据中的幽灵峰（phantom peaks），这些通常是由于非特异性结合或实验操作错误产生的假阳性信号。
  • 峰值检测：Phantompeakqualtools还包含峰值检测功能，能够识别和分析数据中的信号峰，这对于确定蛋白质-DNA相互作用和其他基因组标记非常重要。

Github:

  • https://github.com/kundajelab/phantompeakqualtools

编程语言：R

2来源

题目：ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia
期刊：Genome Research
日期：2012/9DOI：10.1101/gr.136184.111

3如何安装

conda 安装

推荐使用conda安装，可以直接安装运行所需的依赖包，方便快捷

conda activate chipseq
conda install bioconda::phantompeakqualtools
依赖

4使用示例

计算cross correlation

##激活环境
conda activate chipseq

run_spp.R -rf \
-c=/home/data/t020559/chip_seq/GSE205035_PRJNA843319/e_mkdup/SRR19436494_H3K27ac-OXA_mkdup.bam \
-p=8 \
-odir=/home/data/t020559/chip_seq/GSE205035_PRJNA843319/test \
-savp=SRR19436494_H3K27ac-OXA_mkdup.pdf \
-out=SRR19436494_H3K27ac-OXA_mkdup_cross.txt

##参数释义
-c=<ChIP_alignFile>：#ChIP测序数据文件的完整路径和名称，
-i=<Input_alignFile>：#输入对照（input control）文件的完整路径和名称
-p=<num> ：#设置线程，默认为0
-s=<min>:<step>:<max>：#评估交叉相关性的链移位范围，默认为-500:5:1500
-rf : #如果存在同名文件，则替换（覆盖输出结果）
-tmpdir=<tempdir>：#临时目录，如果未指定，则使用R函数 tempdir() 的结果
-filtchr=<chrnamePattern>：#用于移除映射到特定染色体的标签的模式，例如使用"chrM"将移除所有名称中包含"chrM"的染色体的标签
-odir=<outputDirectory>：#输出目录的名称，如果未设置，则使用当前目录
-savn=<narrowpeakfilename>：#保存NarrowPeak文件的名称
-savr=<regionpeakfilename>：#保存RegionPeak文件的名称
-savd=<rdatafile> ：#保存R数据文件
-savp=<plotdatafile> ：#保存交叉相关图的文件
-out=<resultfile>：#将peakshift/phantomPeak结果追加到文件中
运行示例

5结果说明

cat SRR19436494_H3K27ac-OXA_mkdup_cross.txt
SRR19436494_H3K27ac-OXA_mkdup.bam       26755432        0,330,375       0.252146613583201,0.237250183628446,0.237133227855954   150     0.3129999       1500    0.2282759       1.10457   0.2817472       -1

结果共11列，分别为：

  • Filename - 文件名，包含tagAlign或BAM格式的文件名。
  • numReads - 有效测序深度，即输入文件中映射读取的总数。
  • estFragLen - 逗号分隔的链交叉相关峰值，按相关度降序排列。
  • corr_estFragLen - 逗号分隔的链交叉相关值，也是按降序排列（与第3列的顺序相对应）。
  • phantomPeak - 读长/幽灵峰的链位移。
  • corr_phantomPeak - 在幽灵峰处的相关值。
  • argmin_corr - 交叉相关性最低时的链位移。
  • min_corr - 交叉相关性的最小值。
  • NSC - 标准化链交叉相关系数，计算公式为第4列值除以第8列值（NSC = COL4 / COL8）。
  • RSC - 相对链交叉相关系数，计算公式为（第4列值减第8列值）除以（第6列值减第8列值）（RSC = (COL4 - COL8) / (COL6 - COL8)）。
  • QualityTag - 基于阈值化的RSC的质量标记（编码：-2：非常低，-1：低，0：中等，1：高，2：非常高）。

图中的黑线代表实验数据的交叉相关性，蓝色虚线代表片段长度峰，红色虚线表示主峰的真实位移。这些峰的位置说明了测序片段和可能的非特异性信号的分布。NSC值为1.10457，这个值大于1，表明相对于背景信号，片段长度峰有较好的信号强度。RSC值为0.28147，这个值低于常规的阈值1，说明信号与幽灵峰之间的对比度不是很高，信号与噪声的比例较低。

NSC（Normalized Strand Coefficient） ：片段长度交叉相关峰与背景交叉相关的比值，是评估信号与噪声比的一个强有力的度量。

  • 数值范围：NSC的值从1开始，可以达到更高的正数。数值1表示数据中信号与噪声的比例是平衡的，超过1表示信号强于噪声。
  • 临界阈值：1.1被认为是一个关键的临界点。NSC值大于或等于1.1通常表示数据具有相对较高的信号强度。
  • 低于阈值的含义：当NSC值远低于1.1（如小于1.05）时，可能表明数据中的信号与噪声比低，或峰值数量少。这种情况可能是生物学上合理的，例如某些转录因子在特定组织类型中仅绑定少数位点；也可能是由于样本质量不佳。

RSC（Relative Strand Coefficient）：片段长度峰与读长峰的比值，同样用于评估信号的质量

  • 数值范围：RSC的值从0开始，可以增加到更高的正数。数值1是另一个关键的临界点。
  • 临界阈值：RSC值在1或以上通常表示样本中的信号较为清晰和显著。
  • 低于阈值的含义：当RSC值显著低于1（如小于0.8）时，同样可能表示数据中的信号与噪声比低。低RSC分数可能是由于ChIP实验失败、测序质量低、读取映射错误多、测序深度不足等原因造成的。同NSC一样，如果一个数据集中的结合位点数量很少（如少于200），并且这种现象在生物学上是合理的，也可能显示低RSC分数。

6Peak calling

其也有峰值检查的功能，示例代码如下。不过速度很慢。不是很建议使用。推荐使用：MACS 。见MACS3—探索基因组调控的钥匙

run_spp.R -rf \
-c=/home/data/t020559/chip_seq/GSE205035_PRJNA843319/e_mkdup/SRR19436494_H3K27ac-OXA_mkdup.bam \
-i=/home/data/t020559/chip_seq/GSE205035_PRJNA843319/e_mkdup/SRR19436498_Input-OXA_mkdup.bam \
-fdr=0.01 \
-odir=/home/data/t020559/chip_seq/GSE205035_PRJNA843319/test/OXA \
-savn -savr -savp -savd
运行示例



文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：