0.質(zhì)控
trim_galore -q 20 --phred33 --illumina --stringency 1 -e 0.1 \
--length 36 --rrbs --cores 4 --fastqc -o ./cut_adapter \
--paired R1.fastq R2.fastq
-q 設(shè)定phred quality閾值，默認(rèn)20（99%的read質(zhì)量）门岔，如果測(cè)序深度較深，可以設(shè)定25
--phred33 默認(rèn)記分方式烤送，代表Q+33=ASCII碼的方式來(lái)記分方式
--stringency 設(shè)定可以忍受的前后adapter重疊的堿基數(shù)寒随，默認(rèn)是1
--illumina 指定要去除的接頭序列，可先前運(yùn)行一下fastqc來(lái)確定接頭序列類型
--length 設(shè)定長(zhǎng)度閾值帮坚，小于此長(zhǎng)度會(huì)被拋棄
--rrbs 指定輸入文件是MspI消化的RRBS樣本（識(shí)別位點(diǎn)：）CCGG
--fastqc 修剪完成后妻往，以默認(rèn)模式在FastQ文件上運(yùn)行FastQC
-o 輸出目錄
-e 設(shè)定默認(rèn)質(zhì)量控制數(shù)试和，默認(rèn)是0.1讯泣，即ERROR rate大于10%的read 會(huì)被舍棄，如果添加了--paired參數(shù)則會(huì)舍棄一對(duì)reads
--paired 對(duì)于雙端結(jié)果阅悍，一對(duì)reads中若一個(gè)read因?yàn)橘|(zhì)量或其他原因被拋棄晦墙，則對(duì)應(yīng)的另一個(gè)read也拋棄寡痰。
<filename> #如果是采用illumina雙端測(cè)序的測(cè)序文件连躏，應(yīng)該同時(shí)輸入兩個(gè)文件入热。
1.建立亞硫酸鹽基因組index
bismark_genome_preparation --bowtie2 基因組所在目錄（目錄下應(yīng)含有fa或fasta格式的文件）
2. 比對(duì)
bismark --bowtie2 -p 4 -N 0 -L 20 --quiet --un --ambiguous 亞硫酸鹽基因組index所在目錄 -1 fastq1 -2 fastq2 -o 輸出目錄
-p 線程數(shù)
-N 允許錯(cuò)配數(shù)目
-L 比對(duì)時(shí)建立的seed大小
--quiet 不在控制臺(tái)輸出比對(duì)流程信息
--un 輸出不匹配序列信息到_unmapped_reads.fq.gz中
--ambiguous 輸出多處匹配序列信息到_ambiguous_reads.fq中，不指定則輸出到_unmapped_reads.fq.gz中
-o 輸出目錄
--sam 指定輸出的格式為sam蠢箩，無(wú)此參數(shù)則輸出bam格式文件
輸出結(jié)果有兩個(gè)重要文件逻谦，一個(gè)bismark_bt2_PE_report.txt，一個(gè)bismark_bt2_pe.bam
3-1.去重（可選勇婴，建議對(duì)WGBS使用此步驟，但不應(yīng)將其用于代表性較低的文庫(kù)础米，例如RRBS屁桑，擴(kuò)增子或靶標(biāo)富集文庫(kù)）
deduplicate_bismark -p --bam 輸入的bam文件 -o 輸出目錄
-p 代表雙端測(cè)序
--bam 代表輸入文件為bam文件
-o 輸出目錄
3-2.提取甲基化信息（可選）
bismark_methylation_extractor -p --comprehensive --no_overlap --bedGraph --counts --buffer_size 20G --report --cytosine_report --genome_folder 基因組所在目錄 -o 輸出目錄 輸入的bam文件
-s 輸入文件是從單端讀取數(shù)據(jù)生成的Bismark結(jié)果文件蘑斧。如果既沒有-s也不-p設(shè)置實(shí)驗(yàn)的類型將被自動(dòng)確定
-p 輸入文件是從雙端讀取數(shù)據(jù)生成的Bismark結(jié)果文件。如果既沒有-s也不-p設(shè)置實(shí)驗(yàn)的類型將被自動(dòng)確定
--comprehensive 將可能的甲基化信息都輸出，包括CHG,CHH,CpG
--no_overlap 對(duì)于雙端讀取庸论，read_1和read_2理論上是可能重疊的罐农。這個(gè)選項(xiàng)避免了重復(fù)的甲基化計(jì)算了2次。雖然這消除了對(duì)序列片段中心更多甲基化的計(jì)算偏差麸恍，但實(shí)際上可能會(huì)刪除相當(dāng)一部分?jǐn)?shù)據(jù)，默認(rèn)調(diào)用此參數(shù)
--bedGraph 輸出bedGraph文件
--counts 指在bedGraph中有每個(gè)C上甲基化reads和非甲基化reads的數(shù)目(在--cytosine_report指定的情況下)
--buffer_size 使用指定的內(nèi)存去計(jì)算甲基化信息
--report 會(huì)有一個(gè)甲基化的summary
--cytosine_report 指報(bào)道全基因組所有的CpG
--genome_folder 參考基因組的位置
-o 輸出目錄
重要文件有bismark_bt2_pe.CpG_report.txt（可以用來(lái)變化格式使methylkit包可以讀让羝）和bismark_bt2_pe.bedGraph.gz（可以用IGV打開，報(bào)告給定胞嘧啶的位置及其甲基化狀態(tài)）
3-3.總結(jié)甲基化信息（可選）
bismark2report（使用bismark_bt2_PE_report.txt）
bismark2summary（使用bismark_bt2_pe.bam）
4.排序
samtools sort 輸入的bam文件 -o 輸出的排序后的bam文件
(如要指定目錄的話粹淋，目錄一定要先創(chuàng)建好廓啊，其不能在運(yùn)行時(shí)創(chuàng)建)
5.讀入R
my.methRaw=processBismarkAln(location =bam文件全路徑(用list表示),
sample.id=樣品分組(用list表示),
assembly="mm10",# 比對(duì)時(shí)所用的參考基因組
read.context = "none",#要讀入內(nèi)存的內(nèi)容
save.context=c("CpG","CHG","CHH"),#要保存的內(nèi)容
save.folder=getwd())#保存路徑
6.回到terminal第步，進(jìn)入工作目錄樊展，去除非一般染色體
cat CpG文件 | grep chr > 新的CpG文件
7.重新讀入R
myobj=methRead(location=CpG文件全路徑(用list表示),
sample.id=樣品分組(用list表示),
assembly="mm10",
treatment=c(1,1,0,0),
context="CpG",
mincov = 10)
示例：
library(methylKit)
file.list=list(system.file("extdata", "test1.myCpG.txt", package = "methylKit"),
system.file("extdata", "test2.myCpG.txt", package = "methylKit"),
system.file("extdata", "control1.myCpG.txt", package = "methylKit"),
system.file("extdata", "control2.myCpG.txt", package = "methylKit"))
myobj=methRead(file.list,
sample.id=list("test1","test2","ctrl1","ctrl2"),
assembly="hg38",
treatment=c(1,1,0,0),
context="CpG",
mincov = 10)
#理想情況下墩弯，您應(yīng)該首先過濾具有極高（低）覆蓋率的堿基脓魏，可使用filterByCoverage()函數(shù)解決PCR偏倚性珊燎，
#然后運(yùn)行normalizeCoverage()函數(shù)以標(biāo)準(zhǔn)化樣本之間的覆蓋率芦瘾。
#這兩個(gè)功能將有助于減少統(tǒng)計(jì)測(cè)試中由于某些樣本讀數(shù)的系統(tǒng)過采樣而可能產(chǎn)生的偏差祷愉。
getMethylationStats(myobj[[1]],plot=TRUE,both.strands=FALSE)#樣品統(tǒng)計(jì)信息
getCoverageStats(myobj[[1]],plot=TRUE,both.strands=FALSE)#樣品統(tǒng)計(jì)信息
myobj=filterByCoverage(myobj,lo.count=10,lo.perc=NULL,hi.count=NULL,hi.perc=99.9)
myobj <- normalizeCoverage(myobj)
tiles=tileMethylCounts(myobj,win.size=500,step.size=500)
meth=unite(tiles, destrand=FALSE)#合并樣本_DMR
# meth=unite(myobj, destrand=FALSE)#合并樣本_DML
head(meth)
#樣品聚類信息
getCorrelation(meth,plot=TRUE)
clusterSamples(meth, dist="correlation", method="ward", plot=TRUE)
# hc = clusterSamples(meth, dist="correlation", method="ward", plot=F)
PCASamples(meth, screeplot=TRUE)
PCASamples(meth)
##查找差異甲基化堿基或區(qū)域
myDiff=calculateDiffMeth(meth)
myDiff25p=getMethylDiff(myDiff,difference=25,qvalue=0.01)
diffMethPerChr(myDiff,plot=T,qvalue.cutoff=0.01, meth.cutoff=25)
##統(tǒng)計(jì)甲基化位點(diǎn)的基因組分布和cpg島分布(即注釋差異甲基化的堿基或區(qū)域)
library(genomation)
myDiff_GRanges <- as(myDiff25p,"GRanges")
gene.obj=readTranscriptFeatures("c:/Users/ZSY/Desktop/DNA_methylation/hg38.bed.txt")#讀取bed文件躯嫉，有關(guān)啟動(dòng)子/外顯子/內(nèi)含子的注釋
diffAnn=annotateWithGeneParts(myDiff_GRanges,gene.obj)# 差異甲基化位點(diǎn)的基因組分布
#library(TxDb.Hsapiens.UCSC.hg38.knownGene)
#diffAnn <- annotateDMRInfo(myDiff_GRanges, 'TxDb.Hsapiens.UCSC.hg38.knownGene')
cpg.obj=readFeatureFlank("c:/Users/ZSY/Desktop/DNA_methylation/hg38.bed.txt",feature.flank.name=c("CpGi","shores"))#讀取CpG島（CpGi）和CpG海岸注釋
diffCpGann=annotateWithFeatureFlank(myDiff_GRanges,cpg.obj$CpGi,cpg.obj$shores,
feature.name="CpGi",
flank.name="shores")#差異甲基化位點(diǎn)的CpG島分布
head(getMembers(diffAnn))#獲得感興趣的區(qū)域是否與外顯子/內(nèi)含子/啟動(dòng)子重疊
head(getMembers(diffCpGann))#獲得感興趣的區(qū)域是否與CpG島重疊
head(getAssociationWithTSS(diffAnn))#獲得距TSS的距離最近的基因名稱
getTargetAnnotationStats(diffAnn,percentage=TRUE,precedence=TRUE)#與內(nèi)含子/外顯子/啟動(dòng)子重疊的差異甲基化區(qū)域的百分比/數(shù)目
plotTargetAnnotation(diffAnn,precedence=TRUE,main="differential methylation annotation")#上面的數(shù)據(jù)畫圖
getFeatsWithTargetsStats(diffAnn,percentage=TRUE)#與差異甲基化堿基重疊的內(nèi)含子/外顯子/啟動(dòng)子占所有的內(nèi)含子/外顯子/啟動(dòng)子的百分比
plotTargetAnnotation(diffCpGann,col=c("green","gray","white"),main="differential methylation annotation")#顯示CpG島，CpG海岸和其他地區(qū)差異甲基化堿基的百分比
上面所用的gene和cpg bed文件可從ucsc下載
