跟着Briefings in Bioinformatics学数据分析：植物线粒体基因组组装流程GSAT初步尝试

论文

Master graph: an essential integrated assembly model for the plant mitogenome based on a graph-based framework

https://academic.oup.com/bib/article-abstract/24/1/bbac522/6854450?redirectedFrom=fulltext&login=false

bbac522.pdf

软件的github链接 https://github.com/hwc2021/GSAT

论文作者在CGM基因组沙龙做了学术报告，录屏在B站可以找到，链接是 https://www.bilibili.com/video/BV1v24y1h7Nw/?spm_id_from=333.999.0.0

论文中作者利用自己开发的流程组装了拟南芥和水稻的图形化的线粒体基因组，我找到了论文中拟南芥的数据，重复了一下整个过程，这篇推文做一个记录

拟南芥的全基因组测序数据，包括ngs hifi ont 来源于论文 https://www.sciencedirect.com/science/article/pii/S1672022921001741

数据可以从国家基因组科学数据中心下载，之前也写了推文进行介绍

这个流程的第一步是使用spades这个软件利用二代测序数据进行一个初步的组装，这一步如果是利用完整的全基因组测序数据去做的话时间会非常长，而且得到的结果也不理想（不知道是不是我自己的做法存在问题）。

这一步我去下载了NCBI已经有的拟南芥线粒体基因组序列，把二代全基因组测序数据比对到这个线粒体基因组，然后提取比对上的reads，最后用比对上的reads去做GSAT流程的第一步

比对（bowtie2）提取reads的代码

bowtie2-build at_mt.fasta at_mt.index/atmito
bowtie2 -p 8 --local -x at_mt.index/atmito -1 ../ngs/CRR302670_f1.fastq.gz -2 ../ngs/CRR302670_r2.fastq.gz -S atMito.sam
samtools sort -@ 8 atMito.sam -O bam -o atMito.sorted.bam
samtools index atMito.sorted.bam
grep '>' at_mt.fasta | awk '{print $1}' | awk 'gsub(">","")' > chr.list
python get_mapped_paired_end_reads_from_bam.py --bam atMito.sorted.bam --chrlist chr.list --or1 r1.fq --or2 r2.fq
bgzip r1.fq
bgzip r2.fq

然后按照GSAT 这个流程的github主页的文档去运行，运行这个流程需要提供config文件，config文件里的参数除了文件路径其他的我都没有改，因为还不明白每个参数起到的作用

这个流程是用perl语言写的，第一次使用可能会提示有一些perl的模块没有安装，需要对应的安装好相应的模块，因为我不太熟悉perl语言，perl模块的安装对于我来说一直是一个难点，还好这次没有遇到太多的报错

流程第一步

~/biotools/GSAT/bin/gsat graphShort -conf example.conf.short.reads

得到结果

image.png

这一步用到的是二代测序数据，接下来的内容全部用到的是三代测序数据

作者在论文里写道是推荐用Hifi数据的，但是我hifi数据还没有下载，暂时用ont的数据做下面的三个步骤

首先是选取一部分数据，因为总的数据量很大

zcat ../ont/CRR302667.fastq.gz | ~/biotools/seqkit/seqkit sample -p 0.2 -s 1234 -o at.ont.long.reads.fq.gz

~/biotools/seqkit/seqkit fq2fa at.ont.long.reads.fq.gz -o at.ont.long.reads.fa

第二步流程

~/biotools/GSAT/bin/gsat graphLong -conf example.conf.long.reads.01

得到结果

image.png

第三步流程

~/biotools/GSAT/bin/gsat graphSimplification -conf example.conf.long.reads.02

得到结果

image.png

第四步流程

~/biotools/GSAT/bin/gsat graphCorrection -conf example.conf.long.reads.03

得到结果

image.png

大体的流程能够走下来，但是后续如何处理还不明白，论文中提到的分析也还有好多看不懂，还需要多看几遍

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！

image.png