ChIP-seq数据比对实战

豆豆写于2020.5.7
首先是ChIP-seq分析的前言介绍部分：
1：了解ChIP-seq的实验流程
2：继续了解ChIP-seq
3：关于ChIP-seq的实验对照与偏差来源
4：ChIP-seq的实验设计补充
5：ChIP-seq数据库及实战数据介绍
然后开始实战部分：
6：ChIP-seq计算资源准备与实战数据下载
7：ChIP-seq数据质控和过滤
8：ChIP-seq数据比对注意事项
这一次将介绍bowtie2比对的方法和结果检查

1 利用bowtie2比对

bowtie软件开发的初衷就是：基于Burrows-Wheeler算法，针对大基因组的物种，进行短序列的比对工作。相比之前的版本，bowtie2不再设定错配、多重比对的阈值参数，而采用了一种打分机制，分值越高说明比对越相似。bowtie和bowtie2，分别对处理50bp以下和50bp以上的数据。

比对的步骤一般是：构建参考基因组索引 + 将fq文件比对到索引

比对后的结果还要根据比对质量值进行过滤，MAPQ>10 一般就可以保留那些唯一比对的reads

1.1 首先下载参考基因组

WKD=~/public/mm_nrf1
CONFIG=$WKD/sra/0-config.txt
CLEAN_DIR=$WKD/clean
REF_DIR=$WKD/reference
ALIGN_DIR=$WKD/align
wget http://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz -O $REF_DIR/chromFa.tar.gz
tar -zxvf $REF_DIR/chromFa.tar.gz -C $REF_DIR
rm $REF_DIR/chromFa.tar.gz $REF_DIR/*random* $REF_DIR/chrUn*

1.2 构建索引

for i in 1 $(seq 11 19) $(seq 2 9) M X Y;
do 
    #echo $i
    cat $REF_DIR/chr${i}.fa >> $REF_DIR/mm10.fasta
done
# 10线程(约36min)
bowtie2-build --threads 10 $REF_DIR/mm10.fasta $REF_DIR/mm10

1.3 比对--先看一步一步的拆分运行

# 首先是最原始的单端（-U）比对=》fq to sam
bowtie2 -x $REF_DIR/mm10 -U $CLEAN_DIR/${name}.fastq.gz > ${ALIGN_DIR}/${name}.sam 
# 然后sam过滤(-q 10)并转为bam
samtools view -Sb -q 10 ${ALIGN_DIR}/${name}.sam > ${ALIGN_DIR}/${name}.nonSorted.bam
# 将bam进行排序（sort）
samtools sort ${ALIGN_DIR}/${name}.nonSorted.bam > ${ALIGN_DIR}/${name}.sorted.bam
# 删除中间文件，节省空间
rm ${ALIGN_DIR}/${name}.sam ${ALIGN_DIR}/${name}.nonSorted.bam
# 构建index
samtools index ${ALIGN_DIR}/${name}.sorted.bam

1.4 比对--一步运行（使用管道）

bowtie2 -x $REF_DIR/mm10 -U $CLEAN_DIR/${name}.fastq.gz  | samtools view -Sb -q 10 | samtools sort > ${ALIGN_DIR}/${name}.sorted.bam

1.5 可视化一下结果

# 看bam结果，例如
samtools view align/H3K27AC_CHIP_TKO_1.sorted.bam | head -1
# H3K27AC_CHIP_TKO_1.26986706    0   chr1    3000563 42  50M *   0   0   TAGATTTGCTGTCAGGCTGCTAGTGTATACTCTAGTTTCCTTTTGGAGGC  GAAAGGGGIIGIIGGGGGGIIIIIIIIIGIIIIIIIIIIIGIIIGGIIIG  AS:i:0  XN:i:0XM:i:0    XO:i:0  XG:i:0  NM:i:0  MD:Z:50 YT:Z:UU

# 看bed结果，例如
bamToBed -i  align/H3K27AC_CHIP_TKO_1.sorted.bam | head -1
# chr1    3000562 3000612 H3K27AC_CHIP_TKO_1.26986706 42  +

再看一下坐标系统

• 1-based coordinate system ：序列的第一个碱基设为数字1，如SAM,VCF,GFF,wiggle格式

• 0-based coordinate system ：序列的第一个碱基设为数字0，如BAM, BCFv2, BED, PSL格式

• 虽然bam是0-based，但使用samtools view来看bam，实际上是看到的sam文件，因此它是1-based（来自：https://www.biostars.org/p/282999/）。因此看bam结果时坐标是：3000563，而bed中坐标是：3000562

2 比对完成后，可以进行一些基础的统计

• 比对率（The fraction of aligned reads）：有多少原始的fq中的reads能够比对到基因组上，这个值应该尽可能高，一般大于70%。如果比对率小于50%，就可能存在文库污染或者其他问题

• 唯一比对reads的起始位点（The fraction of unique read start positions）：可以作为衡量文库复杂度的指标。在input样本中，reads应该“铺满”整个基因组，而ChIP样本应该在某些区域富集，导致某些区域存在很多相同的比对起始位点，但再多也只能计算一次。因此相比于input，ChIP样本中唯一比对reads的起始位点所占比例应该更低；另外相比于组蛋白标记，转录因子相关的比例也会更低，因为转录因子特异性结合区域更集中。一般ChIP-seq样本能得到80%以上的唯一比对起始位点（而RNA-seq只有15%左右）

2.1 统计比对的fq中有多少reads

# 我们这里使用的是clean fq
# 因为之前进行过滤的时候指定了：--fastqc（过滤完直接再次进行fastqc），所以会对每个样本单独创建一个目录，然后里面再分别存放过滤后的fq以及质控结果。因此看到下面的代码中路径会有点奇怪
gunzip -c $CLEAN_DIR/${name}.fastq.gz/${name}_trimmed.fq.gz | awk '(NR%4==2)' | wc -l

2.2 统计比对上了多少reads

bamToBed -i ${ALIGN_DIR}/${name}.sorted.bam | wc -l

2.3 统计唯一比对reads的起始位点

这里因为存在比对到基因组正负链的问题，所以需要根据比对结果对正链、负链分别计算

例如：

bamToBed -i H3K27AC_CHIP_TKO_1.sorted.bam | head -3
# 最后一列就是标注正负链
chr1    3000562 3000612 H3K27AC_CHIP_TKO_1.26986706 42  +
chr1    3000638 3000688 H3K27AC_CHIP_TKO_1.40152327 30  -
chr1    3000974 3001024 H3K27AC_CHIP_TKO_1.46142248 37  -

然后使用awk来进行统计，这里需要一些awk数组的知识

推荐一本书：Sed-awk-101Hacks，看Chapter 12. Awk Associative Arrays 就能了解awk中数组的知识；还有：A muggle's guide to AWK arrays 【https://www.datafix.com.au/BASHing/2019-06-07.html 】

这里简单介绍一下

• 数组存放键值对，格式是：arrayname[string]=value

• arrayname是数组的名称

• string是数组中的索引index（或者理解成其他编程语言中的键key也可以）

• 每个index都是唯一的，awk中的index不允许有重复，因此即使有好几个同样的index，存到awk中也只有一个【所以我们这里就是利用这个特性来对起始位点去重复】

• value就是对应string的值了

需要注意的是，其中的index很随意，可以是数值可以是字符串，数值也可以无序排列，因为即使是数值，awk也会转换成字符串，例如下面两个其实是一样的；

item[101]="HD Camcorder"
item["101"]="HD Camcorder"

 

awk中变量需要提前定义，而数组不需要。随想随用，当然更不需要提前定义数组有多大

了解完上面的基础以后，下面来看代码

解释：管道之前的命令很简单，就是把bam变成bed（上面有演示bed的格式）；然后如果第6列是+，就用第二列start position当起始位点；如果第6列是-，就用第三列end position当起始位点（正链的起点就是bed文件中的起点；负链的起点就是bed文件中的终点）

把染色体:位置 信息存为position，然后把它放在一个数组total中，保存起来，相同的染色体:位置 信息会只记录一次，也就是唯一的起始位点信息；但如何统计有多少个呢？因为前面提到数组记录的都是字符串（即使存入的是数值，也会被转为字符串），因此这里用了一个巧妙的办法，我现在有了total这个数组，我不管其中的内容，只要有一个index，就记为1，然后把total数组中的所有1加起来，就反映了有多少的index，也即是有多少的唯一起始位点

bamToBed -i ${ALIGN_DIR}/${name}.sorted.bam | awk '{if($6=="+"){position=$1":"$2}else if($6=="-"){position=$1":"$3};total[position]=1}END{print length(total)}'

2.4 最后上一个单行命令--多方位统计

bamToBed -i ${ALIGN_DIR}/${name}.sorted.bam | awk -v OFS="\t" -v sample=$name -v clean=$(gunzip -c $CLEAN_DIR/${name}.fastq.gz/${name}_trimmed.fq.gz | awk -v OFS="\t" '(NR%4==2)' | wc -l) 'BEGIN{max=0}{map++;if($6=="+"){position=$1":"$2}else if($6=="-"){position=$1":"$3};count[position]++; if(count[position]>max){max=count[position];maxPos=position}} END{totalPos=length(count); print sample,clean,map,map*100/clean,totalPos,totalPos*100/map, maxPos,count[maxPos],count[maxPos]*100/map}' >${ALIGN_DIR}/align_stat.txt

看似很复杂，它实际上统计了：样本名称（sample）、比对时clean fq中reads数（clean）、比对上的reads数（map）、比对率（map*100/clean）、唯一比对的起点总数（totalPos）、唯一比对的起点比例（totalPos*100/map）、哪个位点出现重复比对次数最多（maxPos）、该位点重复了多少次（count[maxPos]）、该位点重复出现的比例（count[maxPos]*100/map）

看结果，input样本唯一的起始位点占比最高（97%以上），其次是H3K27AC的样本（92%以上），最后是NRF1转录因子（80%多），和之前介绍唯一比对reads的起始位点是一致的