2021-01-01 更新
簡單記錄一下平時(shí)的分析過程
和RNA-seq前期流程類似 -- 質(zhì)控、去接頭、比對參考基因組、排序
后期就是要提取甲基化位點(diǎn)，包括CpG、CHG、CHH三種context，H代表非G位點(diǎn)（A、C、T）。得到bedgraph文件后將個(gè)樣本匯總為一個(gè)GR (GenomicRanges)文件，便于后續(xù)分析
一、質(zhì)控、去接頭
1. fastqc 質(zhì)控
1.1 主要參數(shù)
- -h
打印幫助文檔 - -o
結(jié)果輸出文件夾，必須存在，因?yàn)椴粫詣?dòng)被創(chuàng)建 - --extract
將結(jié)果壓縮文件解壓，不需要解壓則不需使用該參數(shù) - -t
線程數(shù)，設(shè) 2 就足夠了 - -a
指定需要額外去除的接頭文件 - -q
靜默運(yùn)行，不報(bào)告標(biāo)準(zhǔn)輸出，只報(bào)告錯(cuò)誤信息 - fastq1 (fastq2)
最后就是fastq文件了，不需要加參數(shù)名，可以單端可以雙端
1.2 示例
fastqc -q -o out_dir fastq_R1 fastq_R2
更多信息需要你自己查看幫助文檔和 FastQC 官方手冊.pdf
另外，官方網(wǎng)頁版 的有對每一模塊進(jìn)行詳盡地解釋，并對給出警告或錯(cuò)誤的可能原因，針不戳！
-q: 靜默運(yùn)行，不打印處理過程
-o: 結(jié)果輸出文件夾
- 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
額，去接頭好像還沒寫，改天一定。
質(zhì)控 -- fastqc； 去接頭 -- trim-galore
二、比對基因組
- 安裝 bismark
# 從github下載
git clone https://github.com/FelixKrueger/Bismark.git
cd Bismark
make
# 添加環(huán)境變量
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
# 安裝成功
- 構(gòu)建參考基因組
# 進(jìn)入保存基因組的文件，先下載基因組文件
cd Genome
mkdir hg19
cd hg19
wget -c http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz &
# 構(gòu)建甲基化比對基因組文件
bismark_genome_preparation hg19
另一個(gè)軟件是一個(gè) BSMAP，目前用的是后者，兩者關(guān)系也有很多人介紹，我也不知道。
- 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: 質(zhì)量閾值，低于該值的舍棄
- -f: 去除低質(zhì)量的reads，將包含Ns > n (n指設(shè)的值)的reads舍棄
- -p: 參與處理的處理器個(gè)數(shù)
- -r: 如何報(bào)告重復(fù)的Hits，不知道啥意思
- -v: 所容許的錯(cuò)配率
- -s: seed大小，WGBS模式為16，RRBS模式為12，最小8最大16
- -S: 用于隨機(jī)數(shù)生成的種子在選擇多個(gè)命中時(shí)使用其他種子值根據(jù)讀取的索引編號生成偽隨機(jī)數(shù)，以允許再現(xiàn)映射結(jié)果。默認(rèn)值為0。（從系統(tǒng)時(shí)鐘獲取種子，映射結(jié)果不可重現(xiàn)。）翻譯的，同不知道啥意思
- -n: set mapping strand information
- 2> 標(biāo)準(zhǔn)錯(cuò)誤重定向，在這里是保存BSMAP的報(bào)告內(nèi)容
- samtools view -b : 將SAM文件輸出為BAM文件
- -o: 輸出的BAM文件名
這樣就得到了BAM文件了
三、排序、去重
- 排序
首先按照比對的基因組坐標(biāo)進(jìn)行排序
sambamba sort \
-m 8GB \
--tmpdir tmp \
-t 5 \
-o output_sorted.bam input_file.bam
- -m: 所有線程的存儲上限
- --tmpdir: 臨時(shí)文件夾
- -t: 線程數(shù)
- -o: 輸出文件的路徑名
- 去重
去除多重比對、重復(fù)、未比對上的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
輸入為上文創(chuàng)建的bismark基因組文件和排序去重后的BAM文件
- --opref: 要保存文件的前綴
- 如果只是這樣后面不再加其他參數(shù)，則默認(rèn)提取CpG位點(diǎn)，輸出文件為'sample_CpG.bedGraph'
- 如果后面再加上'--CHG'，則為提取CHG位點(diǎn)，輸出文件為'sample_CHG.bedGraph'
同理：'--CHH' 對應(yīng) 'sample_CHH.bedGraph'
至此，所有CpG位點(diǎn)就全部被提取出來了。
五、將CpG位點(diǎn)保存為GR文件
由于測序是區(qū)分正負(fù)鏈的，而在分析的時(shí)候不區(qū)分，所以需要合并正負(fù)鏈的信息。
還需要將與基因組CpG位點(diǎn)不匹配的位點(diǎn)去除，因此需要load一個(gè)全基因組CpG位點(diǎn)文件。
細(xì)節(jié)我就不寫了，只寫主要操作，即將每個(gè)樣本循環(huán)保存為GR文件放入一個(gè)list里面，最后再unlist一下，就變成了一個(gè)包含所有樣本的所有CpG位點(diǎn)信息的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"))
