chip-seq生信笔记

Posted 2023-04-27

参考技术A 染色质免疫共沉淀技术（Chromatin Immunoprecipitation，ChIP）也称结合位点分析法，是研究体内蛋白质与DNA相互作用的有力工具，通常用于转录因子结合位点或组蛋白特异性修饰位点的研究。

知乎文章： https://zhuanlan.zhihu.com/p/90180058
流程文章： https://www.jianshu.com/p/21e8c51fca23
从数据到igv可视化分析： https://blog.csdn.net/qq_29300341/article/details/54811085

使用fastqc软件

fastqc file1 file2

使用multiqc软件进行多个qc结果的合并

multiqc <analysis directory>

bowtie主要适用于将短序列比对到参考genome上，速度快。

mapping序列到genome上，首先要建立genome的index。command需要待建立的genome文件，和输出index的文件夹。

bowtie2-build [options]* <reference_in> <bt2_index_base>

bowtie将reads对比到genome上，生成sam文件；sam文件是序列比对到基因组上的结果展示，或者展示多重比对结果。sam文件包括比对的注释信息（header section）和比对结果部分（alignment section）。注释信息是比对操作的说明，包括参考序列，程序说明等；比对结果事对每一个片段（segment）的说明，包括比对到参考序列的位置，mapping的质量等的说明。

bowtie2 [options]* -x <bt2-idx> -1 <m1> -2 <m2> | -U <r> | --interleaved <i> | -b <bam> [-S <sam>]

bowtie2对比序列，首先要看是单端测序还是双端测序，双端测序需要将两端测序的链匹配起来。还可以选定输出文件的后缀和存放的目录。

使用samtools将sam文件转换为bam文件（view），将bam序列进行排序（sort）。使用samtools对bam文件的对比结果，并输出统计结果，检验测序read的质量。

samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]

options -b 输出bam文件

samtools sort [options...] [in.bam]

对bam文件进行sort排序能进一步减小文件体积，加快运算速度。

chip-seq的可以检测的富集方式包括两种：1.broad domains和narrow peak。broad domain是组蛋白在整个基因组的修饰，narrow peaks是特定的突出指，如转录因子的结合。当需要对特定的目标靶点进行研究时，可以利用treat组和control组进行对比，找出二者的不同。所谓callpeak，是指寻找基因组上的表达峰peak，chip-seq是对蛋白结合的DNA进行测序，每个read都意味着有一个蛋白结合到基因组的该处上，基因组的peak就是read表达量最高的地方，调控因子一般都在gene的上游或者下游，离gene越近的调控因子与gene表达的相关性越高，所以要callpeak，寻找不同gene之间以及gene和转录因子之间的关系等。

macs2 callpeak [-h] -t TFILE [TFILE ...] [-c [CFILE [CFILE ...]]]

-t 为treatment组数据， -c 为control组数据， -g 选择genome，还可以设置输出的目录和文件名， --bdg 可以bdg文件，用于igv查看peaks。

macs检验值设立： https://www.jianshu.com/p/390f6d57488d

将callpeak生成的.bdg文件直接放入IGV(intergrative genomic viewer)。

motif分析。寻找peak序列的共同模式序列。motif的输入文件为call的".bed"文件。“.bed”文件包含的信息为peak（summit）所在染色体和具体位置。

annotatePeaks.pl <peak file | tss> <genome version> [additional options...]
分为两种，一种只对peak的信息进行注释，展现peak相关的gene和到geneTSS
（transcription start site）的距离。

另一种是对peak和read的信息进行annotate，显示reads在peak summit 两侧的分布情况。
首先要用homer的 makeTagDirectory 对.bam文件进行处理生成tag文件夹。生成tag文件的处理，包括对bam文件的排序，质控，以及生成一些后续分析需要的重要参数。
然后再进行annotatePeaks.pl分析，加上参数
-d <tag directory 1> [tag directory 2] ... (list of experiment directories to show tag counts for)
再使用生成的annotation的文件，使用R进行绘图。

deeptools分析画图。

只需要macs得出的peaks的bed文件，和选择参考gene组就可以。
findMotifsGenome.pl <pos file> <genome> <output directory> [additional options]

Chip-seq peak annontation

Chip-seq peak annontation

Chip-seq peak annontation

PeRl

narrowPeak/boardPeak

narrowPeak/boardPeak 是ENCODE可提供下载的两种 Chip-seq 经过参考人类基因组mapping后的关于peak的数据.
其他类型的seq数据储存个数可以参看FAQformat

narrowPeak

数据按照以下规则储存:

1. string chrom: "Reference sequence chromosome or scaffold"
2. uint   chromStart: "Start position in chromosome"
3. uint   chromEnd: "End position in chromosome"
4. string name: "Name given to a region (preferably unique). Use . if no name is assigned"
5. uint   score: "Indicates how dark the peak will be displayed in the browser (0-1000) "
6. char[1]  strand: "+ or - or . for unknown"
7. float  signalValue: "Measurement of average enrichment for the region"
8. float  pValue: "Statistical significance of signal value (-log10). Set to -1 if not used."
9. float  qValue: "Statistical significance with multiple-test correction applied (FDR -log10). Set to -1 if n-ot used."
10. int   peak: "Point-source called for this peak; 0-based offset from chromStart. Set to -1 if no point-sour-ce called."

boardPeak

数据按照以下规则储存:

1. string chrom: "Reference sequence chromosome or scaffold"
2. uint   chromStart: "Start position in chromosome"
3. uint   chromEnd: "End position in chromosome"
4. string name: "Name given to a region (preferably unique). Use . if no name is assigned"
5. uint   score: "Indicates how dark the peak will be displayed in the browser (0-1000) "
6. char[1]  strand: "+ or - or . for unknown"
7. float  signalValue: "Measurement of average enrichment for the region"
8. float  pValue: "Statistical significance of signal value (-log10). Set to -1 if<BR> not used."
9. float  qValue: "Statistical significance with multiple-test correction applied (FDR -log10). Set to -1 if n-ot used."

在接下去的peak annotation中,只演示narrowPeak.bed格式数据.

示例数据

演示数据来自ENCODE 的 H3K9me3 的Chip-seq,样本ID为ENCFF199BLM.

样本的基本信息:

• Homo sapiens liver male adult (32 years)
  • Target: H3K9me3
  • Lab: Bing Ren, UCSD
  • Project: Roadmap

narrowPeak 数据基本信息:

##   seqnames     start       end       name score strand signalValue  pValue
## 1    chr10 100134639 100134831  Peak_7974   178      .     4.41261 6.68330
## 2    chr10 100446376 100446664  Peak_4189   244      .     5.13248 8.41215
## 3    chr10 100779568 100779699 Peak_32369   101      .     3.34436 4.50936
## 4    chr10  10088147  10088346 Peak_28509   112      .     4.05830 4.84890
## 5    chr10 101149252 101149594  Peak_6146   211      .     4.89252 7.53305
## 6    chr10 101173156 101173424 Peak_32985    94      .     3.45278 4.33257
##    qValue peak
## 1 2.96490  132
## 2 4.37217  165
## 3 1.47374  105
## 4 1.69566   90
## 5 3.67439  174
## 6 1.30993  237

构建注释文件

在peak annotation中需要一个用于参考的注释文件,需要包括一下信息: 1. 染色体名 2. 转录本起始位点 3. 转录本终止位点 4. gene的名字 5. 正反链

注释文件可以直接从外部导入,也可以利用biomaRt 生成

library(ChIPpeakAnno)
library(biomaRt)
mart <- useMart("ensembl")
datasets <- listDatasets(mart)
mart <- useDataset("hsapiens_gene_ensembl",mart)
# 需要筛选的特征
props <- c("ensembl_gene_id", "external_gene_name", "transcript_biotype", "chromosome_name", "start_position", "end_position", "strand")

# 筛选染色体号
lincRNA <- subset(
  getBM(attributes=props, mart=mart, filters = "chromosome_name", values = c(1:22,"X","Y")),
  transcript_biotype == "lincRNA"
  )
# 在染色体前加"chr"保持和narrowPeak数据一致
lincRNA[,4] <- paste0("chr", lincRNA[,4])

得到的数据包含以下信息:

##   ensembl_gene_id external_gene_name transcript_biotype chromosome_name
## 1 ENSG00000276255       RP5-881P19.7            lincRNA            chr1
## 2 ENSG00000234277          LINC01641            lincRNA            chr1
## 3 ENSG00000238107      RP11-495P10.5            lincRNA            chr1
## 4 ENSG00000274020          LINC01138            lincRNA            chr1
## 5 ENSG00000225620      RP11-569A11.2            lincRNA            chr1
## 6 ENSG00000237520       RP11-443B7.2            lincRNA            chr1
##   start_position end_position strand
## 1      228073909    228076550     -1
## 2      227393591    227431035      1
## 3      148295180    148297556      1
## 4      148290889    148519604     -1
## 5      202632428    202632911      1
## 6      234957231    234959989      1

进行peak annotation

有了这个注释文件以后,我们就可以根据我们想要的筛选规则对peak进行注释,主要用到的包是 ChIPpeakAnno.用到的函数为annotatePeakInBatch.

# 将前面得到的注释文件转换为RangedData对象
library(ChIPpeakAnno)
myCustomAnno <- RangedData(
  IRanges(
    start=lincRNA[,"start_position"], 
    end=lincRNA[,"end_position"],
    names=lincRNA[,"ensembl_gene_id"]),
  space=lincRNA[,"chromosome_name"],
  strand=lincRNA[,"strand"])
# 读入需要注释的narrowPeak数据
bed_file <- read.table("ENCFF199BLM.bed", header = T, sep = "\t", stringsAsFactors = F)
#　将peak数据转换为GRanges
peaks <- toGRanges(bed_file, format="narrowPeak", colNames = colnames(bed_file))
#　根据需要进行筛选
anno <- annotatePeakInBatch(peaks, AnnotationData=myCustomAnno, 
                            output="overlapping", 
                            FeatureLocForDistance="TSS",
                            bindingRegion=c(-2000, 2000))

从anno中提取我们需要的lincRNA的ID

result_lincRNA <- [email protected]$feature
head(result_lincRNA)

## [1] "ENSG00000237579" "ENSG00000235180" "ENSG00000232259" "ENSG00000204365"
## [5] "ENSG00000226578" "ENSG00000260137"