在WGBS測序中，我們選用bismark（http://www.bioinformatics.bbsrc.ac.uk/projects/bismark/）軟件對fq文件進(jìn)行mapping，該軟件基于bowtie或者bowtie2，將BS-seq reads C→T G→A分別轉(zhuǎn)化。再分別mapping到BS轉(zhuǎn)化過的基因組。得到的四個(gè)alignment結(jié)果來判斷最合適的unique alignment。同時(shí)軟件還可以統(tǒng)計(jì)出甲基化的類型如CpG、CHG或者CHH等。

1、軟件安裝

下載地址(github):https://github.com/FelixKrueger/Bismark.git

關(guān)聯(lián)軟件：samtools，bowtie/bowtie2

安裝方法：

git clonehttps://github.com/FelixKrueger/Bismark.git

tar-zxvf bismark_v0.15.0.tar.gz

2、軟件使用方法

bismark軟件分析BSSeq數(shù)據(jù)主要分為三個(gè)步驟：構(gòu)建基因組并創(chuàng)建bowtie2索引，4次DNAmapping，統(tǒng)計(jì)bam文件中的信息

※構(gòu)建基因組創(chuàng)建索引

選擇bismark軟件中的bismark_genome_preparation工具，需要給定bowtie2的路徑以及參考基因的路徑（包含fa和fai文件），操作代碼如下：

bismark_genome_preparation--path_to_bowtie2/usr/local/bowtie2/--verbose/data/genomes/homo_sapiens/GRCh37/

※mapping

選擇bismark工具進(jìn)行mapping，需要給出基因組路徑（第一步中--verbose路徑），用法如下

bismark[options]{-1-2|}

例如：bismark--bowtie2--path_to_bowtie/home/novelbio/software/bowtie2/../GRCH37/-1 filtered.1.fq.gz-2 filtered.2.fq.gz-o result/

雙端數(shù)據(jù)需要輸入-1與-2，單端數(shù)據(jù)直接輸入即可

※統(tǒng)計(jì)bam文件信息

選擇bismark_methylation_extractor工具進(jìn)行統(tǒng)計(jì)，用法如下：

用法：bismark_methylation_extractor[options]

測試使用代碼./bismark_methylation_extractor SRR534203_filtered.fq.gz_bismark_bt2.bam-s--gzip--bedGraph--genome_folder../ath_tair10/

其中輸入文件為第二部生成的bam文件，-s代表單端bam文件，-p代表雙端bam文件--gzip代表對結(jié)果文件進(jìn)行壓縮--bedGragh代表生成帶有甲基化率的bed文件

3、結(jié)果展示

※創(chuàng)建參考基因組

bismark將基因組的fa文件轉(zhuǎn)化為兩份，并分別使用bowtie2構(gòu)建索引

※Mapping結(jié)果

Mapping的結(jié)果中提供了reads的Mapping率，uniqueMapping情況，以及不同種類的甲基化程度

bam文件記錄展示：SRR534203.2_SN608_VA028:5:1101:24.50:89.20_length=50 16 chr3 6025118 42 49M*0 0

CTCACATCAATAAAATCTAATTCAATCCTCACCTCATCTTCAAAATAAA

FGIIIHDEJHDJIIGGIGIHHCHHGCJHFGCIHFJIHHHHHHFFDDD=1

NM:i:8 MD:Z:9G1G1G0G3G0G23G0G4

XM:Z:.........x.h.hh...xh.......................hh....XR:Z:CT XG:Z:GA

在每條reads記錄中提供了該位點(diǎn)的甲基化情況，在XM:Z:記錄中，"."代表不是甲基化位點(diǎn)，"z/Z"代表CpG位點(diǎn)，其中z代表未發(fā)生甲基化位點(diǎn)，Z代表發(fā)生甲基化的位點(diǎn)，

“x/X"代表CHG位點(diǎn)，"h/H"代表CHH位點(diǎn)，“u/U"代表CN或CHN位點(diǎn)

在生成的mapping Report結(jié)果匯總提到的發(fā)生CpG甲基化的位點(diǎn)個(gè)數(shù)其實(shí)就是全部reads中出現(xiàn)"Z"的數(shù)量總和，其他種類甲基化的算法也是一樣，甲基化率則根據(jù)（發(fā)生甲基化數(shù)量/（發(fā)生+未發(fā)生））計(jì)算

※統(tǒng)計(jì)結(jié)果

結(jié)果統(tǒng)計(jì)信息截圖

其中.bismark.cov.gz文件記錄了每個(gè)甲基化位點(diǎn)的覆蓋度，包含發(fā)生甲基化的reads數(shù)，未發(fā)生甲基化的reads數(shù)以及甲基化頻率，截圖如上右

bedGraph.gz文件以bed文件的格式記錄了甲基化位點(diǎn)的甲基化頻率

CpG_report.txt.gz文件記錄了位點(diǎn)，覆蓋度以及附近的位點(diǎn)信息

CpG_OT/OB文件記錄了每一條reads的CpG甲基化情況，OT代表original top strand，OB代表original bottom strand，文件截圖如下：