HIC测序数据生信分析——第二节，预处理HIC数据

HIC原始数据——HICclean数据

这一步是所有二代数据都要做的，使用软件Trim_galore和fastqc，修剪低质量的碱基和接头序列，过滤低质量的Reads。
可以使用seqkit软件做fastq文件的统计，很方便快速。
结果文件：Unknown_BM942-001H0001_1_val_1.fq和Unknown_BM942-001H0001_2_val_2.fq
seqkit stat -j 36 -T Unknown_BM942-001H0001_1_val_1.fq
seqkit fx2tab -n -l genome.fasta > len.txt

HICclean数据——预处理后的HIC数据

这里预处理的意思是：对HIC数据中嵌合体序列进行修剪，与contig水平基因组比对（通常单端比对，再合并），去除噪声序列（再连（Re-liagation），自连（Same Circularised）等），以及去除因PCR重复的序列等。
这里提供三种方法。方法一：Hicup软件；方法二：ALLHiC软件；方法三：juicer软件。前两者都会拿到bam文件，而后者结果是拿到一个仅用于3D-DNA的merge_nodups.txt文件

方法一：Hicup软件

有四步，分别是hicup_truncater、hicup_mapper、hicup_filter和hicup_deduplicator

安装Hicup软件

官网：https://github.com/StevenWingett/HiCUP/
#可以直接使用conda安装，先conda search hicup，查看有哪些版本可以进行安装。下载最新版conda install hicup=0.9.2
conda create -n hicup
conda activate hicup
conda install hicup=0.9.2
conda deactivate

准备工作

准备文件：contig的基因组序列文件genome.fa和HICcleandata：Unknown_BM942-001H0001_1_val_1.fq和Unknown_BM942-001H0001_2_val_2.fq

#准备一：构建基因组的索引文件
conda activate hicup
cd /home/zhaohuiyao/Genome_hic/baimaike/03hicup
ln -s …/Genome/genome.fa ./
bowtie2-build genome.fa genome.fa

#准备二：参考基因组的酶切位点图谱
hicup_digester --re1 ^GATC,DpnII --genome O_sis_digester_db genome.fa
#结果文件Digest_O_sis_digester_db_DpnII_None_18-03-19_31-08-2023.txt

#准备三：hicup的配置文件
hicup --example #生成一个配置文件的模板文件hicup_example.conf
#对模板配置文件进行修改
#修改1：Outdir: /home/zhaohuiyao/Genome_hic/baimaike/03hicup。指定输出目录
#修改2：Threads: 36。线程数
#修改3：Bowtie2: /home/zhaohuiyao/miniconda3/envs/hicup/bin/bowtie2。指定Bowtie2软件位置
#修改4：Index: /home/zhaohuiyao/Genome_hic/baimaike/03hicup/genome.fa。指定参考开基因组的索引文件
#修改5：
Digest: /home/zhaohuiyao/Genome_hic/baimaike/03hicup/Digest_O_sis_digester_db_DpnII_None_16-26-48_11-08-2023.txt。指定参考基因组的酶切位点图谱文件
#修改6：Format: Sanger。指定fastq格式（从这四个中选择，‘Sanger’, ‘Solexa_Illumina_1.0’, ‘Illumina_1.3’，‘Illumina_1.5’）。如果没有指定，HiCUP会自己进行分析。当然你可以在fastqc结果中查看到这个信息
#修改7：/home/zhaohuiyao/Genome_hic/baimaike/01Cleandata/Unknown_BM942-001H0001_1_val_1.fq和/home/zhaohuiyao/Genome_hic/baimaike/01Cleandata/Unknown_BM942-001H0001_2_val_2.fq
#还可以修改的位置有Quiet:0、Keep:0、Zip:1、Longest: 700、Shortest:

运行hicup

hicup --config hicup_example.conf
#生成一个html文件，包括四个模块的结果的结果以及可视化。Unknown_BM942-001H0001_1_val_1.Unknown_BM942-001H0001_2_val_2.hicup.bam.HiCUP_summary_report_fGOPpcmbtt_18-07-56_31-08-2023.html

#高质量的HIC数据文件：Unknown_BM942-001H0001_1_val_1.Unknown_BM942-001H0001_2_val_2.hicup.bam
samtools view Unknown_BM942-001H0001_1_val_1.Unknown_BM942-001H0001_2_val_2.hicup.bam | wc
conda deactivate

方法二：ALLHiC软件

安装ALLHiC软件

#使用安装包安装
#安装依赖软件samtools v1.9+、bedtools、matplotlib v2.0+

cd /home/zhaohuiyao/Biosoft
git clone https://github.com/tangerzhang/ALLHiC
cd ALLHiC/
chmod +x bin/*
chmod +x scripts/* 

#使用Docker安装
docker search ALLHiC
docker pull wangnan9394/allhic

准备工作

准备文件：contig的基因组序列文件genome.fa和HICcleandata：Unknown_BM942-001H0001_1_val_1.fq和Unknown_BM942-001H0001_2_val_2.fq

运行ALLHiC

docker run -it -d -v /home/zhaohuiyao/Genome_hic:/home/zhaohuiyao/Genome_hic --name allhic001 wangnan9394/allhic /bin/bash
docker exec -it allhic001 /bin/bash（进入容器，使用exit退出，使用docker kill allhic001 干掉容器）
export PATH=/home/ALLHiC-master/scripts/:/home/ALLHiC-master/bin/:$PATH

#准备一：构建基因组的索引文件
cd /home/zhaohuiyao/Genome_hic/baimaike/04ALLHiC
ln -s …/Genome/genome.fa ./
bwa index genome.fa genome.fa
samtools faidx genome.fa

#将HIC数据比对到基因组上
bwa aln -t 36 genome.fa …/01Cleandata/Unknown_BM942-001H0001_1_val_1.fq.gz > hic_clean_R1.sai
bwa aln -t 36 genome.fa …/01Cleandata/Unknown_BM942-001H0001_2_val_2.fq.gz > hic_clean_R2.sai
bwa sampe genome.fa hic_clean_R1.sai hic_clean_R2.sai …/01Cleandata/Unknown_BM942-001H0001_1_val_1.fq.gz …/01Cleandata/Unknown_BM942-001H0001_2_val_2.fq.gz > hic_clean.bwa_aln.sam

#过滤
PreprocessSAMs.pl hic_clean.bwa_aln.sam genome.fa MBOI
#MBOI和DpnII识别序列一致
filterBAM_forHiC.pl hic_clean.bwa_aln.REduced.paired_only.bam hic_clean.sam
samtools view -bt genome.fa.fai hic_clean.sam > hic_clean.bam
#准备工作完成，拿到预处理文件hic_clean.bam

方法三：juicer软件

安装juicer软件

#使用Docker安装
docker search juicer
docker pull rnakato/juicer

准备工作

准备文件：contig的基因组序列文件genome.fa和HICcleandata：Unknown_BM942-001H0001_1_val_1.fq和Unknown_BM942-001H0001_2_val_2.fq

运行juicer

#创建容器
docker run -it -d -v /home/zhaohuiyao/Genome_hic:/home/zhaohuiyao/Genome_hic --name juicer001 rnakato/juicer /bin/bash
docker exec -it juicer001 /bin/bash（进入容器，使用exit退出，使用docker kill juicer001 干掉容器）

#准备一：构建基因组的索引文件
cd /home/zhaohuiyao/Genome_hic/baimaike/05juicer/Genome
ln -s …/…/Genome/genome.fa ./
baw index genome.fa genome.fa
samtools faidx genome.fa

#创造酶切位点文件
python /opt/juicer/misc/generate_site_positions.py DpnII genome genome.fa
#创建contig长度文件
awk ‘BEGIN{OFS=“\t”}{print $1, $NF}’ genome_DpnII.txt > genome.chrom.size

#准备二：HIC数据
cd /home/zhaohuiyao/Genome_hic/baimaike/05juicer/fastq
ln -s …/…/01Cleandata/Unknown_BM942-001H0001_1_val_1.fq.gz hic_clean_R1.fastq.gz
ln -s …/…/01Cleandata/Unknown_BM942-001H0001_2_val_2.fq.gz hic_clean_R2.fastq.gz

#运行juicer
cd /home/zhaohuiyao/Genome_hic/baimaike/05juicer/
/opt/juicer/scripts/juicer.sh -g genome -s DpnII -z …/Genome/genome.fa -y …/Genome/genome_DpnII.txt -p …/Genome/genome.chrom.size -D /opt/juicer -t 36
#结果文件aligned/merged_nodups.txt
#juicer在生成merged_nodups.txt文件后，还会跑inter_30之类文件，是后续3D-DNA不需要的，如果这里报错，则可以选择忽略。

以上三种方法，选择一种，完成HIC数据预处理工作。