使用Bsmap处理RRBS数据

BSMAP是非常流行的DNA甲基化分析软件。本流程将介绍使用该流程的基本过程。使用的数据来为SRX1869508， 自NCBI GEO。

确定Trimming的参数

如下图所示，RRBS在Read 1的5-prime需要去掉T，在Read 2的5-prime要去掉TCA。

Drawing

Read 1的5-prime前几个碱基是T[CT]GG, 第一个T是建库的时候添加的，需要去除。

使用FastQC

我们可以使用Fastqc来检查每个测序位置的碱基分布，以此来决定需要Trimming的碱基个数。

In [4]:
# fastqc /sibcb2/bioinformatics/JupyterLab/dnaMethylation/Fastq/SRX1869508_1.fastq.gz

Drawing

可有看出，第一个碱基T已经被去除了，前三个碱基的模式为[TC]GG，不要去除更多碱基，我们可以是用Fastqc来检查Read 2.

Drawing

RRBS的Read 2测序结果的前三个碱基为TCA，在这个数据中第一个T已经去掉了，呈现出CAA的特征，还需要去除三个碱基。

使用Bismark进行预比对

Bismark是非常流行的DNA甲基化分析软件，能够产生丰富的质控信息，详细见Bismark示例。我们可以使用Bismark来分析这个数据中的一小部分序列（一百万条）进行质控。

In [6]:
# zcat dnaMethylation/Fastq/SRX1869508_1.fastq.gz |head -4000000 | gzip -c > dnaMethylation/FastqSubset/SRX1869508_1.fastq.gz
# zcat dnaMethylation/Fastq/SRX1869508_2.fastq.gz |head -4000000 | gzip -c > dnaMethylation/FastqSubset/SRX1869508_2.fastq.gz

根据FastQC的结果Read 1的5-prime去除0个碱基，Read 2的5-prime去除3个碱基；在3-prime我们默认统一去掉3个碱基，配置文件如下：

In [4]:
! cat dnaMethylation/json/SRX1869508Subset.json

{
	"bismarkBS.file_basename":"SRX1869508",
	"bismarkBS.masterFolder":"/sibcb2/bioinformatics/TBC/dnaMethylation/BismarkQC",
	"bismarkBS.Fastq_R1":"/sibcb2/bioinformatics/TBC/dnaMethylation/FastqSubset/SRX1869508_1.fastq.gz",
	"bismarkBS.Fastq_R2":"/sibcb2/bioinformatics/TBC/dnaMethylation/FastqSubset/SRX1869508_2.fastq.gz",
	"bismarkBS.paired":"true",
	"bismarkBS.RRBS":"true",
	"bismarkBS.clip_R1":"0",
	"bismarkBS.three_prime_clip_R1":"3",
	"bismarkBS.clip_R2":"3",
	"bismarkBS.three_prime_clip_R2":"3",	
	"bismarkBS.genome_folder":"/sibcb2/bioinformatics/iGenome/Bismark/hg19",
	"bismarkBS.threadN":"12"
}

将该任务提交到服务器：

java -Xmx1g -jar /sibcb2/bioinformatics/software/cromwell/cromwell-50.jar \
    run /sibcb2/bioinformatics/code/WDL/bismarkBS.wdl \
    --inputs /sibcb2/bioinformatics/JupyterLab/dnaMethylation/json/SRX1869508Subset.json

一般只需要几分钟的时间，Bismark产生的报告文件如下：1M reads的Bismark报告。可以看出，这个结果和所有reads的Bismark报告在很大程度上吻合，也说明了采用1M的reads进行质量控制的可行性。

Bsmap比对

真正的数据分析步骤利用BSMAP来完成，我们使用python编写了一个方便的流程。

In [9]:
! python /sibcb2/bioinformatics/code/Run_BSMAP.py -h

usage: Run_BSMAP.py [-h] -R1 FASTQ_R1 [-R2 FASTQ_R2] [-C1 CLIP_R1]
                    [-C2 CLIP_R2] [-T1 THREE_PRIME_CLIP_R1]
                    [-T2 THREE_PRIME_CLIP_R2] -G GENOMEFILE [-N THREADN]
                    [-n N] -O OUT -tag TAG -tmp TMP [-cgGR CGGR] [-CGI CGI]
                    [-rrbs] [-amplicon] [-runQC]

optional arguments:
  -h, --help            show this help message and exit
  -R1 FASTQ_R1, --Fastq_R1 FASTQ_R1
                        Read 1 fastq file.
  -R2 FASTQ_R2, --Fastq_R2 FASTQ_R2
                        Read 2 fastq file.
  -C1 CLIP_R1, --clip_R1 CLIP_R1
                        Number of bases to be trimmed on 5' of Read 1.
  -C2 CLIP_R2, --clip_R2 CLIP_R2
                        Number of bases to be trimmed on 5' of Read 2.
  -T1 THREE_PRIME_CLIP_R1, --three_prime_clip_R1 THREE_PRIME_CLIP_R1
                        Number of bases to be trimmed on 3' of Read 1.
  -T2 THREE_PRIME_CLIP_R2, --three_prime_clip_R2 THREE_PRIME_CLIP_R2
                        Number of bases to be trimmed on 3' of Read 2.
  -G GENOMEFILE, --genomeFile GENOMEFILE
                        Genome fasta file for alignment.
  -N THREADN, --threadN THREADN
                        Number of threads.
  -n N, --n N           Set mapping strand information.
  -O OUT, --out OUT     Output folder.
  -tag TAG, --tag TAG   Sample ID.
  -tmp TMP, --tmp TMP   tmpFolder.
  -cgGR CGGR, --cgGR CGGR
                        CpG site R data.
  -CGI CGI, --CGI CGI   CpG Island BED file.
  -rrbs, --rrbs         RRBS data
  -amplicon, --amplicon
                        Amplicon data
  -runQC, --runQC       whether run QC script

就这个具体的数据而言，具体参数如下：

python /sibcb2/bioinformatics/code/Run_BSMAP.py \
    -R1 dnaMethylation/Fastq/SRX1869508_1.fastq.gz \
    -R2 dnaMethylation/Fastq/SRX1869508_2.fastq.gz \
    -C1 0 \
    -T1 3 \
    -C2 3 \
    -T2 3 \
    -G /sibcb2/bioinformatics/iGenome/Bismark/hg19/hg19.fa \
    -N 16 \
    -O dnaMethylation/BSMAP \
    -tag SRX1869508 \
    -tmp tmp \
    --cgGR /sibcb2/bioinformatics/iGenome/CpG/hg19_CpG_sites.RData \
    --CGI /sibcb2/bioinformatics/iGenome/cpgIsland/out/hg19_cpgIsland.bed \
    --rrbs \
    --runQC
In [3]:
! tree dnaMethylation/BSMAP/SRX1869508

dnaMethylation/BSMAP/SRX1869508
|-- RData
|   |-- SRX1869508_CHG.RData
|   |-- SRX1869508_CHH.RData
|   `-- SRX1869508_CpG.RData
|-- SRX1869508.bam
|-- SRX1869508.bam.bai
|-- mCall
|   |-- SRX1869508_CHG.bedGraph
|   |-- SRX1869508_CHH.bedGraph
|   `-- SRX1869508_CpG.bedGraph
`-- report
    |-- SRX1869508_1.fastq.gz_trimming_report.txt
    |-- SRX1869508_2.fastq.gz_trimming_report.txt
    |-- SRX1869508_BSMAP_report.txt
    `-- SRX1869508_QC_summary.txt

3 directories, 12 files

其中比对的QC信息：

In [14]:
! cat dnaMethylation/BSMAP/SRX1869508/report/SRX1869508_BSMAP_report.txt

[bsmap] @Thu Sep 10 21:14:06 2020 	loading reference file: /sibcb2/bioinformatics/iGenome/Bismark/hg19/hg19.fa 	(format: FASTA)
[bsmap] @Thu Sep 10 21:14:33 2020 	93 reference seqs loaded, total size 3137161264 bp. 27 secs passed
[bsmap] @Thu Sep 10 21:15:22 2020 	create seed table. 76 secs passed
[bsmap] @Thu Sep 10 21:15:22 2020 	Pair-end alignment(16 threads),
	Input read file #1: tmp/21e61e22-f38c-11ea-b14b-6c92bfc12f02/trimmed/SRX1869508_R1_val_1.fq.gz 	(format: gzipped FASTQ)
	Input read file #2: tmp/21e61e22-f38c-11ea-b14b-6c92bfc12f02/trimmed/SRX1869508_R2_val_2.fq.gz 	(format: gzipped FASTQ)
	Output: STDOUT	 (format: SAM)
[bsmap] @Thu Sep 10 23:41:21 2020 	total read pairs: 167236501 	total time consumed:  8835 secs
	aligned pairs: 114133411 (68.2%), unique pairs: 114133411 (68.2%), suppressed non-unique pairs: 25997063 (15.5%)
	unpaired read #1: 4674665 (2.8%), unique reads: 4674665 (2.8%), suppressed non-unique reads: 29585193 (17.7%)
	unpaired read #2: 7169678 (4.3%), unique reads: 7169678 (4.3%), suppressed non-unique reads: 30185514 (18.0%)

综合QC信息：

In [2]:
! cat dnaMethylation/BSMAP/SRX1869508/report/SRX1869508_QC_summary.txt

Feature	Value
COV_CGI_N_5x	1396717
COV_CGI_N_10x	1347682
COV_CGI_median	150
COV_CGI_mean	173.8
COV_nonCGI_N_5x	3037958
COV_nonCGI_N_10x	2531038
COV_nonCGI_median	3
COV_nonCGI_mean	45.1
CpG_CGI_methylation	0.004518
CpG_nonCGI_methylatoin	0.8506
CHG_CGI_methylation	0.0006193
CHG_nonCGI_methylation	0.001729
CHH_CGI_methylation	0.0005513
CHH_nonCGI_methylation	0.0009982