以PacBio为代表的三代测序公司开发的第三代测序技术既可以直接测得DNA碱基序列，又可以测得其碱基上的修饰。
SMRT(Single-Molecule Real-Time)测序即单分子实时测序。变合成边测序的过程中，携带修饰的碱基所发出的荧光信号滞留时间有所差异，通过检测信号脉冲的长短即可区分普通碱基和携带修饰的碱基。
SMRT Sequencing由于信号信息都记录在原始文件中（subreads，BAM格式），故基于一定深度（往往很深，50x、100x等，具体修饰类型不一样，目前费用自然比较高）的subreas信息即可检测全基因组范围内的碱基修饰和motif。
那么，如何基于现有的软件进行剪辑修饰的检测呢 ？
一、subreads比对
使用官方SMRTlinks（已经更新到6.0版本）中的pbalign对subreads BAM文件进行比对
Mapping PacBio sequences to references using an algorithm selected from a selection of supported command-line alignment algorithms. Input can be a fasta, pls.h5, bas.h5 or ccs.h5 file or a fofn (file of file names). Output can be in SAM or BAM format. If output is BAM format, aligner can only be blasr and QVs will be loaded automatically.
比对命令如下：
pbalign subreads.bam hg19.fasta raw.bam --tmpDir tmp --nproc 9 --concordant
其中比对的算法支持blasr、bowtie、gmap，默认为blastr
# 若一个样本测得多个cell数据，各个cell数据分别比对后得到bam文件再进行合并：
例如：rawbam.txt记录了多个bam文件的路径
$path/cell1.raw.bam
$path/cell2.raw.bam
$path/cell3.raw.bam
接着对BAM文件进行文件上合并：
samtools merge -b rawbam.txt raw.merged.bam -@ 8
得到合并后的比对文件，即raw.merged.bam
或者，使用SubreadSet批量进行（比较推荐这种方式）：
dataset create bam.xml bam.list --type SubreadSet
pbalign $PWD/bam.xml hg19.fasta raw.bam --tmpDir $PWD/tmp --nproc 8 --concordant
二、检测甲基化位点
得到比对结果bam文件后，利用SMRTlink自带的的ipdSummary软件进行甲基化位点的检测。
# index genome file
samtools faidx hg19.fasta
# modification identify
ipdSummary raw.bam \
--reference hg19.fasta \
--gff basemods.gff \
--csv kinetics.csv \
--bigwig IpdRatio.bw \
--identify m6A,m4C \
--methylFraction \
--numWorkers 10 \
--pvalue 0.01 \
--minCoverage 3
其中一些关键性的阈值参数设置描述如下：
--identify： Specific modifications to identify (comma-separated list). Currrent options are m6A, m4C, m5C_TET. Using --control overrides this option. (default: m6A,m4C)
--minCoverage： Minimum coverage required to call a modified base (default: 3)
--methylMinCov： Do not try to estimate methylFraction unless coverage is at least this. (default: 10)
GFF结果示例（basemods.gff）：
chrM kinModCall modified_base 16569 16569 36 + . coverage=256;context=CCCTTAAATAAGACATCACGATGNNNNNNNNNNNNNNNNNN;IPDRatio=1.60
CSV结构所示例（kinetics.csv）：
refName,tpl,strand,base,score,tMean,tErr,modelPrediction,ipdRatio,coverage,frac,fracLow,fracUp
"chrM",1,0,G,27,1.364,0.190,0.773,1.765,148,,,
参考：
https://github.com/PacificBiosciences/kineticsTools/blob/master/doc/manual.rst
网友评论