DSS是一個(gè)可用于做RNA-seq差異表達(dá)分析或甲基化差異分析的R包绝骚，在做差異甲基化分析時(shí)，DSS對每個(gè)CpG進(jìn)行統(tǒng)計(jì)檢驗(yàn)垦沉，然后根據(jù)我們指定的閾值可以篩選出差異甲基化位點(diǎn)（differential methylation loci, DML）或差異甲基化區(qū)域（differential methylation regions, DMR）迫皱。
DSS的輸入文件格式如下：
chr：染色體妄痪；
pos：CpG位點(diǎn)；
N：所有reads數(shù)目嘁傀；
X：檢測到甲基化的reads數(shù)目兴蒸；
使用python腳本從bismark輸出的report文件中提取上述信息并構(gòu)建符合上述格式的文件，python代碼如下：
import sys
infile = sys.argv[1]
def writelis(lis, fil):
? ? with open(fil, "w") as out_f:
? ? ? ? firstline = "chr" + "\t" + "pos" + "\t" + "N" + "\t" + "X" + "\n"
? ? ? ? out_f.write(firstline)
? ? ? ? for it in lis:
? ? ? ? ? ? line1 = it
? ? ? ? ? ? out_f.write(line1)
? ? out_f.close()
def report2dss(reportfile):
? ? with open(reportfile, "r") as ref:
? ? ? ? CG = []
? ? ? ? CHG = []
? ? ? ? CHH = []
? ? ? ? file1=reportfile.split("clean")[0] + "CGout.txt"
? ? ? ? file2 = reportfile.split("clean")[0] + "CHGout.txt"
? ? ? ? file3 = reportfile.split("clean")[0] + "CHHout.txt"
? ? ? ? f = ref.readlines()
? ? ? ? for line in f:
? ? ? ? ? ? lin = line.strip().split()
? ? ? ? ? ? chr = lin[0]
? ? ? ? ? ? pos = lin[1]
? ? ? ? ? ? type = lin[-2]
? ? ? ? ? ? numc = int(lin[3])
? ? ? ? ? ? numn = int(lin[4])
? ? ? ? ? ? allc = numn + numc
? ? ? ? ? ? DDSline = chr + "\t" + pos + "\t" + str(allc) + "\t" + str(numc) + "\n"
? ? ? ? ? ? if type == "CG":
? ? ? ? ? ? ? ? CG.append(DDSline)
? ? ? ? ? ? elif type == "CHG":
? ? ? ? ? ? ? ? CHG.append(DDSline)
? ? ? ? ? ? elif type == "CHH":
? ? ? ? ? ? ? ? CHH.append(DDSline)
? ? ? ? ? ? else:
? ? ? ? ? ? ? ? print(line)
? ? ? ? print(len(CG), len(CHG), len(CHH))
? ? writelis(CG, file1)
? ? writelis(CHG, file2)
? ? writelis(CHH, file3)
report2dss(infile)
輸出文件會(huì)將CG细办，CHG橙凳，CHH分開，并符合DSS輸入需求：
隨后參考http://www.reibang.com/p/a81c3176238b做DMR分析笑撞，代碼如下：
if (!requireNamespace("BiocManager", quietly = TRUE))
? install.packages("BiocManager")
BiocManager::install("DSS")
library(DSS)
require(bsseq)
##CG
data1.1 <- read.table("xiao-F-4-1_P_1_CGout.txt", header = T)? ##sex-reverse
data1.2 <- read.table("xiao-F-4-2_P_1_CGout.txt", header = T)? ##sex-reverse
data1.3 <- read.table("xiao-F-4-3_P_1_CGout.txt", header = T)? ##sex-reverse
data2.1 <- read.table("xiao-F-50-1_P_1_CGout.txt", header = T)
data2.2 <- read.table("xiao-F-50-2_P_1_CGout.txt", header = T)
data2.3 <- read.table("xiao-F-50-3_P_1_CGout.txt", header = T)
data3.1 <- read.table("xiao-M-44-1_P_1_CGout.txt", header = T)
data3.2 <- read.table("xiao-M-44-2_P_1_CGout.txt", header = T)
data3.3 <- read.table("xiao-M-44-3_P_1_CGout.txt", header = T)
head(data1.1)
head(data1.2)
head(data2.1)
head(data2.2)
head(data3.1)
head(data3.2)
BSobj <- makeBSseqData(list(data1.1, data1.2, data1.3, data2.1, data2.2, data2.3,
? ? ? ? ? ? ? ? ? ? ? ? ? ? data3.1, data3.2, data3.3),
? ? ? ? ? ? ? ? ? ? ? c("F4-1", "F4-2", "F4-3", "F50-1", "F50-2",
? ? ? ? ? ? ? ? ? ? ? ? "F50-3", "M44-1", "M44-2", "M44-3"))
snow <- SnowParam(workers = 9)
dmlResult <- DMLtest(BSobj, group1 = c("F4-1", "F4-2", "F4-3"),
? ? ? ? ? ? ? ? ? ? group2 = c("F50-1", "F50-2", "F50-3"),smoothing=TRUE,BPPARAM=snow)
dmlResult2 <- DMLtest(BSobj, group1 = c("F4-1", "F4-2", "F4-3"),
? ? ? ? ? ? ? ? ? ? ? group2 = c("M44-1", "M44-2", "M44-3"),smoothing=TRUE,BPPARAM=snow)
dmlResult3 <- DMLtest(BSobj, group1 = c("F50-1", "F50-2", "F50-3"),
? ? ? ? ? ? ? ? ? ? ? group2 = c("M44-1", "M44-2", "M44-3"),smoothing=TRUE,BPPARAM=snow)
##DML1
dmls <- callDML(dmlResult,delta=0.25,p.threshold = 0.01)
write.table(dmls, "Fre4_f50_CG0.25.bed", sep="\t",row.names=FALSE, quote=FALSE)
dmrs <- callDMR(dmlResult,delta=0.25,p.threshold = 0.01)
write.table(dmrs, "Fre4_f50_CG_DMR0.25.bed", sep="\t",row.names=FALSE, quote=FALSE)
##DML2
dmls <- callDML(dmlResult2,delta=0.25,p.threshold = 0.01)
write.table(dmls, "Fre4_male_CG0.25.bed", sep="\t",row.names=FALSE, quote=FALSE)
dmrs <- callDMR(dmlResult2,delta=0.25,p.threshold = 0.01)
write.table(dmrs, "Fre4_male_CG_DMR0.25.bed", sep="\t",row.names=FALSE, quote=FALSE)
##DML3
dmls <- callDML(dmlResult3,delta=0.25,p.threshold = 0.01)
write.table(dmls, "female_male_CG0.25.bed", sep="\t",row.names=FALSE, quote=FALSE)
dmrs <- callDMR(dmlResult3,delta=0.25,p.threshold = 0.01)
write.table(dmrs, "female_male_CG_DMR0.25.bed", sep="\t",row.names=FALSE, quote=FALSE)
showOneDMR(dmrs[1,], BSobj)?
