1飞袋、RRBS簡(jiǎn)介
RRBS即reduced representation bisulfite sequencing：
第一步：使用MSPI內(nèi)切酶處理碴犬，其識(shí)別序列是CCGG陕靠，切完之后，雙鏈序列如下：
CGG--------------------------------C
C--------------------------------GGC
第二步：自動(dòng)補(bǔ)齊
CGG--------------------------------CCG
GCC--------------------------------GGC
第三步：bisulfite conversion
將雙鏈中沒有甲基化的C變?yōu)閁嘉熊，甲基化的C保持不變夏块，這樣轉(zhuǎn)化之后疏咐，原來互補(bǔ)的兩條雙鏈就不互補(bǔ)了
2、數(shù)據(jù)預(yù)處理
2.1 fastqc質(zhì)量檢測(cè)
fastqc -t 2 -o ./00_fastqc R1.fastq R2.fastq
主要看Q30的百分比脐供，以及reads中有無接頭
從圖中可看出浑塞，reads中只有illumina universal adapter，因此在后續(xù)的處理中只需處理這一種接頭
2.2 trim_galore去除接頭
trim_galore --phred33 --illumina --stringency 3 -e 0.1 \
--gzip --length 35 --rrbs --fastqc -o ./cut_adapter \
--paired R1.fastq R2.fastq
詳細(xì)參數(shù)使用trim_galore -h查看政己，這里特別解釋一下--rrbs這個(gè)參數(shù)酌壕，是專門用于用MspI處理的RRBS數(shù)據(jù)，加上這個(gè)參數(shù)歇由，在切除接頭之后會(huì)自動(dòng)切除2bp的堿基卵牍，即上面說的自動(dòng)補(bǔ)齊的那兩個(gè)堿基。另外--illumina這個(gè)參數(shù)是用來切除illumina universal adapter（AGATCGGAAGAGC）的沦泌，另外參數(shù)--nextera切除nextera adapter（CTGTCTCTTATA）糊昙，--small_rna切除Illumina Small RNA 3' Adapter（GATCGTCGGACT）
3、使用bismark轉(zhuǎn)換基因組序列
參考基因組下載：
https://support.illumina.com/sequencing/sequencing_software/igenome.html
bismark_genome_preparation --path_to_bowtie dir_bowtie --genomic_composition \
--verbose dir_ref dir_hg38.fa
4谢谦、序列比對(duì)
bismark dir_ref -1 $R1 -2 $R2 --path_to_bowtie $bowtie \
-o ./01_align --rg_tag --rg_id $fname --rg_sample $fname \
--unmapped --prefix $fname --basename $fname --nucleotide_coverage
比對(duì)好的結(jié)果為bam格式的释牺，其內(nèi)容與BWA比對(duì)后的bam文件稍微有些不同
上面是一條bam文件的記錄，共17列回挽，下面分別解釋一下各列的含義：
第1列：
E00500:282:HN3M3CCXY:7:1101:8613:2628_1:N:0:ATCACG
#reads編號(hào)没咙，_1代表來自原來的R1，后面的ATCACG為index序列</pre>
第2列:
此列為一數(shù)值厅各，具體解釋如下圖
相同數(shù)值對(duì)于R1和R2來說代表不同的鏈镜撩，此外，這些數(shù)值也是由多個(gè)值相加所得队塘，下面這幾個(gè)標(biāo)簽都是2的n次方袁梗，隨機(jī)挑選其中的幾個(gè)，他們的和是唯一的：
| 1 | 代表這個(gè)序列采用的是PE雙端測(cè)序 |
| 2 | 代表這個(gè)序列和參考序列完全匹配憔古，沒有插入缺失 |
| 4 | 代表這個(gè)序列沒有mapping到參考序列上 |
| 8 | 代表這條序列的另一端序列沒有比對(duì)到參考序列上遮怜，如果這條是R1，它對(duì)應(yīng)的R2端沒有比對(duì)上 |
| 16 | 代表這個(gè)序列比對(duì)到參考序列的負(fù)鏈上 |
| 32 | 代表這個(gè)序列對(duì)應(yīng)的另一端序列比對(duì)到參考序列的負(fù)鏈上 |
| 64 | 代表這個(gè)序列是R1端序列 |
| 128 | 代表這個(gè)序列是R2端序列 |
| 256 | 代表這個(gè)序列不是主要的比對(duì)鸿市，一條序列可能比對(duì)到了多個(gè)位置锯梁，只有一個(gè)是首要的比對(duì)位置 |
| 512 | 代表這個(gè)序列在QC時(shí)失敗了，被過濾不掉了 |
| 1024 | 代表這個(gè)序列是PCR重復(fù)序列 |
| 2048 | 代表這個(gè)序列是補(bǔ)充的比對(duì)焰情，不常用 |
第3-4列：
chr4 184651014
代表這條reads比對(duì)到4號(hào)染色體的184651014這個(gè)位置
第5列：
為一數(shù)值陌凳，Mapping quality，比對(duì)的質(zhì)量分?jǐn)?shù)内舟，越高說明該read比對(duì)到參考基因組上的位置越唯一
第6列：
CIGAR值合敦，堿基匹配上的堿基數(shù)。
| M | match/mismatch | I | insert | |
| D | deletion | N | skipped（跳過這段區(qū)域） | |
| S | soft clipping（被剪切的序列存在于序列中） | H | hard clipping（被剪切的序列不存在于序列中） | |
| P | padding(填充) |
135M
代表135個(gè)堿基在比對(duì)時(shí)完全比對(duì)
3S6M1P1I4M
代表前三個(gè)堿基被剪切去除了验游，然后6個(gè)比對(duì)上了充岛，然后打開了一個(gè)缺口，有一個(gè)堿基插入耕蝉，最后是4個(gè)比對(duì)上了崔梗，是按照順序的
第7列：
MRNM(chr)，mate的reference sequence name垒在，實(shí)際上就是mate比對(duì)到的染色體號(hào)蒜魄，若是沒有mate，則是*场躯，若為=則代表與此條比對(duì)到相同的染色體上谈为。
第8列：
mate position，mate比對(duì)到參考序列上的第一個(gè)堿基位置推盛，若無mate,則為0
第9列：
insert size
第10列：
Sequence峦阁，就是read的堿基序列，如果是比對(duì)到負(fù)鏈上則對(duì)read進(jìn)行了reverse completed
第11列：
ASCII耘成，read質(zhì)量的ASCII編碼
第12列：
NM-tag榔昔，與bowtie中NM:i相同。read string轉(zhuǎn)換成reference string需要的最少核苷酸的edits:插入/缺失/替換
第13列：
XX-tag瘪菌，對(duì)錯(cuò)配的描述撒会，不包括indel(插入和缺失)
XX:Z:2C4C2C6C1C1AC18C3C15C4C2CC14
表示2個(gè)堿基完全匹配，一個(gè)C替換师妙，接著4個(gè)堿基完全匹配诵肛，一個(gè)C發(fā)生替換
第14列：
XM-tag，methylation call string
XM:Z:........xz.z....hx............x....h.hh....xz....xz...x....h.hh..hhhx....z...h......h
CHG默穴、CHH怔檩、CG代表甲基化的三種模式
CG：C后面接的是G
CHG：C后面跟了任意一個(gè)堿基褪秀，再接上一個(gè)G
CHH：C后面跟的是除了G以外的任何堿基
第15列：
XR-tag，read conversion state for the alignment 薛训。
共兩種轉(zhuǎn)換：CT和GA媒吗，GA就是指將read里的所有G轉(zhuǎn)換成A
第16列：
XG-tag，genome conversion state for the alignment乙埃。
共兩種：GA和CT闸英。CT是指將全基因組上所有的C轉(zhuǎn)換成T
第17列：
表示樣本名稱，是在做比對(duì)時(shí)通過--rg_tag自己設(shè)置的
5介袜、去除重復(fù)序列
deduplicate_bismark --paired --bam bam_file --output_dir dir_bamfile
注：RRBS其實(shí)是不需要做這一步的
6甫何、抽提出甲基化的統(tǒng)計(jì)信息
bismark_methylation_extractor --paired-end --report --output dir --gzip --bedGraph bam_file
產(chǎn)生的結(jié)果文件如下：
CHG_OB_FF01N_ATCACG_S9_L007_pe.deduplicated.txt.gz
CHG_OT_FF01N_ATCACG_S9_L007_pe.deduplicated.txt.gz
CHH_OB_FF01N_ATCACG_S9_L007_pe.deduplicated.txt.gz
CHH_OT_FF01N_ATCACG_S9_L007_pe.deduplicated.txt.gz
CpG_OB_FF01N_ATCACG_S9_L007_pe.deduplicated.txt.gz
CpG_OT_FF01N_ATCACG_S9_L007_pe.deduplicated.txt.gz
FF01N_ATCACG_S9_L007_pe.deduplicated.bedGraph.gz
FF01N_ATCACG_S9_L007_pe.deduplicated.bismark.cov.gz
FF01N_ATCACG_S9_L007_pe.deduplicated.M-bias.txt
FF01N_ATCACG_S9_L007_pe.deduplicated_splitting_report.txt
