本文详细介绍了在Linux环境下处理ChIA-PET数据的过程，包括构建bwa基因组索引文件、下载和解压数据、运行ChIA-PET命令以及解决遇到的报错问题。在构建索引时遇到问题，通过合并基因组文件并建立bwa索引来解决。在运行ChIA-PET时，针对找不到命令的错误，采用修改环境变量的方法进行修复。文章还提及了peak calling的概念及其在基因组数据分析中的作用。
文章目录
ChIA-PET2 -g hg19.fa -b human.hg19.genome -d 1 -f ENCFF000KYG.fastq -r ENCFF000KYK.fastq -o OUTdir12 -n index-1
一、构建bwa的基因组索引文件
1、一步到位下载hg19基因组文件
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19*
放入服务器
bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>#官网提供的命令
bwa index -a bwtsw hg19.fa#实际操作
最后生成文件：hg19.fa.amb、hg19.fa.ann、hg19.fa.bwt、hg19.fa.pac和hg19.fa.sa。
amb是ambiguous的缩写，也就是模棱两可的意思，也就是除了ATCG/atcg以外的字符. amb和ann用来记录基因组中除了ATCG以外碱基的信息。而pac文件则是碱基信息高度压缩。
构建索引时需要注意的问题：bwa构建索引有两种算法
-a bwtsw对于短的参考序列是不工作的，必须要大于等于10Mb
-a is是默认参数，这个参数不适用于大的参考序列，必须要小于等于2G
2、遇到报错
(base) [root@node01 chia-pet1]# ChIA-PET2 -g ucsc.hg19.fasta.pac -b human.hg19.genome -f
[05-21 11:40:30] Start ChIA-PET2 V0.9.3 ...
[05-21 11:40:30] Running Step 1: Trim Linker ...
trimLinker -t 1 -m 0 -k 0 -e 0 -l 15 -o OUTdir -n index-1 -A GTTGGATAAG -B GTTGGAATGT ENC
thread is 1
mode is 0
keepempty is 0
Error allowed is 0
min length of trimmed read is 15
Output dir is OUTdir
Output name is index-1
linkerA is GTTGGATAAG
linkerB is GTTGGAATGT
Reads 1: ENCFF000KYB.fastq.gz
Reads 2: ENCFF000KYQ.fastq.gz
Processed 2000000 pair reads
（省略一部分）
1Empty PETs: 38719182
2Empty PETs: 3761387
Valid PETs: 142457250
[05-21 12:14:00] Running Step 2: BWA ...
bwa_wrap ucsc.hg19.fasta.pac OUTdir/index-1_1.valid.fastq 1 OUTdir/index-1_1.valid.sam 0
Running BWA on trimmed reads ...
bwa mem -t 1 ucsc.hg19.fasta.pac OUTdir/index-1_1.valid.fastq | samtools view -h -F 2048
[E::bwa_idx_load_from_disk] fail to locate the index files
Exit: Please check input or Rerun step 2
(base) [root@node01 chia-pet1]# bwa
Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.17-r1188
Contact: Heng Li <lh3@sanger.ac.uk>
3、解决：
①将24个参考基因组.fa文件进行合并，手动合称为一个
cat *.fa > hg19.fa
得到hg19.fa文件
②hg19.fa文件建bwa索引
bwa index -a bwtsw hg19.fa
得到这五个文件
二、找hg19的genome文件
在bedtools中，genomes文件夹下
三、在ENCODE下载K562细胞系中，CTCF的ChIA-PET数据
1、下载数据
2、解压数据
由于解压后文件二十多G，于是为了试验程序，又下载了一个较小的数据
得到这样两个文件
解压
这是PORL2A解压后的数据
四、跑ChIA-PET
1、命令
ChIA-PET2 -g hg19.fa -b human.hg19.genome -f ENCFF000KYB.fastq -r ENCFF000KYQ.fastq -o OUTdir -n index-1
-g bwa的基因组索引文件
-b 为bedtools提供的染色体大小文件（UCSC上下载）
-f，-r 输入的两个fastq（.gz）文件
-A，-B 两个linker序列，默认为GTTGGATAAG 和 GTTGGAATGT
-o 输出目录，默认为output
-n 输出文件的前缀名
这里是引用
$ ChIA-PET2 -h
usage : ChIA-PET2 -g genomeindex -b bedtoolsgenome -f fq1 -r fq2 -A linkerA -B linkerB -o OUTdir -n prefixname
Use option -h|--help for more information
ChIA-PET2 0.9.3 2017.11.07
----------------------------
OPTIONS
-s|--start: start from which step(1:8): 1:Trim Linkers; 2:Map Reads; 3:Build PETs; 4:Call Peaks; 5:Find Interactions
6:Plot QC; 7:Estimate statistical confidence; 8:Phase PETs(optional), default=1
-g|--genome: genome index for bwa
-b|--bedtoolsgenome: chromsomes size file for bedtools
-f|--forward: one fastq(.gz) file
-r|--reverse: the other fastq(.gz) file
-A|--linkerA: one linker sequence, default=GTTGGATAAG
-B|--linkerB: the other linker sequence, default=GTTGGAATGT
-o|--output: output folder, default=output
-n|--name: output prefix name, default=out
-m|--mode: 0,1,2; 0: A/B linkers; 1: bridge liker; 2: Enzyme site, default=0
-e|--err: Maximum mismatches allowed in linker sequence, default=0
-k|--keepempty: 0,1,2; 0:No linker-empty reads; 1:keep 1 linker-empty read; 2:keep 2 linker-empty reads. default=0
-t|--thread: threads to run, default=1
-d|--short: short reads (0 or 1), default=0 for reads >70bp. If the read length is about 20bp, set d=1
-M|--macs2 parameters, default="-q 0.05"
-Q|--mapq: mapq cutoff, default=30
-C|--cutoffPET: PET count cutoff before running MICC, default=2
-S|--slop: slop length, default=100
-E|--extend: extend length on both sides, default=500
-l|--length: min length of reads after linker trimming. default=15
-P|--phased: optional phased genotype file: 'chr1\tstart\tend\tA\tC'
[-h|--help]: help
[-v|--version]: version
2、报错
找不到命令
3、解决
①三种查找路径的方法
which ChIA-PET2
whereis ChIA-PET2
find / -iname "ChIA-PET2" #运行以下命令以在系统中查找给定文件
②修改环境变量
将找到的路径添加到环境变量里
vim /etc/profile
#把下面的命令加入末尾
export PATH=$PATH:/root/bin/ChIA-PET2_0.9.3/bin
注意，一定在后缀后加/bin
③修改变量之后，需要马上生效变量设置。
source /etc/profile
4、把所需数据都放到一个文件夹下
通过cd，到指定目录下，输入命令开始跑数据
第一步大概半小时
第二步大概
五、报错
1、请重复第三步，或检查输入数据
[main] CMD: bwa mem -t 1 hg19.fa OUTdir9/index-1_2.valid.fastq
[main] Real time: 4573.145 sec; CPU: 4626.485 sec
[05-27 18:06:30] Running Step 3: Build Bedpe ...
buildBedpe OUTdir9/index-1_1.valid.sam OUTdir9/index-1_2.valid.sam OUTdir9/index-1.bedpe 30 1 1
Running buildbedpe...
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 1353876 PET
All 33353876
Reads1 Low MAPQ 0 0%
Reads1 Unmapped 33353876 100%
Reads2 Low MAPQ 0 0%
Reads2 Unmapped 33353876 100%
Output PETs 0 0%
Exit: Please check input or Rerun step 3
(base) [root@node01 chia-pet1]#
1、1处理
ChIA-PET2 -g hg19.fa -b human.hg19.genome -d 1 -f ENCFF000KYG.fastq -r ENCFF000KYK.fastq -o OUTdir12 -n index-1
2、masc2 peak calling报错
[05-29 02:43:13] Running Step 3: Build Bedpe ...
buildBedpe OUTdir12/index-1_1.valid.sam OUTdir12/index-1_2.valid.sam OUTdir12/index-1.bedpe 30 1 1
Running buildbedpe...
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 2000000 PET
Processed 1353876 PET
All 33353876
Reads1 Low MAPQ 22004694 65.97%
Reads1 Unmapped 8067789 24.19%
Reads2 Low MAPQ 22752786 68.22%
Reads2 Unmapped 7406492 22.21%
Output PETs 1980267 5.937%
removeDup OUTdir12/index-1.bedpe OUTdir12/index-1.rmdup.bedpe 1
Running removeDups...
Duplicates 360379
Uniques 1619888
Intra 433597
Inter 29182
OneEndMapped 1157109
buildTagAlign OUTdir12/index-1.rmdup.bedpe OUTdir12/index-1.rmdup.bedpe.tag
Running buildTagAlign...
[05-29 02:49:29] Running Step 4: MACS2 ...
macs2_wrap OUTdir12/index-1.rmdup.bedpe.tag OUTdir12/index-1 "-q 0.05"
Running Macs2...
macs2 callpeak -t OUTdir12/index-1.rmdup.bedpe.tag -f BED --keep-dup all -n OUTdir12/index-1 -q 0.05
Traceback (most recent call last):
File "/home/fengyongxian/shm/tartar/miniconda/yes/bin/macs2", line 4, in <module>
__import__('pkg_resources').run_script('MACS2==2.2.7.1', 'macs2')
File "/home/fengyongxian/shm/tartar/miniconda/yes/lib/python3.8/site-packages/pkg_resources/__init__.py", line 3243, in <module>
def _initialize_master_working_set():
File "/home/fengyongxian/shm/tartar/miniconda/yes/lib/python3.8/site-packages/pkg_resources/__init__.py", line 3226, in _call_aside
f(*args, **kwargs)
File "/home/fengyongxian/shm/tartar/miniconda/yes/lib/python3.8/site-packages/pkg_resources/__init__.py", line 3255, in _initialize_master_working_set
working_set = WorkingSet._build_master()
File "/home/fengyongxian/shm/tartar/miniconda/yes/lib/python3.8/site-packages/pkg_resources/__init__.py", line 568, in _build_master
ws.require(__requires__)
File "/home/fengyongxian/shm/tartar/miniconda/yes/lib/python3.8/site-packages/pkg_resources/__init__.py", line 886, in require
needed = self.resolve(parse_requirements(requirements))
File "/home/fengyongxian/shm/tartar/miniconda/yes/lib/python3.8/site-packages/pkg_resources/__init__.py", line 772, in resolve
raise DistributionNotFound(req, requirers)
pkg_resources.DistributionNotFound: The 'MACS2==2.2.7.1' distribution was not found and is required by the application
ERROR in running Macs2 ...
(base) [root@node01 chia-pet1]#
peak calling
从功能上看，可以把peak calling看成一种对高维基因组数据的数据降维（dimensionality reduction）或特征选取（feature selection）
peak calling是找被较多reads所覆盖的区域，找到基因组上开放程度高的区域，这些区域一般开放信号较强，更有可能有生物学意义（比如转录因子结合）
扫一扫
6401
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
