新冠病毒数据分析(8):你也可以拼一条病毒基因组

前面我们介绍过,病毒的基因组虽然很小,但确是最难拼接的,这看起来似乎不合逻辑,按理说只有几十K应该很容易拼接才是,但如果你真正尝试去拼接就知道了。因为病毒基因组突变率高,并且测序覆盖度高,一般都1000X以上,在加上测序错误,这几个原因综合起来,就让两个片段连接出现太多可能性,导致无法连接,也就连接不起来了。举个简单的例子,我不怕妖魔鬼怪,但是却被一直小小的蚊子搞得一夜难眠。


病毒基因组拼接方法

由于不能直接使用原始测序数据进行拼接,目前病毒基因组主要采用基于参考序列指导(reference guide)的方法,所谓生成一致性序列(consuses)的方法。首先将测序数据与最近源的参考序列进行比对,得到每个位点的覆盖情况,然后利用软件,挑选每个位点最佳的碱基。比如参考序列为A,比对到这个位点上有1000个碱基,其中900个为T,100个为A,那么我们就选择T作为待拼接基因组的碱基类型。最终我们会得到一条与参考基因组长度一致,只是部分位点有所差异的基因组,来作为最终得到的病毒基因组序列。例如下图中,参考序列的基因组为GTCTG,一致性序列的结果为GACTC。

这种生成一致性序列的方法最开始是用于变异检测分析中的方法,后来用于病毒基因组拼接。这种方法显然有很大的问题。这种方法有个前提,就是参考序列一定要非常近源,比如参考序列是29903个碱基,新测序的样品是否有可能是30000个碱基呢,但是利用这种方法得到的基因组序列必然还是29903个碱基,剩下的97个碱基信息就没有了。所以你看到目前拼接出来的基因组都是29903,虽然后面通过碱基Polish可以插入或者删除部分序列,但整体长度变化不大。如果第一个参考序列拼接质量有问题,后面基于这个参考序列得到的新数据,都存在问题。

但是,但是,因为目前只能使用这种方法,没有更好的方法,只能这样。

你也可以拼接基因组

案例流程:本案例使用的流程来自于Erin Young/Kelly Oakeson基于illumina测序的Artic Network工作流程
案例地址

https://github.com/CDCgov/SARS-CoV-2_Sequencing/tree/master/protocols/BFX-UT_ARTIC_Illumina

安装软件

conda install -c bioconda entrez-direct sra-tools bwa samtools ivar quast

1、创建工作目录

mkdir ncov
cd ncov

2、下载数据

#非常关心美国情况,这里下载一个美国的测序数据efetch -db=nuccore -format=fasta -id=MN908947.3 >MN908947.3.fa
prefetch SRR11247077 -O ./
fasterq-dump SRR11247077.sra 

3、bwa比对

bwa index MN908947.3.fa
bwa mem -t 12 MN908947.3.fa SRR11247077.sra.fastq  | samtools sort | samtools view -F 4 -o usa.sorted.bam

4、切除引物

#下载引物列表
wget https://raw.githubusercontent.com/artic-network/artic-ncov2019/master/primer_schemes/nCoV-2019/V1/nCoV-2019.bed
#转换适合ivar格式文件
perl -ne 'my @x=split m/t/; print join("t",@x[0..3], 60, $x[3]=~m/LEFT/?"+":"-"),"n";' nCoV-2019.bed >ARTIC-V1.bed 
ivar trim -i usa.sorted.bam -b ARTIC-V1.bed -p usa.primertrim
samtools sort usa.primertrim.bam -o usa.primertrim.sorted.bam

5、生成一致性序列

samtools mpileup -A -d 6000000 -B -Q 0 --reference  MN908947.3.fa usa.primertrim.sorted.bam | ivar consensus -p usasus -n N

6、quast评估

#需要下载MN908947.3的基因列表文件gff3格式
quast.py usa.consensus.fa -r MN908947.3.fa --featuresMN908947.3.gff3 --ref-bam usa.sorted.bam --output-dir quast

7、计算覆盖深度

samtools bedcov nCoV-2019.bed.txt usa.sorted.bam
生物信息学

新冠病毒数据分析(6):下载新冠病毒原始测序数据

2020-3-30 0:12:00

生物信息学

新冠病毒数据分析(7):ARTIC Network工作流

2020-4-1 23:33:40

声明 本网站部分文章源于互联网,出于传递更多信息和学习之目的转载,并不保证内容正确或赞同其观点。
如转载稿涉及失效、版权等问题,请立即联系管理员;我们会予以修改、删除相关文章,请留言反馈
Notice: When your legal rights are being violated, please send an email to: [email protected]
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索