使用Bismark处理RRBS数据

Bismark是非常流行的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是建库的时候添加的，需要去除。如果是网上下载的数据，这个T可能已经被去除了，我们可以使用Fastqc来检查。

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

Drawing

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

Drawing

RRBS的Read 2测序结果的前三个碱基为TCA，在这个数据中第一个T已经去掉了，呈现出CAA的特征，还需要去除三个碱基。当我们需要批量处理很多RRBS的时候，也可以编写程序，自动预测需要的Trimming的碱基数据。下面以Read 2为例，简要示例。

In [12]:
! python /sibcb2/bioinformatics/code/Run_count_ACGTN.py -h

usage: Run_count_ACGTN.py [-h] -I FASTQGZ -O OUTFILE [-E TOPREADS]
                          [-L NUMBASE]

optional arguments:
  -h, --help            show this help message and exit
  -I FASTQGZ, --fastqGZ FASTQGZ
                        Input Fastq file.
  -O OUTFILE, --outFile OUTFILE
                        Output file.
  -E TOPREADS, --topReads TOPREADS
                        Number of reads used.
  -L NUMBASE, --NumBase NUMBASE
                        Number of bases used.

计算Read 2前20个碱基的分布：

python /sibcb2/bioinformatics/code/Run_count_ACGTN.py \
    -I dnaMethylation/Fastq/SRX1869508_2.fastq.gz \
    -O dnaMethylation/ACGTN/SRX1869508_2.txt

结果如下：

In [14]:
! cat dnaMethylation/ACGTN/SRX1869508_2.txt

A	C	G	T	N
12168	971415	4896	11522	0
985084	6679	4530	3708	0
979526	7615	4610	8250	0
608353	287701	2687	101260	0
384088	336536	16642	262469	266
454065	252679	66992	226265	0
408912	372513	15006	203570	0
491834	264010	48071	195720	366
509770	271352	24041	194670	168
423796	263545	37075	275129	456
482526	295701	23435	198282	57
514312	279645	30466	175427	151
428419	358126	21848	191479	129
404622	248514	36424	310441	0
405629	332348	17445	243951	628
478512	322871	19981	178278	359
369914	333226	19065	276799	997
415960	326218	62915	194684	224
371414	409879	18869	199188	651
368298	358408	18750	254463	82

可以清楚看到前三个碱基的特征为CAA，也可以用程序自动判断：

In [17]:
! Rscript /sibcb2/bioinformatics/code/Run_RRBS_find_trimming.R --countFile dnaMethylation/ACGTN/SRX1869508_2.txt --Read 2 --outFile /dev/stdout

SRX1869508_2	2	2

目前该程序只能判断Read 2的TCA的去除策略。在输出的第三列，程序自动判断需要Trimming的碱基数为2.

流程

我们编写了一个WDL的流程，运行这个样本的配置文件如下：

In [1]:
! cat dnaMethylation/json/SRX1869508.json

{
	"bismarkBS.file_basename":"SRX1869508",
	"bismarkBS.masterFolder":"/sibcb2/bioinformatics/JupyterLab/dnaMethylation/Bismark",
	"bismarkBS.Fastq_R1":"/sibcb2/bioinformatics/JupyterLab/dnaMethylation/Fastq/SRX1869508_1.fastq.gz",
	"bismarkBS.Fastq_R2":"/sibcb2/bioinformatics/JupyterLab/dnaMethylation/Fastq/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/SRX1869508.json

生成的报告文件：Bismark report of SRX1869508。生成的结果文件目录如下：

In [3]:
! tree /sibcb2/bioinformatics/JupyterLab/dnaMethylation/Bismark/SRX1869508

/sibcb2/bioinformatics/JupyterLab/dnaMethylation/Bismark/SRX1869508
|-- Call
|   |-- CHG_context_SRX1869508_nSorted.txt.gz
|   |-- CHH_context_SRX1869508_nSorted.txt.gz
|   |-- CpG_context_SRX1869508_nSorted.txt.gz
|   |-- SRX1869508_nSorted.bedGraph.gz
|   `-- SRX1869508_nSorted.bismark.cov.gz
|-- SRX1869508.bam
|-- SRX1869508.bam.bai
|-- SRX1869508_report.html
`-- report
    |-- SRX1869508_1.fastq.gz_trimming_report.txt
    |-- SRX1869508_2.fastq.gz_trimming_report.txt
    |-- SRX1869508_M-bias.txt
    |-- SRX1869508_PE_report.txt
    `-- SRX1869508_splitting_report.txt

2 directories, 13 files