国外的BSseq数据分析教程


Epigenomics Workshop 2021

Tutorials（目录）

Babraham Institute training Courses

开发了bismark、trim-galore等知名生信软件的机构
Babraham Institute： https://www.bioinformatics.babraham.ac.uk/training.html
注：该网站不止包含BSseq数据分析教程，还包括一堆的生信相关教程，不过课程开设的年份有点早，有些软件可能已有更好的替代品。
image.png
截图上的资料我下载了，现传上来并简单分类一下：

  • BS-Seq theory and QC lecture，BS-seq基础知识

BS-Seq theory and QC lecture.pptx
Basic BS-Seq processing Exercises.docx

  • Differential Methylation Analysis，如何计算差异甲基化，以及相应的统计学方法选择。

Differential Methylation lecture.pptx

  • Visualising and Exploring BS-Seq Data，讨论了关于如何计算位点和区域的甲基化水平

Visualising and Exploring Lecture.pptx

  • 偏向于介绍seqmonk的使用，当然也介绍了其他的知识点
  1. Differential methylation Exercises.docx
  2. Visualising and exploring Exercises.docx
  3. SeqMonk tools for methylation analysis.pptx

RRBS数据分析流程

RRBS甲基化分析流程：https://blog.csdn.net/qq_29300341/article/details/80551131
我对这个教程最感兴趣的地方是他使用R包来进行基因组区域的注释并计算其甲基化水平，代码如下：

  1. library(GenomicRanges)
  2. library(IRanges)
  3. library(TxDb.Mmusculus.UCSC.mm10.ensGene)
  5. txdb_mm10 <- TxDb.Mmusculus.UCSC.mm10.ensGene
  6. trans <- as.data.frame(transcripts(txdb_mm10))
  8. # 取出正链上的基因
  9. transUp <- trans[trans$strand=="+",]
  11. # 取出Tss上游200到下游1500bp内的范围
  12. grPromoterUp <- GRanges(seqnames = transUp$seqnames,strand = "+",
  13.                  ranges = IRanges(start = transUp$start-1500,end = transUp$start+200))
  15. # # 取出负链上的基因
  16. # transDown <- trans[trans$strand=="-",]
  17. # # 取出Tss上游200到下游1500bp内的范围
  18. # grPromoterDown <- GRanges(seqnames = transDown$seqnames,strand = "-",
  19. #                       ranges = IRanges(start = transDown$end-200,end = transDown$end+1500))
  20. grPromoterUp$value <- NaN #value存储每个promoter区域的甲基化程度
  21. wt01Cpg <- read.table("C:/Users/f/Desktop/wt01_CpG_filter.txt",header = TRUE)
  22. grCpG <- GRanges(seqnames = wt01Cpg$chr,
  23.                  ranges = IRanges(start = wt01Cpg$base, width = 1))
  24. grCpG$meth <- wt01Cpg$freqC
  25. grCpG$coverage <- wt01Cpg$coverage
  27. #---------最重要的就是GenomicRanges包的根据基因组坐标找交集的findOverlaps函数
  28. hitObj<- findOverlaps(grPromoterUp,grCpG)
  30. promoterid<- unique(hitObj@from)
  31. for(i in 1:length(promoterid)){
  32.   CpGid<- hitObj[hitObj@from==promoterid[i],]@to
  33.   promoterCpG <- grCpG[CpGid]
  34.   grPromoterUp[promoterid[i]]$value<- weighted.mean(promoterCpG$meth,promoterCpG$coverage) #加权平均值
  35. }
  36. grPromoterUp <- as.data.frame(grPromoterUp)
  37. grPromoterUp <- cbind(grPromoterUp,transUp$tx_id,transUp$tx_name)
  38. write.table(grPromoterUp,"wt01UpStrandPromoterMeth.txt",sep = "\t",quote = FALSE,row.names = FALSE)
  39. #Further :把负链上的基因也做一遍

表观组学 - 甲基化分析
这篇博文里面的

注意细节

是否去重

WGBS在比对完后需要去重；RRBS不需要去重；
去重指的是去除PCR引入的重复扩增的reads，但这些重复的reads其实都是来源于同一条被超声波或机械或酶学打断的序列。

bismark软件的输出文件格式

非常值得注意的一个问题是：
bismark coverage 文件格式里面记录的CpG位点数目 比 全基因组甲基化胞嘧啶报告，即后缀为.”CpG_report.txt”文件格式（就算去掉没有覆盖到的位点）数目要多一点点，这个是正常的，bismark的GitHub的常见问题里面有解答。

文库是否具有方向性/链特异性

image.png
参考：单细胞甲基化测序（scBS-seq）比对率奇低？你的文库方向性参数可能没选对！

DNA甲基化检测手段.png

关于RRBS使用**MspI**（甲基化敏感的限制性核酸内切酶）酶切C^CGG位点的生物学小知识点：
下面的文字来自于：https://www.thermofisher.cn/order/catalog/product/cn/zh/ER0541
image.png
image.png