DSS是一个可用于做RNA-seq差异表达分析或甲基化差异分析的R包，在做差异甲基化分析时，DSS对每个CpG进行统计检验，然后根据我们指定的阈值可以筛选出差异甲基化位点（differential methylation loci, DML）或差异甲基化区域（differential methylation regions, DMR）。
DSS的输入文件格式如下：
chr：染色体；
pos：CpG位点；
N：所有reads数目；
X：检测到甲基化的reads数目；
使用python脚本从bismark输出的report文件中提取上述信息并构建符合上述格式的文件，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)
输出文件会将CG，CHG，CHH分开，并符合DSS输入需求：
随后参考https://www.jianshu.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)
网友评论