DNA甲基化数据分析全流程

2021-01-01 更新

简单记录一下平时的分析过程

和RNA-seq前期流程类似 -- 质控、去接头、比对参考基因组、排序

后期就是要提取甲基化位点，包括CpG、CHG、CHH三种context，H代表非G位点（A、C、T）。得到bedgraph文件后将个样本汇总为一个GR (GenomicRanges)文件，便于后续分析

一、质控、去接头

1. fastqc 质控

1.1 主要参数

-h
打印帮助文档
-o
结果输出文件夹，必须存在，因为不会自动被创建
--extract
将结果压缩文件解压，不需要解压则不需使用该参数
-t
线程数，设 2 就足够了
-a
指定需要额外去除的接头文件
-q
静默运行，不报告标准输出，只报告错误信息
fastq1 (fastq2)
最后就是fastq文件了，不需要加参数名，可以单端可以双端

1.2 示例


fastqc -q -o out_dir fastq_R1 fastq_R2

更多信息需要你自己查看帮助文档和 FastQC 官方手册.pdf

另外，官方网页版的有对每一模块进行详尽地解释，并对给出警告或错误的可能原因，针不戳！

-q: 静默运行，不打印处理过程
-o: 结果输出文件夹

trim-galore 去接头

trim_galore \
--paired \
-o path_to_out_folder/trimmed \
--rrbs \
--basename demo \ 
fastq_R1 \
fastq_R2

trim_galore --clip_R1 5 --three_prime_clip_R1 2 --rrbs -o trimmed --basename SRX1635022 .fastq.gz

RNA-seq 数据分析完整流程

额，去接头好像还没写，改天一定。

质控 -- fastqc；去接头 -- trim-galore

二、比对基因组

安装 bismark

# 从github下载
git clone https://github.com/FelixKrueger/Bismark.git
cd Bismark
make
# 添加环境变量
echo export PATH = .../software/Bismark:$PATH >> ~/.bashrc
# bismark --version


          Bismark - Bisulfite Mapper and Methylation Caller.

                       Bismark Version: v0.22.1
        Copyright 2010-19 Felix Krueger, Babraham Bioinformatics
              www.bioinformatics.babraham.ac.uk/projects/
                https://github.com/FelixKrueger/Bismark


# 安装成功

构建参考基因组

# 进入保存基因组的文件，先下载基因组文件
cd Genome
mkdir hg19
cd hg19
wget -c http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz &
# 构建甲基化比对基因组文件
bismark_genome_preparation hg19

另一个软件是一个 BSMAP，目前用的是后者，两者关系也有很多人介绍，我也不知道。

BSMAP的用法

bsmap \
-a fastq_R1 \
-b fastq_R2 
-d path_to_genome_file/hg19.fa \
-q 20 \
-f 5 \
-p 5 \
-r 0 \
-v 0.05 \
-s 16 \
-S 1 \
-n 0 \
2> path_to_report_folder/BSMAP_report.txt | \
samtools view -b \
-o path_to_aligned_folder/aligned.bam

-a: read1

-b: read2

-d: 基因组文件

-q: 质量阈值，低于该值的舍弃

-f: 去除低质量的reads，将包含Ns > n (n指设的值)的reads舍弃

-p: 参与处理的处理器个数

-r: 如何报告重复的Hits，不知道啥意思

-v: 所容许的错配率

-s: seed大小，WGBS模式为16，RRBS模式为12，最小8最大16

-S: 用于随机数生成的种子在选择多个命中时使用其他种子值根据读取的索引编号生成伪随机数，以允许再现映射结果。默认值为0。（从系统时钟获取种子，映射结果不可重现。）翻译的，同不知道啥意思

-n: set mapping strand information

2> 标准错误重定向，在这里是保存BSMAP的报告内容

samtools view -b : 将SAM文件输出为BAM文件

-o: 输出的BAM文件名

这样就得到了BAM文件了

三、排序、去重

排序

首先按照比对的基因组坐标进行排序

sambamba sort \
-m 8GB \
--tmpdir tmp \
-t 5 \
-o output_sorted.bam input_file.bam

-m: 所有线程的存储上限

--tmpdir: 临时文件夹

-t: 线程数

-o: 输出文件的路径名

去重

去除多重比对、重复、未比对上的reads

sambamba markdup \
--overflow-list-size 1000000 \
--tmpdir tmp \
-t 5 \
input_sorted.bam \
output.bam \
2> MarkDup_report.txt

--overflow-list-size: size of the overflow list where reads, thrown from the hash table, get a second chance to meet their pairs (default is 200000 reads); increasing the size reduces the number of temporary files created
还是那句话，看不懂

其他的同上

最后就得到了排序且去重的BAM文件了

四、提取甲基化信息

MethylDackel extract \
path_to_genome_file/Bismark/hg19/hg19.fa \
demo.bam \
--opref path_to_save_folder/sample

输入为上文创建的bismark基因组文件和排序去重后的BAM文件

--opref: 要保存文件的前缀

如果只是这样后面不再加其他参数，则默认提取CpG位点，输出文件为'sample_CpG.bedGraph'

如果后面再加上'--CHG'，则为提取CHG位点，输出文件为'sample_CHG.bedGraph'
同理：'--CHH' 对应 'sample_CHH.bedGraph'

至此，所有CpG位点就全部被提取出来了。

五、将CpG位点保存为GR文件

由于测序是区分正负链的，而在分析的时候不区分，所以需要合并正负链的信息。
还需要将与基因组CpG位点不匹配的位点去除，因此需要load一个全基因组CpG位点文件。

细节我就不写了，只写主要操作，即将每个样本循环保存为GR文件放入一个list里面，最后再unlist一下，就变成了一个包含所有样本的所有CpG位点信息的GR对象了

for(i in 1:Ns){

    tag = vFolder[i]
    fileIn = paste0(tag, "_CpG.bedGraph")
    cat("Processing file", fileIn, "\n")

    Tx = read.table(file = fileIn, row.names = NULL, header = FALSE, sep = "\t", skip = 1)
    chr = as.character(Tx[,1])
    ST = as.numeric(Tx[,2]) + 1
    ED = as.numeric(Tx[,3])
    methyCount = as.numeric(Tx[,5])
    readCount = methyCount + as.numeric(Tx[,6])
    gr = GRanges(seqnames = Rle(chr),  ranges = IRanges(ST,ED),  methyCount = methyCount, readCount = readCount)

    gr = gr[countOverlaps(gr, cgGR, type = "any", ignore.strand = TRUE) == 1]
    x = findOverlaps(gr, cgGR, type = "any", ignore.strand = TRUE)
    cM = data.frame(mcols(gr))
    agT = aggregate(cM, by = list(subjectHits(x)),FUN="sum")
    mgr = cgGR[as.integer(agT[,1])]
    mcols(mgr) = agT[,2:3]
    mgr$SID = tag
    grList[[i]] = mgr
}

CpG_Batch = unlist(GRangesList(grList))
save(CpG_Batch, file = paste0("RData/", args$tag, "_CpG_Batch.RData"))