R-loop数据分析之R-ChIP(数据预处理)

文件重命名

我们需要对下载的SRRXXXXX文件进行重命名，毕竟有意义的命名才能方便后续展示。那么，应该如何做呢？

首先，你需要将GSE97072页面的中Samples这部分的内容复制到一个文本文件中（我将其命名为sample_name.txt)，分为两列，第一列是GSM编号，第二列是样本的命名。

sample name

注：这里面有一个希腊字符在不同系统表示有所不同，所以我在复制之后手动修改。

此外，你还需要将里面的(-)和(+)进行替换，因为括号在shell里有特殊含义， 为了保证命名的连贯性，我将(-)替换成-neg，将(+)替换成-pos.

 sed -i 's/(-)/-neg/; s/(+)/-pos/' sample_name.txt

随后，我们需要从download_table.txt中提取出SRR编号和GSM编号的对应的关系，这个需要用到Linux的文本处理命令

paste <(egrep -o 'GSM[0-9]{6,9}' download_table.txt ) <(egrep -o 'SRR[0-9]{6,9}' download_table.txt) >  gsm_srr.txt
join <(sort gsm_srr.txt ) <(head -n 24 sample_name.txt | sort) > gsm_srr_sample_name.txt

注: 样本一共有25行，而我们只下载了24个数据，所以删除了最后一行的"K562-WKKD-V5chip"

最后，就是根据生成的文件对样本进行重命名了。

awk '{print "rename " $2 " " $3   " analysis/0-raw-data/" $2 "*.gz"}' gsm_srr_sample_name.txt |bash

这行代码看起来有点复杂，但是其实做的事情就是构建了一系列rename的命令行，然后在bash下执行。

R-ChIP分析

数据预处理

这一部分主要是将序列比对后原来的参考基因组上，标记重复，并且去掉不符合要求的比对。

让我们先写一个比对脚本将序列比对到参考基因组上，脚本命名为03.r_chip_align.sh，存放在scripts下

#!/bin/bash

set -e
set -u
set -o pipefail

# configuration
threads=8
index=index/hg19
FQ_DIR="analysis/0-raw-data"
ALIGN_DIR="analysis/2-read-align"
LOG_DIR="analysis/log"
TMP_DIR="analysis/tmp"

mkdir -p ${ALIGN_DIR}
mkdir -p ${LOG_DIR}
mkdir -p ${TMP_DIR}

samples=${1?missing sample file}

exec 0< $samples
# alignment
while read id;
do
    if [ ! -f ${FQ_DIR}/$id.bt2.done ]
    then
    echo "bowtie2 --very-sensitive-local --mm -p $threads -x $index -U ${FQ_DIR}/$id.fastq.gz 2> ${LOG_DIR}/$id.bt2.log | \
    samtools sort -@ 2 -m 1G -T ${TMP_DIR}/${id} -o ${ALIGN_DIR}/${id}.sort.bam \
    && touch ${ALIGN_DIR}/$id.bt2.done" | bash
    fi
done

这个脚本的前半部分都在定义各种变量，而#alignment标记的后半部分则是实际运行的比对命令

你会发现我的实际运行部分脚本有点奇怪，我没有直接运行比对，而是用echo通过管道传递给了bash执行。

这样做的原因是, 当我用sed -i 's/| bash/#| bash/ 03.r_chip_align.sh'将| bash替换成#| bash后，运行这个脚本就可以检查代码是否正确，然后再用sed-i s/#| bash/| bash/03.r_chip_align.sh修改回来 。 后面的脚本同理。

接着再写一个脚本用于标记重复04.markdup.sh，也是放在scripts下

#!/bin/bash

set -e
set -u
set -o pipefail

# configuration
threads=8
index=index/hg19
FQ_DIR="analysis/0-raw-data"
ALIGN_DIR="analysis/2-read-align"
LOG_DIR="analysis/log"
TMP_DIR="analysis/tmp"

mkdir -p ${ALIGN_DIR}
mkdir -p ${LOG_DIR}
mkdir -p ${TMP_DIR}

samples=${1?missing sample file}

exec 0< $samples
# mark duplication
while read id;
do
    if [ ! -f ${ALIGN_DIR}/${id}.mkdup.done ]
    then
    echo "sambamba markdup -t $threads ${ALIGN_DIR}/${id}.sort.bam ${ALIGN_DIR}/${id}.mkdup.bam \
    && touch ${ALIGN_DIR}/${id}.mkdup.done" | bash
    fi
done

再准备一个脚本用于去除不符合要求的比对，命名为05.bam_filter.sh

#!/bin/bash

set -e
set -u
set -o pipefail

# configuration
threads=8
index=index/hg19
FQ_DIR="analysis/0-raw-data"
ALIGN_DIR="analysis/2-read-align"
LOG_DIR="analysis/log"
TMP_DIR="analysis/tmp"

mkdir -p ${ALIGN_DIR}
mkdir -p ${LOG_DIR}
mkdir -p ${TMP_DIR}

samples=${1?missing sample file}

exec 0< $samples
# filter
while read id;
do
    if [ ! -f ${ALIGN_DIR}/${id}.flt.done ]
    then
    echo "
    samtools view -@ threads -bF 1804 -q 30 ${ALIGN_DIR}/${id}.mkdup.bam -o ${ALIGN_DIR}/${id}.flt.bam\
    && samtools index ${ALIGN_DIR}/${id}.flt.bam \
    && touch ${ALIGN_DIR}/${id}.flt.done "| bash
    fi
done

最后新建一个samples01.txt用于存放将要处理的样本名

HKE293-D210N-V5ChIP-Rep1
HKE293-D210N-Input-Rep1
HKE293-D210N-V5ChIP-Rep2
HKE293-D210N-Input-Rep2
HKE293-D210N-V5ChIP-Rep3
HKE293-D210N-Input-Rep3
HKE293-WKKD-V5ChIP
HKE293-WKKD-Input
HKE293-delta-HC-V5ChIP

这样子就可以依次运行脚本了

• 比对: bash scripts/03.r_chip_align.sh samples01.txt
• 标记重复：bash scripts/04.bam_markdup.sh samples01.txt
• 去除不符合要求的比对: bash scripts/05.bam_filter.sh samples01.txt

处理完之后可以对每个样本都进行一次统计，包括如下信息：

• 处理前的原始reads数
• 处理后对唯一比对reads数
• 唯一比对reads数所占原始reads数的比例

这个工作同样可以用shell脚本完成, 脚本为06.sample_align_stat.sh

#!/bin/bash

set -e
set -u
set -o pipefail

samples=${1?missing sample file}
threads=8
ALIGN_DIR="analysis/2-read-align"

echo -e "Experiment \t Raw Reads \t Uniquely mapped Reads \t ratio"
exec 0< $samples

while read sample;
do
     total=$( samtools view -@ ${threads} -c ${ALIGN_DIR}/${sample}.sort.bam )
     unique=$( samtools view -@ ${threads} -c ${ALIGN_DIR}/${sample}.flt.bam )
     ratio=$( echo "scale=2; 100 * $unique / $total " | bc )
     echo -e "$sample \t $total \t $unique \t $ratio %"
done

运行方法是bash scripts/06.sample_align_stat.sh samples01.txt > results/library_stat.txt，运行结果如下，和原本的Table S2对比，你会发现结果基本一致，有出入的地方我推测是标记重复这一步所用软件不同。