WGBS上游分析

WGBS的分析全流程：

主要參考資料：
WGBS甲基化分析
Bismark軟體使用入門
沉浸式體驗WGBS(上游)
甲基化流程淺析
聽說你不會處理WGBS資料？安排上
全基因組甲基化分析簡述：使用BS-Seeker2
Bismark Bisulfite Mapper學習筆記（二）甲基化資訊提取以及檔案解讀
DNA甲基化分析流程詳解

檢視使用者quota：quota -uvs

資料清洗

Trim galore可以自動檢測adapter對NGS資料進行質量過濾，還可以檢視質量分佈圖。也就是等於fastqc+trimmomatic。trim_galore適用於所有高通量測序，包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera和smallRNA測序平臺的雙端和單端資料
Trim_galore的工作流程：
第一步 首先去除低質量鹼基，然後去除3'末端的adapter(如果沒有指定具體的adapter，程式會自動檢測前1million的序列）
第二步 對比前12-13bp的序列是否符合以下型別的adapter:

• Illumina: AGATCGGAAGAGC # 如果輸入引數 --Illumina,就會預設trimmed前13bp的adapter
• Small RNA: TGGAATTCTCGG # 同上 如果輸入引數 --small_rna,就會預設trimmed前12bp的adapter
• Nextera: CTGTCTCTTATA # 同上 如果輸入引數 --nextera,就會預設trimmed前12bp的adapter

在trim_galore執行中，會先生成trimmed檔案，隨後生成val檔案，可以看官方講解這二者的區別：Can you explain the difference between _val_1/_val_2 and _trimmed?

# 安裝Trim galore
conda install -c bioconda trim-galore -y

fastqc --noextract -t 10 -f fastq -o ./fastqc ./rawdata/*_1.fastq.gz
fastqc --noextract -t 10 -f fastq -o ./fastqc ./rawdata/*_2.fastq.gz

# 雙端測序資料
cleandata=./1_trim_galore
rawdata=./0_rawdata
cat sample.txt | while read id
do
echo "trim_galore --fastqc --paired -o ./1_trim_galore ./0_rawdata/${id}_1.fastq.gz ./0_rawdata/${id}_2.fastq.gz" 
done > trim_galore.sh

引數詳解：

• --phred33：指定測序資料的Phred質量值，33或者64，預設為33。如何判斷呢，可以見下面兩個文章，但是最近的資料基本上都是33了：

fastq格式檔案及phred33的判斷
Fastq 格式說明 & (Phred33 or Phred64)

• -j 8：指定執行緒數，8為最大。它的使用基於python3版本，如果報錯，可以嘗試更新python環境
• -q 20：指定質量值，切除質量得分低於20的序列
• --length：設定輸出reads長度閾值小於設定值會被拋棄，預設值20bp；即小於20bp的被去除。注意，在pe150下，不要設定太高，可以50或36
• --stringency 3：可以忍受的前後adapter重疊的鹼基數，預設為1（非常苛刻）。可以適度放寬，因為後一個adapter幾乎不可能被測序儀讀到
• --fastqc：當分析結束後，使用預設選項對結果檔案進行fastqc分析
• --paired：對於雙端測序結果，一對reads中，如果有一個被剔除，那麼另一個會被同樣拋棄，而不管是否達到標準
• --max_n 3：read中允許包含的N的總數（整數），否則將完全刪除該read。 在配對末端設定中，任一read超過此限制將導致整個配對從修剪後的輸出檔案中刪除。
• -o：指定輸出目錄

資料質控

# 安裝MultiQc
conda install -c bioconda multiqc -y
# 整合FastQC結果
multiqc *.zip

WGBS比對

準備參考基因組

# 下載NCBI上豬的參考基因組
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/GCF_000003025.6_Sscrofa11.1_genomic.fna.gz
# 解壓並重新命名
gunzip GCF_000003025.6_Sscrofa11.1_genomic.fna.gz && mv GCF_000003025.6_Sscrofa11.1_genomic.fna pig.fa

# 下載Ensembl上豬的參考基因組
wget ftp://ftp.ensembl.org/pub/release-102/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz

# 解壓並重新命名
gunzip Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz && mv Sus_scrofa.Sscrofa11.1.dna.toplevel.fa pig.fa

探索參考基因組是否有NW_018085246.1

# 檢視參考基因組的染色體編號
grep ">" pig.fa | cut -d " " -f 1 | cut -d ">" -f 2 | sort | uniq -c > pig.fa.chr.txt

# 檢視包含NW_018085246對應的行數
grep -n NW_018085246 pig.fa
# 31088069:>NW_018085246.1 Sus scrofa isolate TJ Tabasco breed Duroc unplaced genomic scaffold, Sscrofa11.1 Contig556, whole genome shotgun sequence

# 輸出31088069-31088070行
sed -n '31088069,31088070p' pig.fa

轉化參考基因組，並建立索引

在比對之前，需要加入lambda序列資料，lambda序列是在重硫酸鹽實驗前人為加入的一段序列用來評估重硫酸鹽的轉化效率，這個一般由為你測序的公司提供，所以我們這一步就將其合併到參考基因組中

看到了有這種說法，但是我暫時並不知道這個lambda序列

# 安裝bismark和bowtie2
conda install -c bioconda bismark bowtie2 -y

# 檢視bowtie2的路徑
which bowtie2

# 轉化參考基因組
nohup bismark_genome_preparation \
    --bowtie2 \
    --parallel 10 \
    /data/wangji/wgbs/genome > bismark_genome_preparation.log 2>&1 &

重要引數說明：
--bowtie2：呼叫bowtie2/hisat2建立基因組索引
--path_to_aligner：比對軟體所在資料夾的全路徑
--verbose：輸出詳情
--parallel：設定執行緒，索引建立是並行執行，因此實際執行緒要×2
--large-index：大基因組索引建立
--yes：如果有安全類問題則自動選擇yes，比如覆蓋某個已存在的檔案
<path_to_genome_folder>：基因組所在資料夾路徑，即/data2/zhangzipeng/wang/wgbs/genome

比對

比對的中間檔案詳解見：bismark 識別甲基化位點-比對篇

cat sample.txt | while read id
do
echo "bismark --bowtie2 -p 8 -o ./2_bamfile --genome /work/home/wklgroup01/wgbs/genome -1 ./1_trim_galore/${id}_1_val_1.fq.gz -2 ./1_trim_galore/${id}_2_val_2.fq.gz &"
done > bismark.sh

nohup sh bismark.sh &

bismark -N 0 -L 20 \
    --un --ambiguous --sam \
    --bowtie2 \
    --path_to_bowtie2 ~/miniconda3/envs/wgbs/bin/ \
    --parallel 10 \
    -o ./bamfile \
    --fastq \
    --prefix SRR6457892 \
    --genome /data/wangji/wgbs/smk/pig_genome \
    -1 ./SRR6457892_1_val_1.fq.gz \
    -2 /data2/zhangzipeng/wang/wgbs/cleandata/SRR6457892_2_val_2.fq.gz

nohup bismark --bowtie2
--multicore 20
--gzip
-p 10
-o bamfile_ensembl
--genome ./genome/ncbi
-1 ./trim_reads/SRR5582006_1_val_1.fq.gz
-2 ./trim_reads/SRR5582006_2_val_2.fq.gz > bismark_ensembl.log 2>&1 &
重要引數說明：

--multicore ：這個決定了以多少個執行緒同時執行，執行的時候會產生相同數量的temp檔案，temp檔案會有正常的和C_to_T轉換的
--gzip：輸出檔案壓縮
-p:這個引數是bowtie2的引數，指定執行緒數

並行的一些方法：

cat srr | while read id
do
echo "bismark --bowtie2 --path_to_bowtie2 ~/miniconda3/envs/wgbs/bin/ -o ./mapping --fastq --prefix ${id} --genome /data2/zhangzipeng/wang/wgbs/genome -1 ${id}_1.fastq.gz -2 ${id}_2.fastq.gz"
done > bismark.sh

cat srr | while read id; do
  (bismark --bowtie2 --path_to_bowtie2 ~/miniconda3/envs/wgbs/bin/ -o ./mapping --fastq --prefix ${id} --genome /data2/zhangzipeng/wang/wgbs/genome -1 ${id}_1.fastq.gz -2 ${id}_2.fastq.gz) &
done
wait

檢視bam檔案

# 檢視bam檔案
samtools view -h SRR5582006_1_val_1_bismark_bt2_pe.bam | head -n 500
# 能夠看到有@SQ	SN:NW_018085246.1	LN:362670
# 其中SN表示染色體編號，LN表示染色體長度

輸出檔案：

(a) test.file.R1_bismark_bt2_pe.sam 所有對齊和甲基化的資訊，也即儲存了成功比對的reads資訊

甲基化呼叫字串包含一個點“.”代表 BS-read 中不涉及胞嘧啶的每個位置，或者包含以下三個不同胞嘧啶甲基化上下文的字母之一（大寫 = 甲基化，小寫 = 未甲基化）：

字母	含義
z	unmethylated C in CpG context
Z	methylated C in CpG context
x	unmethylated C in CHG context
X	methylated C in CHG context
h	unmethylated C in CHH context
H	methylated C in CHH context
u	unmethylated C in Unknown context (CN or CHN)
U	methylated C in Unknown context (CN or CHN)
對於CpG, 採用字母Z的大小寫來表徵甲基化狀態；對於CHG, 採用字母X的大小寫來表徵甲基化狀態；對於CHH， 採用字母H 的大小寫來表徵甲基化狀態。

(b) test.file.R1_bismark_bt2_PE_report.txt 對齊和甲基化的主要資訊概括

生成報告

執行bismark2report生成Testpaired_PE_report.html報告

bismark2report

比對後的資料處理

去除重複序列

RRBS一定不要做這一步

cat sample.txt | while read id
do
echo "deduplicate_bismark --bam -p ./2_bamfile/${id}_1_val_1_bismark_bt2_pe.bam --output_dir ./3_duplicate/ &"
done > bismark.sh

nohup deduplicate_bismark -p --bam SRR5582006_1_val_1_bismark_bt2_pe.bam --output_dir ./deduplicate/ > deduplicate.log 2>&1 &

重要引數說明：
--sam/--bam：刪除 Bismark 對齊產生的SAM/BAM 檔案中的重複資料
-p/--paired：前一步雙端資料產生的結果檔案
-s/--single：前一步單端資料產生的結果檔案
--samtools_path：samtools所在資料夾的全路徑
--output_dir：輸出資料夾路徑
--multiple：指定輸入檔案都作為一個樣本處理，連線在一起進行重複資料刪除

得到輸出檔案：
(a) test.file.R1_bismark_bt2_pe.deduplicated.bam
(b) test.file.R1_bismark_bt2_pe.deduplication_report.txt：去重複後的結果描述檔案

甲基化資訊提取

執行命令

# 後續的plot檔案需要Perl模組GD::Graph安裝
# conda install -c bioconda perl-gdgraph -y
bismark_methylation_extractor -p --gzip --no_overlap --comprehensive --bedGraph --count --CX_context --cytosine_report --buffer_size 20G --genome_folder /work/home/wklgroup01/wgbs/genome/ -o testfile_report ~/bismark_example/03duplicate/test.file.R1_bismark_bt2_pe.deduplicated.bam &

bismark_methylation_extractor --pair-end --comprehensive --no-overlap --bedGraph --counts --report --cytosine_report --genome_folder /work/home/wklgroup01/wgbs/genome ../3_duplicate/*.bam -o ../4_methylation_result

ls *.sam | awk '{print "\"" $0 "\""}' | paste -sd ","
# Define an array of input files
files=(
"SRR25022243_1_val_1_bismark_bt2_pe.deduplicated.sam"
"SRR15132626_1_val_1_bismark_bt2_pe.deduplicated.sam"
)

# Get the input file corresponding to the current job array index
input_file=${files[$SLURM_ARRAY_TASK_ID-1]}

# Run the bismark_methylation_extractor command for the input file
bismark_methylation_extractor \
    -p --gzip --no_overlap \
    --comprehensive \
    --bedGraph \
    --count \
    --CX_context \
    --cytosine_report \
    --buffer_size 20G \
    --genome_folder /work/home/wklgroup01/wgbs/genome/ \
    -o testfile_report \
    "/work/home/wklgroup01/wgbs/deduplicate/$input_file"

nohup bismark_methylation_extractor -p --ignore_r2 6 --multicore 10 --bedgraph -o ./methylation --cytosine_report --genome_folder ./genome/ncbi ./bamfile/*.deduplicated.bam > bismark_methylation_extractor.log 2>&1 &
重要引數說明：
--comprehensive ：合併所有四個可能的特定鏈，將甲基化資訊轉換為context-dependent的輸出檔案
--no_overlap：避免因雙端讀取read_1和read_2理論上的重疊，導致甲基化重複計算。可能會刪去相當大部分的資料，對於雙端資料的處理，預設情況下此選項處於啟用狀態，可以使用--include_overlap禁用。
-p/--paired-end：前一步雙端資料產生的結果檔案
--bedGraph：指將產生一個BedGraph檔案儲存CpG的甲基化資訊（不與--CX聯用，僅產生CpG資訊）
--CX/--CX_context：與--bedGraph聯用，產生一個包含三種型別甲基化資訊的BedGraph檔案;與--cytosine_report聯用，產生一個包含三種型別甲基化資訊的全基因組甲基化報告
--cytosine_report：指將產生一個儲存CpG的甲基化資訊的報告（不與--CX聯用時，僅產生CpG資訊），預設情況座標系基於 1
--buffer_size：甲基化資訊進行排序時指定記憶體，預設為2G
--zero_based：使用基於 0 的基因組座標而不是基於 1 的座標，預設值：OFF
--split_by_chromosome：分染色體輸出
--genome_folder：包含未修改的參考基因組和bismark_genome_preparation建立的子資料夾(CT_conversion/和GA_conversion/）的資料夾的路徑，即~/bismark_example/01index/
:BAM 格式的 Bismark 結果檔案

輸出檔案解讀：

1、CHG/CHH/CpG_context_test.file.R1_bismark_bt2_pe.deduplicated.txt

OT – original top strand原始top鏈
CTOT – complementary to original top strand原始top鏈的互補鏈
OB – original bottom strand原始bottom鏈
CTOB – complementary to original bottom strand原始bottom鏈互補鏈

來自OT和CTOT的甲基化calls將提供原始top鏈上胞嘧啶甲基化位置的資訊，來自OB和CTOB的甲基化calls將提供原始bottom鏈上胞嘧啶甲基化位置的資訊。請注意，在Bismark比對步驟中指定--directional(預設模式)選項不會對CTOT或CTOB鏈生成任何報告，所以往往只得到OT和OB的甲基化報告資訊。

col1 : 比對上的序列ID
col2 : 基因組正負鏈：+ -
col3 : 染色體編號
col4 : 染色體位置
col5 : 甲基化C的狀態（XxHhZzUu）

X 代表CHG中甲基化的C
x  代筆CHG中非甲基化的C
H 代表CHH中甲基化的C
h  代表CHH中非甲基化的C
Z  代表CpG中甲基化的C
z  代表CpG中非甲基化的C
U 代表其他情況的甲基化C(CN或者CHN)
u  代表其他情況的非甲基化C (CN或者CHN)

CpG：甲基化C下游是個G鹼基
CHH：甲基化C下游的2個鹼基都是H（A、C、T）
CHG：甲基化的C下游的2個鹼基是H和G

# 輸出NW_018085246的甲基化資訊
grep -n NW_018085246 CHG_OB_SRR5582006_1_val_1_bismark_bt2_pe.deduplicated.txt > NW_018085246.txt
# 18395:SRR5582006.5464_FCC6MP9ACXX:7:1101:20571:36789#/1	-	NW_018085246.1	42109	x
# 18404:SRR5582006.5464_FCC6MP9ACXX:7:1101:20571:36789#/1	+	NW_018085246.1	42081	X

# 檢視18395-18404行
sed -n '18395,18404p' CHG_OB_SRR5582006_1_val_1_bismark_bt2_pe.deduplicated.txt

2、test.file.R1_bismark_bt2_pe.deduplicated.bedGraph.gz

# 使用zcat進行檢視
zcat test.file.R1_bismark_bt2_pe.deduplicated.bedGraph.gz | head -n 10

col1 : 染色體編號
col2 : 胞嘧啶（c）起始位置資訊
col3 : 胞嘧啶（c）終止位置資訊
col4 : 甲基化率

3、test.file.R1_bismark_bt2_pe.deduplicated.bismark.cov.gz

由於甲基化百分數本身並不能提供在某一位置檢測到的甲基化或非甲基化read的實際read覆蓋率，bismark2bedGraph還會輸出一個coverage檔案(使用基於1的基因組座標)，該檔案具有兩個附加列:

• col5 : 甲基化數目
• col6 : 非甲基化數目

# 使用zcat進行檢視
zcat test.file.R1_bismark_bt2_pe.deduplicated.bismark.cov.gz | head -n 10

# 使用解壓後SRR5582006_1_val_1_bismark_bt2_pe.deduplicated.bismark.cov檔案輸出NW_018085246的甲基化資訊
grep -n NW_018085246 SRR5582006_1_val_1_bismark_bt2_pe.deduplicated.bismark.cov > NW_018085246.cov.txt
# head NW_018085246.cov.txt
# 48814275:NW_018085246.1	749	749	0	0	1
# 48814276:NW_018085246.1	802	802	100	1	0
# 48814277:NW_018085246.1	821	821	100	1	0

# 簡單處理檔案
## 依據冒號分割，取第二列
cut -d ":" -f 2 NW_018085246.cov.txt > NW_018085246.cov.txt.1 && mv NW_018085246.cov.txt.1 NW_018085246.cov.txt
## 增加行名：Chromosome	Start	End	Feature	methy_level

4、test.file.R1_bismark_bt2_pe.deduplicated.CX_report.txt

預設情況下只考慮CpG背景中的胞嘧啶，但是可以透過使用--CX 擴充套件到任何序列背景中的胞嘧啶

col1 : 染色體編號
col2 : 位置
col3 : 正負鏈資訊
col4 : 甲基化鹼基數目
col5 : 非甲基化鹼基數目
col6 : 型別
col7 : 三核苷酸背景

5、test.file.R1_bismark_bt2_pe.deduplicated.M-bias.txt/...R1.png/...R2.png

#reads中每個可能位置的甲基化比例
CpG context (R1)
col1 : read position
col2 : 甲基化count（count methylated）
col3 : 非甲基化count（count unmethylated）
col4 : 甲基化率（beta）
col5 : 總coverage

6、test.file.R1_bismark_bt2_pe.deduplicated_splitting_report.txt：甲基化總體報告

使用IGV檢視比對結果

參考資料：

乾貨分享 | IGV軟體的安裝和常用操作介紹
常見設定
保姆級 IGV 基因組瀏覽器使用指南（圖文詳解）
使用bismark進行甲基化分析並使用IGV對結果進行視覺化

有非常多種針對甲基化結果進行處理，以實現IGV視覺化的方法，這主要是因為IGV支援多種檔案格式，包括：

• 複製數檔案：seg
• 序列比對檔案：bam、cram
• 基因組註釋檔案；bed、gtf、gff3、psl、bigbed
• 定量檔案：wig、bedgraph、bigwig、tdf
  要想自己生成對應的檔案格式，可以參考此文件：File Formats

首先對IGV軟體匯入參考基因組及基因註釋檔案

構建索引:IGV匯入參考基因組和基因註釋檔案都需要先構建索引

igv裡面對參考基因組進行構建索引

IGV 工具欄，tools-Run igvtools；選擇index
igv裡面給註釋檔案排序，構建索引
IGV 工具欄，tools-Run igvtools；選擇sort；輸入註釋檔案；生成sort；接著，IGV 工具欄，tools-Run igvtools；選擇index；輸入剛剛的sort檔案；生成index

匯入檔案

參考基因組

IGV 工具欄，Genomes → load genome from
或者Genomes → Create genome File

gff/gtf註釋檔案
File → Load from File→找到註釋檔案即可

IGV中gtf檔案的各種樣式含義：

對bam進行視覺化

對bam檔案進行視覺化，需要先對bam檔案進行排序，然後進行索引才可以匯入IGV中，這個過程可以使用samtools進行，也可以使用IGV自帶的工具進行

幾乎所有的reads中都有很多紅色和綠色的部分，這個是因為預設的顏色模式下，未發生甲基化的C在正義鏈和反義鏈中會被轉化為T和A鹼基，被誤判為突變鹼基。可以在BAM檔案右鍵選擇"Color alignments by"→ "bisulfite mode" → "CG"，轉換為甲基化模式。這裡考慮到大多數情況下關注CG的甲基化，因此選擇了CG模式，如果關注其他甲基化，可以將"CG"替換為對應的甲基化模式。

將bedGraph檔案轉化為tdf格式並匯入

甲基化提取的輸出可以使用選項--bedGraph轉換為一個bedGraph和cov檔案。這個步驟也可以使用獨立指令碼bismark2bedGraph從甲基化提取輸出中完成，會獲得兩個檔案：bedGraph和cov
BedGraph檔案記錄各個位點的甲基化率。BedGraph檔案會比較大，不太適合直接匯入IGV中，可以：

• 按照染色體分別提取bedGraph

# 例如直接提取SRR15132626_1_val_1_bismark_bt2_pe.deduplicated.bedGraph檔案中NW_018085246.1
grep "NW_018085246.1" SRR15132626_1_val_1_bismark_bt2_pe.deduplicated.bedGraph > NW_018085246.1.bedGraph
# 如果是gz檔案，可以使用
zcat SRR15132626_1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz | grep "NW_018085246.1" > NW_018085246.1.bedGraph

# 讓目錄下的所有bedGraph.gz檔案只保留NW_018085246.1染色體的資訊
ls *.bedGraph.gz | while read id
do
echo "zcat ${id} | grep \"NW_018085246.1\" > NW_018085246.1.bedGraph.${id}"
done > NW_018085246.1.bedGraph.sh

# 讓目錄下的所有NW_018085246.1檔案只保留第二列大於139027小於144527的行
ls *.NW_018085246.1.bedGraph | while read id
do
echo "awk '{if(\$2>139027 && \$2<144527) print \$0}' ${id} > 1${id} & rm ${id}"
done > NW_018085246.1.bedGraph.sh

- 使用igvtools轉化為tdf格式
>選單欄中選擇"Tools" → "Run igvtools"，在彈出視窗中，"Command"選擇"To TDF"，"Input File"選擇需要轉換的BedGraph檔案，其他引數保持預設，點選"Run"可以進行轉化

### 將bedGraph檔案轉化為bw格式
[滑動視窗展示甲基化程度](https://www.jianshu.com/p/4533d37a605e)
```bash
# 排序
awk '{ print $1"\t"$2"\t"$3"\tid-"NR"\t"$4; }' GSM1032058_1_val_1_bismark_bt2_pe.deduplicated.bedGraph | sort-bed - > input.bed

# 下載各個染色體的長度檔案並排序
fetchChromSizes hg38 > hg38.bounds.unsorted.txt
awk '{ print $1"\t0\t"$2; }' hg38.bounds.unsorted.txt | sort-bed - > hg38.bounds.bed

# 獲取最終的bed檔案
bedops --chop 1000  hg38.bounds.bed | bedmap --faster --echo --mean --delim "\t" --skip-unmapped - input.bed > final_no_stagger.bed

# 轉換為bw檔案
bedGraphToBigWig final_no_stagger.bed hg38.bounds.unsorted.txt final.bw

最後將final_no_stagger.bed和 final.bw匯入IGV檢視

針對deduplicated.bismark.cov檔案進行處理，獲取格式

deduplicated.bismark.cov檔案的第一列到第三列是甲基化位點的座標，第四列是位點甲基化水平，第五列是methylated counts, 第六列是unmethylated counts。有了它，我們就可以知道全基因組CHG或CpG或CHH甲基化水平的分佈

## 下載包 
install.packages("BiocManager")
library(BiocManager)
BiocManager::install("tidyverse")
library(tidyverse)
## 檔案匯入
SRR5582006 <- read.delim("SRR25022243_1_val_1_bismark_bt2_pe.deduplicated.bismark.cov", header=FALSE)
## 檔案處理
SRR5582006<-SRR5582006[,1:4]
SRR5582006$Feature<-"SRR5582006"
SRR5582006<-select(SRR5582006, V1, V2, V3, Feature, V4)
SRR5582006$V4<-SRR5582006$V4/100
colnames(SRR5582006)<-c("Chromosome", "Start", "End", "Feature", "methy_level")
## 檔案匯出
write.table(SRR5582006, file="SRR5582006.igv", sep = "\t", row.names = F, quote = F)

載入得到的.igv檔案：
“File”--"Load from File"--選擇.igv檔案。如果.igv檔案超過240 MB會出現以下提示，但是我們忽略它，選擇"Continue"

針對deduplicated.CX_report.txt檔案進行處理，獲取格式

bismark_methylation_extractor會生成一個總的bedgraph檔案：*_P_1_bismark_bt2_pe.deduplicated.bedGraph.gz。但該檔案是整個基因組所有甲基化C位點，沒有分CG/CHG/CHH資訊，也不利於後續視覺化，因此根據下面程式碼可以獲取的CX_report.txt（按照染色體拆分）檔案拆分CG，CHG，CHH三類甲基化後轉換成bedgraph格式檔案。

nohup bismark_methylation_extractor -p --parallel 40 --comprehensive --no_overlap--bedGraph --cytosine_report --CX_context --split_by_chromosome --counts--buffer_size 20G --report --samtools_path yourpath/smrtlink_v9/smrtcmds/bin --genome_folder refdir bam    -o outdir &

CX_report.txt各列分別為：
col1：染色體；
col2：染色體上具體位點（以1開始）；
col3：正負鏈；
col4：比對到該位點發生甲基化的reads數；
col5：比對到該位點未發生甲基化的read數；
col6：甲基化型別，有三種：CG, CHG, CHH；
col7：該位點的三核苷酸序列；

CX_report.txt需要使用bismarkCXmethykit.pl指令碼進行轉換，轉換後出現以CG_methykit.txt，CHG_methykit.txt，CHH_methykit.txt為字尾的拆分甲基化型別的檔案，檔案格式如下：

再使用R指令碼轉換成bedgraph檔案，只需保留chr列，base列和freqC列，R指令碼如下：

Args <- commandArgs(T)
Args[1]
myfile=read.table(Args[1],header=T)
outfile=sub("_P_1_bismark_bt2_pe.deduplicated.CX_report.txt.chrChr19.CX_report.txt", "",Args[1])
file_out=sub("_methykit.txt", ".bedgraph", outfile )   ###輸出檔案
chr<-as.character(myfile$chr)
pos<-as.integer(myfile$base)
pos2<-as.integer(myfile$base)
freq<-as.numeric(myfile$freqC)
outf<-data.frame(chr,pos,pos2,freq, check.names = F)   ###合併以上四列，check.names = F可以防止將字串轉化為其他內容。
outf$chr<-factor(outf$chr)
dim(myfile)
dim(outf)  ###檢查一下新檔案和原檔案行數是否相等
write.table(outf, file=file_out, quote=F,col.names = F, row.names = F)  ##輸出

轉換後生成以CG.bedgraph/CHG.bedgraph/CHH.bedgraph為字尾的檔案，檔案內容如下：

符合IGV需求的格式，可以用於視覺化