2021-10-27 [WGS] pacbio third generation methylation modification process

Abstract

I was very busy some time ago , One is that there are many projects , The other is a personal matter , Near the end of the month, I finally have time to continue coding articles . This article introduces the basic process analysis of the third generation methylation . After analyzing the methylation modification of the sequence during sequencing , Use SMRT Official tools for analysis , obtain m4C,m6A,m5C_TET Notes .

Methods and tools

Sequencing instrument ：Pacbio
Analysis tools ：
assemble ：Canu;flye
comparison ：pbmm2;samtools(SMRTlink Bring their own )
notes ：ipdSummary,motifmaker(SMRTlink Bring their own )
SMRTlink Installation document
SMRTlink Reference documents

Operation process

assemble

for bam_path in */1.rawdata/*.bam; # Read bam File path and sample name 
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 It can also be assembled by itself , Here we only use its correction function , Use flye Carry out assembly （ Former colleagues pipeline Reserved ）. Of course, I didn't go into the comparison between the two tools , If you have a better process, you can leave a message to discuss .

comparison

for bam_path in */1.rawdata/*.bam; # Read bam File path and sample name 
do
sample_id=${bam_path%/1.rawdata/*};
echo ${bam_path} ${sample_id};
mkdir ../02.Assemble/ref ../05.align# Create reference and comparison folders 
pbmm2 index ../02.Assemble/flye/${sample_id}_assembly.fasta ../02.Assemble/ref/${sample_id}_assembly.mmi # Index the reference sequence 
pbmm2 align ../02.Assemble/flye/${sample_id}_assembly.fasta ${bam_path} ../05.align/${sample_id}_tmp.bam # comparison 
samtools sort ../05.align/${sample_id}_tmp.bam  ../05.align/${sample_id} # Arrange files 
done

notes

for bam_path in */1.rawdata/*.bam; # Read bam File path and sample name 
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

Result display

You can get 4 results ,gff and csv various 2 File （M. The prefix is added by yourself , Used to distinguish samples ）

basemods

basemods.gff The document provides m4C,m6A,m5C_TET Note information for ,csv Provides information about the mutation site

motif

motif.gff Comment files and basemods.gff The front is almost the same , There are only some differences in the last column of annotation information of some sites , More than a id=BSC;motif.csv Provide information on the position of the methylated fragment sequence , If there are too few sequences , Can be in motifmaker You can set -m Parameter adjustment filter score .

summary

The third generation methylation modification is not complicated , The tools are easy to install and use . I feel that the follow-up is to add some analysis , Such as making statistics , Continue to draw some pictures . Between the current customer needs , Not added yet , Interested friends can join a group to communicate . When the QR code expires, you can add VX：bbplayer2021 , remarks Apply to join the student information exchange group .