甲基化转化率计算

最新推荐文章于 2022-05-21 08:27:27 发布

Keiji1102

最新推荐文章于 2022-05-21 08:27:27 发布

阅读量4.4k

点赞数

分类专栏：表观遗传

本文链接：https://blog.csdn.net/weixin_44649331/article/details/103980347

版权

6 篇文章 4 订阅

订阅专栏

BS转化后需要加入一段已知lambdaDNA序列。
比较sam文件中CT转换情况和lambdaDNA的ref原始碱基进行比较。

perl MethylExtractBSCR.pl 
		seqFile = <sequence file>
		inFile = <alignments input file>
		flagW = <Watson FLAGs>
		flagC = <Crick FLAGs>

单端测序：
0代表单端测序，16代表这个序列比对到参考序列的负链上。
MethylExtractBSCR.pl seqFile=lambdaDNA.fa inFile=input.lambda.sam flagW=0 flagC=16
双端测序：
99代表双端测序、和参考序列完全匹配、比对到正链、first in pair
147代表双端测序、和参考序列完全匹配、比对到负链、second in pair
83代表双端测序、和参考序列完全匹配、比对到负链、first in pair
163代表双端测序、和参考序列完全匹配、比对到正链、second in pair
MethylExtractBSCR.pl seqFile=lambdaDNA.fa inFile=input.lambda.sam flagW=99,147 flagC=83,163

分析用文件：bismark_cov
在这里插入图片描述
在转化率为A的基础上，N个reads中被检测到X个甲基化reads是否可靠。需要用二项分布来证明，FDR来校正。

p = binom.test(x,N,A)$p.value
fdr = p.adjust(p,"fdr")

FDR校正后，只保留FDR<0.01的甲基化位点。再过滤read数目<5的位点。

关注