2021-10-27【WGS】|Pacbio三代甲基化修饰流程

目录

摘要

前段时间特别忙,一个是项目多,另一个是个人私事,临近月底终于有空可以继续码文章。本篇介绍的是三代甲基化的基本流程分析。在测序时分析序列的甲基化修饰后,使用SMRT官方工具进行分析,得到m4C,m6A,m5C_TET的注释。
2021-10-27【WGS】|Pacbio三代甲基化修饰流程

方法与工具

测序仪器:Pacbio
分析工具:
组装:Canu;flye
比对:pbmm2;samtools(SMRTlink自带)
注释:ipdSummary,motifmaker(SMRTlink自带)
SMRTlink安装文档
SMRTlink参考文档

操作流程

组装

bam2fastq转换格式 canu矫正/flye组装 原始数据.bam 原始数据.fastq 组装序列.fasta
for bam_path in */1.rawdata/*.bam; #读取bam文件路径和样品名
do
sample_id=${bam_path%/1.rawdata/*};
bam2fastq ${sample_id}.bam -o ${sample_id}
canu -correct \
     -p ${sample_id} -d ../02.Assemble/${sample_id} \
    genomeSize=6.5m \
    minReadLength=2000 minOverlapLength=500\
    corOutCoverage=300 corMinCoverage=2 \
    -pacbio ${sample_id}.fastq.gz
    flye --plasmids --pacbio-corr ${sample_id}/${sample_id}.correctedReads.fasta.gz -g 6.5m -o flye/${sample_id}_flye -t 16
done

Canu 自身也可以组装,在这里我们只用了他的校正功能,使用flye进行组装(前同事的pipeline被保留)。当然这两个工具的对比我没有去深究,如果大家有更好的流程可以留言讨论。

比对

pbmm2建立索引.mmi -- pbmm2比对 samtools sort排列 组装序列.fasta 数据比对 原始数据.bam 临时文件.tmp.bam 比对文件.bam
for bam_path in */1.rawdata/*.bam; #读取bam文件路径和样品名
do
sample_id=${bam_path%/1.rawdata/*};
echo ${bam_path} ${sample_id};
mkdir ../02.Assemble/ref ../05.align#建立参考和比对文件夹
pbmm2 index ../02.Assemble/flye/${sample_id}_assembly.fasta ../02.Assemble/ref/${sample_id}_assembly.mmi #建立参考序列索引
pbmm2 align ../02.Assemble/flye/${sample_id}_assembly.fasta ${bam_path} ../05.align/${sample_id}_tmp.bam #比对
samtools sort ../05.align/${sample_id}_tmp.bam  ../05.align/${sample_id} #排列文件
done

注释

samtools建立索引.fai samtools建立索引.bai 组装序列.fasta motif.gff 比对文件.bam motifs.csv
for bam_path in */1.rawdata/*.bam; #读取bam文件路径和样品名
do
mkdir -p ../05.align/methylome/${sample_id}
samtools index ../05.align/${sample_id}.bam
samtools faidx ../02.Assemble/flye/${sample_id}_assembly.fasta
ipdSummary ../05.align/${sample_id}.bam --reference ../02.Assemble/flye/${sample_id}_assembly.fasta --gff ../05.align/methylome/${sample_id}/basemods.gff --csv ../05.align/methylome/${sample_id}/basemods.csv --pvalue 0.001 --numWorkers 16 --identify m4C,m6A,m5C_TET 
motifMaker find -f ../02.Assemble/flye/${sample_id}_assembly.fasta -g ../05.align/methylome/${sample_id}/basemods.gff -o ../05.align/methylome/${sample_id}/motifs.csv
motifMaker reprocess -f ../02.Assemble/flye/${sample_id}_assembly.fasta -g ../05.align/methylome/${sample_id}/basemods.gff -m ../05.align/methylome/${sample_id}/motifs.csv -o ../05.align/methylome/${sample_id}/motifs.gff

结果展示

可以得到4个结果,gff和csv各2个文件(M.前缀是自己加的,用于区分样品)
2021-10-27【WGS】|Pacbio三代甲基化修饰流程

basemods

basemods.gff文件提供了m4C,m6A,m5C_TET的注释信息,csv提供了突变位点的信息
2021-10-27【WGS】|Pacbio三代甲基化修饰流程
2021-10-27【WGS】|Pacbio三代甲基化修饰流程

motif

motif.gff注释文件和basemods.gff前面都差不多,只有部分位点最后一列注释信息上有些区别,多一个id=BSC;motif.csv提供甲基化片段序列位置信息,如果序列太少,可以在motifmaker可以设置-m 参数调整筛选分数。
2021-10-27【WGS】|Pacbio三代甲基化修饰流程
2021-10-27【WGS】|Pacbio三代甲基化修饰流程

总结

三代的甲基化修饰其实并不复杂,工具安装也很方便好用。个人感觉后续是加一些分析,比如做统计,继续绘制一些图片的。介于目前客户需求,暂时没有添加,有兴趣的朋友可以加群一起交流。遇见二维码过期可添加VX:bbplayer2021 ,备注 申请加入生信交流群。
2021-10-27【WGS】|Pacbio三代甲基化修饰流程

上一篇:Java 创建 PDF 文件包的两种方法


下一篇:「ABC221」B - typo 题解