官方手冊（2020 Sep更新）
https://github.com/FelixKrueger/Bismark/tree/master/Docs

（一）在開始前你需要了解的：

Bismark是用Perl編寫的，并在命令行里運行。另外，你需要安裝Bowtie 2 或者HISAT2 递雀。 對于0.14或者更高的版本，Bismark可以使用并行運行弊知，同時進(jìn)行比對和甲基化提取步驟握牧。你可以搜索--parallel/--multicore來看更多的細(xì)節(jié)容诬。

首先你需要下載參考基因組，并把它放在一個文件夾里沿腰。參考基因組可以在 Ensembl或者NCBI 網(wǎng)站下載览徒。對于下面講到的例子，你需要下載Homo sapiens基因組颂龙。Bismark支持FastA格式的基因組序列习蓬，文件后綴可以是.fa或.fasta〈肭叮可以輸入單個fasta或者多個fasta文件躲叼。

教程的示例數(shù)據(jù)是test_dataset.fastq，下載地址here(it contains 10,000 reads in FastQ format, Phred33 qualities, 50 bp long reads, from a human directional BS-Seq library)企巢。

（二）Bismark-主要信息

（1）什么是Bismark?

Bismark是一套工具枫慷，可高效率的分析亞硫酸氫鹽序列(BS-Seq)數(shù)據(jù)。Bismark用參考基因組進(jìn)行比對，同時進(jìn)行胞嘧啶甲基化提取流礁。Bismark是用Perl編寫的涕俗，在命令行運行。用亞硫酸氫鹽處理的reads用短Reads對準(zhǔn)軟件Bowtie 2或HISAT2進(jìn)行比對神帅。因此再姑，需要在你的電腦上安裝Bowtie 2(或HISAT2)。

（2）Bismark的安裝

由于我使用的是服務(wù)器里的軟件找御，所以直接調(diào)用就行了元镀。如果想在自己的電腦上安裝Bismark，可以在這個網(wǎng)站下載：here霎桅。
然后解壓：

tar xzf bismark_v0.X.Y.tar.gz

（3）硬件要求

Bismark將參考基因組保存在內(nèi)存中栖疑，除此之外，還運行了多達(dá)四個Bowtie 2的并行實例滔驶。內(nèi)存的使用取決于參考基因組的大小遇革。對于一個大的真核生物基因組(人類或小鼠)，我們的經(jīng)驗大約要用12GB的內(nèi)存揭糕。因此萝快，我們建議在具有5個CPU核和至少12 GB RAM的機(jī)器上運行Bismark。Bowtie 2的內(nèi)存需求有點大(可能是為了允許gapped比對)著角。在使用Bowtie 2運行Bismark時揪漩，我們推薦至少有5個核和> 16GB RAM的系統(tǒng)。
比對速度很大程度上取決于所使用的read長度和比對參數(shù)吏口。允許許多mismatches和使用較短的seed長度往往是相當(dāng)慢的奄容。

（4）這個軟件支持什么樣的BS-Seq文件？

Bismark支持在以下條件下亞硫酸氫鹽處理的reads(全基因組shotgun BS-Seq (WGSBS)产徊、reducated -representation BS-Seq (RRBS)或PBAT-Seq(后亞硫酸氫鹽接頭標(biāo)簽)進(jìn)行比對:

1.序列格式可以是FastQ或FastA
2.單端或雙端reads
3.輸入文件可以被解壓的或被gzip壓縮的(以.gz結(jié)尾)
4.可變reads長度
5.定向或非定向的BS-Seq文庫

下面是一個比對模式的完整列表：

注意：Bismark目前只支持Illumina平臺的數(shù)據(jù)昂勒，尚不支持SOLiD平臺的reads。

（5）Bismark是如何工作的?

序列reads首先由bisulfite進(jìn)行轉(zhuǎn)換(C->T)和反向(G->A轉(zhuǎn)換正向鏈)版本囚痴，然后再與同樣轉(zhuǎn)換過的基因組進(jìn)行比對(也包括C->T和G->A轉(zhuǎn)換)叁怪。這樣比對后就有4種結(jié)果，選取最佳比對的那一種后深滚，然后將其與正侈忍罚基因組序列進(jìn)行比較，從而推斷reads中所有胞嘧啶位置的甲基化狀態(tài)痴荐。如果一個比對具有唯一的最佳比對得分(如as:i字段報告的那樣)血柳，則認(rèn)為該read是惟一比對。如果一個read的比對產(chǎn)生了幾個具有相同數(shù)量的 mismatches或相同比對分?jǐn)?shù)(如AS:i字段)生兆，那么這個read(或一個read-pair)將被完全丟棄难捌。

光看文字可能比較難理解為什么會有4種結(jié)果膝宁，下面的圖應(yīng)該能比較直觀的方便理解：

（6）Bismark比對和提取甲基化報告

完成后，Bismark生成一個包含以下信息的運行報告:

1.使用的比對參數(shù)總結(jié)
2.已分析的序列數(shù)
3.具有唯一最佳比對的序列數(shù)(mapping efficiency)

4.統(tǒng)計總結(jié)唯一最佳比對來自哪一條亞硫酸氫鹽鏈
5.分析的胞嘧啶數(shù)目
6.甲基化和未甲基化胞嘧啶的數(shù)目

在CpG根吁、CHG或CHH環(huán)境中胞嘧啶甲基化百分比(其中H可以是A员淫、T或C)。該百分比根據(jù)公式分別計算:

% methylation (context) = 100 * methylated Cs (context) / (methylated Cs (context) + unmethylated Cs (context)).
需要強(qiáng)調(diào)的是击敌，甲基化百分比值(背景)只是在映射步驟中直接執(zhí)行的一個非常粗略的計算介返。經(jīng)過后處理或過濾后的實際甲基化水平可能有所不同。

（三）運行Bismark

運行Bismark分為三個主要步驟:
首先沃斤，需要對感興趣的基因組進(jìn)行亞硫酸氫鹽轉(zhuǎn)換和并構(gòu)建索引圣蝎，以允許Bowtie進(jìn)行比對。對于每個基因組衡瓶，這一步只需要執(zhí)行一次徘公。注意，Bowtie 2和HISAT2需要不同的索引哮针，因為它們的索引不兼容关面。

Bismark read比對步驟。只需指定一個要分析的文件十厢，參考基因組和比對參數(shù)缭裆。Bismark將生成一個結(jié)合了比對/甲基化call的輸出(默認(rèn)是BAM格式)以及一個運行統(tǒng)計報告。

Bismark甲基化提取寿烟。這個步驟是可選的，它將從Bismark比對輸出中提取甲基化信息辛燥。運行這個額外的步驟允許將甲基化信息分解到不同的背景中筛武，獲得特異鏈的甲基化信息并提供一些過濾選項。你還可以選擇將甲基化信息排序到bedGraph/coverage文件中挎塌，或者進(jìn)一步處理它們到全基因組胞嘧啶甲基化報告中徘六。

下面幾節(jié)將更詳細(xì)地通過示例描述這些步驟。

(I) Bismark 基因組準(zhǔn)備

這個腳本只需要運行一次榴都，以準(zhǔn)備感興趣的亞硫酸氫鹽比對基因組待锈。你需要指定一個目錄，其中包含你想要比對reads的基因組(請注意嘴高，bismark_genome_prepare腳本期望該文件夾中包含F(xiàn)astA文件竿音，擴(kuò)展名為.fa或. FastA，每個文件有單個或多個序列條目)拴驮。Bismark會在這個目錄中創(chuàng)建兩個單獨的文件夾春瞬，一個用于C->T轉(zhuǎn)換基因組，另一個用于G->A轉(zhuǎn)換基因組套啤。在創(chuàng)建了C->T和G->A版本的基因組后宽气，它們將被使用bowtie2-build(或hisat2-build)并行構(gòu)建索引。一旦C->T和G->A基因組索引都被創(chuàng)建，你就不需要再次使用該腳本萄涯。這一步根據(jù)基因組的大小绪氛，有可能需要幾個小時來運行！

注：這一步我使用了服務(wù)器里的128G內(nèi)存涝影，運行大約2小時

請再次注意枣察，Bowtie 2和HISAT2的索引是不兼容的!
這里我們使用的是bowtie2。
先看一下基本格式：

bismark_genome_preparation [options] <path_to_genome_folder>

代碼：

$ bismark_genome_preparation --path_to_aligner /gpfs/share/apps/bowtie2/2.4.1/bin/ /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta

因為我是在服務(wù)器里運行的袄琳，所以我可以指定調(diào)用多少G的內(nèi)存以及多少個cores（使用自己電腦運行的同學(xué)請忽略下面代碼询件，直接跳到下面“(II) Bismark 比對步驟”，這里代碼僅作為我自己的筆記使用）：

$ srun --mem=128G --cpus-per-task=1 --time=8:00:00 --pty bismark_genome_preparation --path_to_aligner /gpfs/share/apps/bowtie2/2.4.1/bin/ /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta

或：

$ sbatch genome_preparation.sh

#!/bin/bash
#SBATCH --job-name=bismark_genome_preparation           # Job name
#SBATCH --mail-type=END,FAIL                            # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=yan.fang@***.org                    # Where to send mail
#SBATCH --ntasks=1                                      # Run on a single CPU
#SBATCH --mem=128gb                                     # Job memory request
#SBATCH --time=08:00:00                                 # Time limit hrs:min:sec
#SBATCH --output=/gpfs/home/fangy04/log_files/bismark_%j.log           # Standard output and error log
#SBATCH -p cpu_short

module load bismark/0.22.1
module load python/cpu/3.6.5
module load bowtie2/2.4.1
module load perl/5.28.0

bismark_genome_preparation --path_to_aligner /gpfs/share/apps/bowtie2/2.4.1/bin/ /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta

這一步運行完成后唆樊，看一下都有哪些輸出文件：

$ tree #這是我存放原始基因組的文件夾
.
`-- WholeGenomeFasta
    |-- Bisulfite_Genome #這個文件夾是運行后生成的宛琅，可以看到里面有兩個子文件夾
    |   |-- CT_conversion #CT轉(zhuǎn)換的基因組文件夾
    |   |   |-- BS_CT.1.bt2
    |   |   |-- BS_CT.2.bt2
    |   |   |-- BS_CT.3.bt2
    |   |   |-- BS_CT.4.bt2
    |   |   |-- BS_CT.rev.1.bt2
    |   |   |-- BS_CT.rev.2.bt2
    |   |   `-- genome_mfa.CT_conversion.fa
    |   `-- GA_conversion #GA轉(zhuǎn)換的基因組文件夾
    |       |-- BS_GA.1.bt2
    |       |-- BS_GA.2.bt2
    |       |-- BS_GA.3.bt2
    |       |-- BS_GA.4.bt2
    |       |-- BS_GA.rev.1.bt2
    |       |-- BS_GA.rev.2.bt2
    |       `-- genome_mfa.GA_conversion.fa
    |-- GenomeSize.xml
    |-- genome.dict
    |-- genome.fa
    `-- genome.fa.fai

4 directories, 18 files

(II) Bismark 比對步驟

這一步代表了亞硫酸氫鹽比對和call甲基化的實際部分。Bismark要求用戶只需指定兩件事:

1.包含感興趣的基因組的目錄逗旁。這個文件夾必須包含未修改的基因組(如.fa或.fasta文件)嘿辟，以及在Bismark基因組準(zhǔn)備步驟中生成的兩個亞硫酸氫鹽基因組子目錄。
2.要分析的序列文件(FastQ或FastA格式)片效。
所有其他參數(shù)都是可選的红伦。

對于每個序列文件或每組配對的paired-end測序文件，Bismark生成一個比對和甲基化call輸出文件淀衣，以及一個報告文件昙读，詳細(xì)說明比對和甲基化call統(tǒng)計信息，用于記錄保存膨桥。

在運行Bismark之前蛮浑，我們建議花些時間使用FastQC對原始序列文件進(jìn)行質(zhì)量控制。FastQC可能能夠發(fā)現(xiàn)與你的BS-Seq文件相關(guān)的異常只嚣，如高堿基calling錯誤率或污染序列沮稚，如PCR引物或Illumina接頭。許多誤差來源會對比對效率和/或比對位置產(chǎn)生不利影響册舞，因此也可能影響call甲基化和得出的結(jié)論蕴掏。有必要的話，需要進(jìn)行Trimming调鲸。但是這里由于是示例文件盛杰，這一步就不做了。在實際的數(shù)據(jù)處理中线得，這一步是要進(jìn)行的饶唤。

如果沒有指定其他選項，Bismark將使用一組默認(rèn)值贯钩，其中一些是:

（1）使用Bowtie 2

1.使用Bowtie 2是默認(rèn)模式
2.如果沒有指定到Bowtie 2的路徑募狂，則假定bowtie2在PATH中
3.標(biāo)準(zhǔn)的比對使用multi-seed長度為20bp且0錯配办素。可以使用選項-L和-N修改這些參數(shù)
4.標(biāo)準(zhǔn)比對使用默認(rèn)的最小比對評分函數(shù)L,0祸穷，-0.2性穿，即f(x) = 0 + -0.2 * x(其中x是讀取長度)。對于75bp的read雷滚，這意味著在比對失效之前需曾，read的最低比對分?jǐn)?shù)可以是-15。這大致等于read中出現(xiàn)2次錯配或1-2個indel(或兩者的組合)祈远。使用--score_min <func>設(shè)置嚴(yán)格度呆万。

即使用戶沒有被要求指定額外的比對選項，但通常建議這樣做(例如车份，當(dāng)默認(rèn)參數(shù)太嚴(yán)格時)谋减。要查看選項的完整列表，請在命令行中鍵入bismark --help扫沼。

（2）定向的BS-Seq文庫 (默認(rèn))

對于一個給定的locus出爹，亞硫酸氫鹽處理DNA和隨后的PCR擴(kuò)增可以產(chǎn)生四條(亞硫酸氫鹽轉(zhuǎn)換)鏈。根據(jù)使用的接頭缎除，BS-Seq庫可以用兩種不同的方式構(gòu)建:

1.如果庫是定向的严就，則只對原top鏈(OT)或原bottom鏈(OB)的(亞硫酸氫鹽轉(zhuǎn)換)版本進(jìn)行測序。即使與OT (CTOT)互補(bǔ)鏈和OB (CTOB)互補(bǔ)鏈在BS-PCR步驟中產(chǎn)生器罐，它們也不會被測序梢为，因為它們在5 '端攜帶了錯誤的接頭。默認(rèn)情況下轰坊，Bismark只對OT和OB鏈執(zhí)行2次read比對抖誉，因此忽略了來自互補(bǔ)鏈的比對，因為它們理論上不應(yīng)該出現(xiàn)在BS-Seq庫中衰倦。
2.或者，也可以構(gòu)建BS-Seq文庫旁理，使在BS-PCR中產(chǎn)生的所有4條不同的鏈樊零，以大致相同的可能性進(jìn)入測序文庫。在本例中孽文，所有四條鏈(OT驻襟、CTOT、OB芋哭、CTOB)都可以產(chǎn)生有效的比對沉衣，這個庫稱為非定向的。指定--non_directional將指示Bismark使用所有四種比對的輸出减牺。

再總結(jié)一下:與原始top鏈的比對或與原始top鏈的互補(bǔ)鏈(OT和CTOT)都將產(chǎn)生top鏈上胞嘧啶的甲基化信息豌习。與原始bottom鏈或與原始bottom鏈互補(bǔ)的鏈(OB和CTOB)比對，都會產(chǎn)生bottom鏈胞嘧啶的甲基化信息，也就是說禀苦，它們似乎會產(chǎn)生參考基因組top鏈G位置的甲基化信息贴唇。

有關(guān)如何提取四種不同比對鏈的甲基化信息的更多信息，請參閱下面關(guān)于Bismark methylation extractor的部分栋艳。

這一步的代碼基本格式：

bismark [options] --genome <genome_folder> {-1 <mates1> -2 <mates2> | <singles>}

在我們這個例子里恰聘，示例數(shù)據(jù)是單端測序，所以代碼是：

$ bismark --genome /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta sample.fastq.gz

（3）Bismark輸出的文件是什么樣的吸占？

從0.6.x版本以后晴叨，輸出的都是BAM/SAM格式。
運行后生成的文件：

（4）Bismark BAM/SAM輸出文件

默認(rèn)情況下矾屯，Bismark為所有比對模式生成SAM輸出兼蕊。請注意，報告的質(zhì)量值是用Sanger格式(Phred 33)編碼的问拘。
每一列如下：

QNAME：序列ID
FLAG:這個FLAG試圖考慮到亞硫酸氫鹽讀取的來源(這不同于普通的DNA比對FLAG!)
RNAME：染色體
POS：起始位置
MAPQ：為Bowtie 2和HISAT2計算
CIGAR
RNEXT
PNEXT
TLEN
SEQ
QUAL：Phred33
NM-tag
MD-tag：每個錯配的堿基(14)XM-tag（甲基化calling字符）
XR-tag：read比對轉(zhuǎn)換狀態(tài)(16)XG-tag(基因組比對轉(zhuǎn)換狀態(tài))

（5）數(shù)據(jù)可視化

要查看比對上的reads的位置遍略，可以使用染色體、開始和結(jié)束位置將Bismark輸出文件導(dǎo)入基因組查看器骤坐，比如SeqMonk(這對于識別基因組中顯示大量比對上的Reads的區(qū)域非常有用)绪杏。比對的輸出也可以用于后續(xù)步驟如去重(每個位置只允許1 條read，去除PCR artefacts)或過濾亞硫酸氫鹽轉(zhuǎn)換相關(guān)non-bisulfite錯配的數(shù)量纽绍。

這里我就先不練習(xí)如何使用SeqMonk可視化了蕾久，如果后續(xù)需要我再深入學(xué)習(xí)

與亞硫酸氫鹽轉(zhuǎn)換相關(guān)的非亞硫酸氫鹽錯配是在基因組中有T，而在BS-read中有C的錯配位置拌夏。這種錯配可能是由于亞硫酸氫鹽read比對的方式造成的僧著。為了不引入甲基化reads的偏倚，包含這種錯配的reads不會自動從比對輸出中刪除障簿。值得注意的是盹愚，這些位置即使沒有執(zhí)行甲基化calls，含有亞硫酸氫鹽轉(zhuǎn)換相關(guān)的非亞硫酸氫鹽錯配可能導(dǎo)致錯誤的比對（如果使用特別寬松的比對參數(shù)）站故。

（6）Methylation call

看一下上一步輸出的bam文件：

$ samtools view test_data_bismark_bt2.bam | head -n 5
SRR020138.15024317_SALK_2029:7:100:1672:902_length=86   16      chr1    57333005        42      50M     *       0   0         TTCTTTCCCATCCCATAAATCCTAAAAATAATAAAAAATCATCCCCAAAT       @@:AC@<=+@?+8)@BCCCA=6BCCCCCCCCCCCCCCCCACB=<88BCCA      NM:i:11       MD:Z:14G2G6G0G0G0G4G1G1G0G10G1   XM:Z:..............z..h......hhhh....h.h.hh..........h. XR:Z:CT XG:Z:GA
SRR020138.15024318_SALK_2029:7:100:1672:137_length=86   0       chr12   129289551       8       50M     *       0   0         AAAAAAAAAAAAAAGAAAAAAAAGAAAAAGAAAAGGAAAAGTAAAAAAAA       =@CAA=@B@CB=98%:AB?>@56/=3<=<)>B@:*=:=61%,<A@@1+12      NM:i:2        MD:Z:41C5G2      XM:Z:.........................................h........ XR:Z:CT XG:Z:CT
SRR020138.15024319_SALK_2029:7:100:1672:31_length=86    0       chr2    10026448        42      50M     *       0   0         ATTTTGTTATAGAGTGGGGTATTTTCGGGAAGAAGGAGGAGGAGTGTATT       BCCCCBCCCCA?:=ACCBCABCCCCCBCCA??5=9@4BB@;??B@BABBA      NM:i:8        MD:Z:1C1C5C9C2C0C22C1C1  XM:Z:.h.x.....x.........h..hh.Z....................h.x. XR:Z:CT XG:Z:CT
SRR020138.15024320_SALK_2029:7:100:1672:1164_length=86  16      chr5    28344365        8       50M     *       0   0         CACAAAATATCAACACCCCTAAACCCCACATTATTCAAAAATCAATTATA       @@@BBBA@A9=A@<?::2:<CB@?=:BBAC??CB@@BBBBC>:ACABCAB      NM:i:11       MD:Z:4G1G1G3G9G9G5G0G0G1T4G2     XM:Z:....x.h.h...x.........h.........h.....hhh......h.. XR:Z:CT XG:Z:GA
SRR020138.15024321_SALK_2029:7:100:1672:433_length=86   0       chr14   38241894        42      50M     *       0   0         TTTTGAGTAGAGAAGTTAGTATTTTAGGGAATTTTTGATTTTTTTAAGTT       BCCBB?B@@A>@-4BBB:7@BBBCBBC@@=A@BCACA;BCBBCBB@@@BB      NM:i:14       MD:Z:0C0C13C0C6C0C6C0C1C3C0C1C0C1C5      XM:Z:hh.............hx......hx......hh.x...hh.hh.h..... XR:Z:CT XG:Z:CT

甲基化call字符串對于BS-read中不涉及胞嘧啶的每個位置都用一個點"."代替皆怕，或包含以下三個不同胞嘧啶甲基化的字母(大寫=甲基化，小寫=未甲基化):

z - C in CpG context - unmethylated #C在CpG里西篓，非甲基化
Z - C in CpG context - methylated #C在CpG里愈腾，甲基化
x - C in CHG context - unmethylated #C在CHG里，非甲基化
X - C in CHG context - methylated #C在CHG里岂津，甲基化
h - C in CHH context - unmethylated #C在CHH里虱黄，非甲基化
H - C in CHH context - methylated #C在CHH里，甲基化
u - C in Unknown context (CN or CHN) - unmethylated #C在未知區(qū)域里吮成，非甲基化
U - C in Unknown context (CN or CHN) - methylated #C在未知區(qū)域里橱乱，甲基化
. - not a C or irrelevant position #該位置不是胞嘧啶

所以根據(jù)這個信息我們就可以知道每一條read的甲基化信息了辜梳。

（7）Bowtie 2 或HISAT2模式的局部比對

(在上一個版本中提過:https://github.com/FelixKrueger/Bismark/releases/tag/0.22.0)

Wu等人發(fā)表的一篇文章進(jìn)一步闡述了我們對單細(xì)胞的BS-seq(通常稱為PBAT文庫)可以生成嵌合reads（generate chimeric read pairs），該文章進(jìn)一步詳細(xì)描述了片段內(nèi)嵌合體阻礙單細(xì)胞BS-seq文庫的高效比對仅醇。在那篇文章里冗美，作者描述了一個pipeline，首先使用paired-end比對析二，然后是單端比對步驟粉洼，使用局部比對來提高分子內(nèi)嵌合體的mapping。為了允許對單細(xì)胞或PBAT庫進(jìn)行這種類型的改進(jìn)叶摄，Bismark還允許使用局部對齊属韧。

請注意，我們?nèi)匀徊唤ㄗh使用局部比對來大幅增加mapping效率(請參見here)蛤吓，但是我們承認(rèn)PBAT / scBSs-seq /scNMT-seq特殊應(yīng)用局部比對可能的確有所不同宵喂。我們還沒為局部比對設(shè)置更合適或更嚴(yán)格的默認(rèn)值，也沒有調(diào)查是否在甲基化提取步驟中需要--ignore標(biāo)簽会傲。

(III) Bismark去重步驟

基本格式：

deduplicate_bismark --bam [options] <filenames>

deduplicate_bismark腳本是為了從Bismark輸出中(包括單端和雙端的SAM/BAM文件)去除基因組中同一位置的比對锅棕，這些重復(fù)由于PCR擴(kuò)增過度而產(chǎn)生。排列在同一基因組位置但在不同鏈上的序列是單獨計分的淌山。

值得注意的是裸燎，RRBS不推薦去重、擴(kuò)增子或其他靶標(biāo)類型豐富的文庫也不建議去重泼疑。

在默認(rèn)模式下德绿，與給定位置的第一次比對將被使用，而不管它的甲基化call退渗。由于比對沒有以任何方式進(jìn)行排序移稳，所以對于每個位置來說，這已經(jīng)足夠接近隨機(jī)讀取了会油。

對于單端比對个粱，只有染色體、起始坐標(biāo)和read鏈將用于去重翻翩。

對于雙端比對几蜻，染色體、read鏈体斩、第一條read的開始坐標(biāo)以及第二條read的開始坐標(biāo)將用于去重。該腳本期望Bismark輸出為SAM/BAM格式颖低。

請注意絮吵，對于paired-end的BAM文件，去重腳本要求Read 1和Read 2在連續(xù)的行中相互跟隨!如果文件是按位置排序的忱屑，請確保根據(jù)read名稱重新排序(例如使用samtools sort -n)蹬敲。

代碼：

$ deduplicate_bismark --bam test_data_bismark_bt2.bam

生成文件：

（1）使用UMI或者barcodes去重

在去重的時候暇昂，除了染色體，起始位置(和結(jié)束位置)和鏈的方向的選項--barcode伴嗡，還要考慮潛在的barcodes或UMIs急波。barcode需要是read ID的最后一個元素，并且必須由一個:分隔瘪校，比如：

MISEQ:14:000000000-A55D0:1:1101:18024:2858_1:N:0:CTCCT

（2）對相同文庫的多個文件進(jìn)行去重

當(dāng)使用選項--multiple時澄暮，所有指定的輸入文件都被視為單個樣本，并連接在一起進(jìn)行去重阱扬。它對SAM文件使用Unix cat泣懊，對BAM文件使用samtools cat。
BAM文件的附加說明:盡管這可以在BAM或CRAM上工作麻惶，但所有輸入文件必須彼此具有相同的格式馍刮。每個輸入文件的序列字典必須是相同的。默認(rèn)情況下窃蹋，header取自要連接的第一個文件卡啰。

下一篇將進(jìn)行甲基化信息提取，以及輸出文件的解讀警没。