基因组重测序的unmapped reads assembly探究 【直播】我的基因组86

在前面的直播基因组系列,我们讲解过那些比对不少我们人类的参考基因组序列的数据,其实可以细致的进行探究。

直播】我的基因组(十五):提取未比对的测序数据

这里主要参考这篇文章的图4:http://www.nature.com/ng/journal/v42/n11/figtab/ng.691F4.html

组装的contig注释到物种

这是2010年发表于nature genetics杂志的Whole-genome sequencing and comprehensive variant analysis of a Japanese individual using massively parallel sequencing 虽然文章选择的是SOAPdenovo,ABySS,Velvet这3款软件来进行组装,但毕竟是2010年的文章了,现在其实有更好的选择,比如Minia

选择Minia工具来组装

Minia软件也是基于de Bruijn图原理的短序列组装工具,优于以前的ABySS和SOAPdenovo,关键是速度非常快,十几分钟就OK了,不消耗计算机资源,所以这里就选择它啦。

下载安装Minia

安装官网的指导说明书下载二进制版本即可,代码如下:

## Download and install Minia# http://minia.genouest.org/cd ~/biosoftmkdir Minia &&  cd Miniawget https://github.com/GATB/minia/releases/download/v2.0.7/minia-v2.0.7-bin-Linux.tar.gz tar -zxvf minia-v2.0.7-bin-Linux.tar.gz ~/biosoft/Minia/minia-v2.0.7-bin-Linux/bin/minia --help ## eg: ./minia -in reads.fa -kmer-size 31 -abundance-min 3 -out output_prefix 

软件使用方法也非常简单,就一行命令,其中最佳 -kmer-size需要用KmerGenie来确定。

使用

step1:提取比对失败的reads

samtools view -f4 jmzeng_recal.bam |perl -alne '{print "\@$F[0]\n$F[9]\n+\n$F[10]" }' >unmapped.fqperl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fastq unmapped.fq -graph_data unmapped.gd -out_good null -out_bad nullperl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i unmapped.gd -png_all -o unmappedperl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i unmapped.gd -html_all -o unmappedcd ~/data/project/myGenome/gatk/jmzeng/unmapped

共31481084/4=7870271,仅仅是7.8M的reads

Input Information

Input file(s):

unmapped.fq

Input format(s):

FASTQ

# Sequences:

7,870,271

Total bases:

1,180,540,650

step2: 用KmerGenie确定kmer值

KmerGenie estimates the best k-mer length for genome de novo assembly.

KmerGenie predictions can be applied to single-k genome assemblers (e.g. Velvet, SOAPdenovo 2, ABySS, Minia).

## http://kmergenie.bx.psu.edu/cd ~/biosoftmkdir KmerGenie &&  cd KmerGeniewget http://kmergenie.bx.psu.edu/kmergenie-1.7044.tar.gztar zxvf kmergenie-1.7044.tar.gzcd kmergenie-1.7044make python setup.py install --user~/.local/bin/kmergenie --help cd ~/data/project/myGenome/gatk/jmzeng/unmapped~/.local/bin/kmergenie unmapped.fq

step3: 运行Minia

cd ~/data/project/myGenome/gatk/jmzeng/unmapped~/biosoft/Minia/minia-v2.0.7-bin-Linux/bin/minia  -in unmapped.fq -kmer-size 31 -abundance-min 3 -out output_prefix

7.8M的reads组装之后有272007条contigs

组装之后:

Prinseq v0.20.4 was used to calculate assembly statistics, including N50 contig size, GC content

cd ~/data/project/myGenome/gatk/jmzeng/unmappedperl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fasta output_prefix.contigs.fa  -graph_data contigs.gd -out_good null -out_bad null perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i contigs.gd -png_all -o contigsperl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i contigs.gd -html_all -o contigsperl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fasta output_prefix.contigs.fa  -stats_assembly

就是给出一些指标,如下;

stats_assembly    N50 176stats_assembly    N75 113stats_assembly    N90 78stats_assembly    N95 70

Input Information

Input file(s):

output_prefix.contigs.fa

Input format(s):

FASTA

# Sequences:

272,007

Total bases:

44,868,011

Length Distribution

Mean sequence length:

164.95 ± 204.44 bp

Minimum length:

63 bp

Maximum length:

10,187 bp

Length range:

10,125 bp

Mode length:

150 bp with 16,461 sequences

然后用RNA-SEQ数据来比对验证! 以后再讲

把组装好的contigs拿去NCBI做blast看看物种分布,Distribution of top nucleotide BLAST hits by species from the NCBI nr database for 1000 random contigs in the assembly!其实上面的prinseq软件也简单的给出了一个污染物种分布情况表,但是这个原理不一样。以后再讲

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-09-08

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏CSDN技术头条

JMeter 怎么学?

大家在网上搜寻许多关于 JMeter 的应用案例时是不是有过这样的遭遇: 明明是按照文档中的内容去做的,但是好些时候总是出错,这个时候会疯狂搜索各类与问题相关的...

4406
来自专栏安恒信息

2013.9.22--9.28 病毒预报

国家计算机病毒应急处理中心通过对互联网的监测发现,近期出现一种恶意木马程序变种Trojan_ Wapomi.A。该变种会使得受感染操作系统中的隐私信息数据、网银...

3134
来自专栏皮振伟的专栏

[virt][concept]虚拟化技术概论--KVM,QEMU,Libvirt

前言: 以作者的经验来看,虚拟化的跨度比较大,很多概念比较难以理解,本来以为“硬件行为,就是这样的”好多概念,都变成虚拟的了。 作者对kernel略懂一二,结合...

2836
来自专栏Kirito的技术分享

选择Kong作为你的API网关

在微服务架构之下,服务被拆的非常零散,降低了耦合度的同时也给服务的统一管理增加了难度。如上图左所示,在旧的服务治理体系之下,鉴权,限流,日志,监控等通用功能需要...

1613
来自专栏数据和云

Oracle 12.2中那些温暖人心的特性

在OOW 2015大会上,Oracle已经发布了12.2的Beta版本,其中的很多亮点新特性引人瞩目,包括在IMO和Multitenant方面,以及在Shard...

3606
来自专栏杨建荣的学习笔记

关于抓取session信息的一个脚本(r3笔记第8天)

关于session的诊断,可以基于动态性能视图,ash,awr.. 自己也写过一些简单的脚本,在平时的工作中也能够完成一些基本的工作。今天在看taner分享的脚...

3196
来自专栏一棹烟波

CUDA共享内存的使用示例

CUDA共享内存使用示例如下:参考教材《GPU高性能编程CUDA实战》。P54-P65 教材下载地址:http://download.csdn.net/down...

1828
来自专栏ChaMd5安全团队

RFIDHacKing频射硬件入门

RFIDHacKing频射硬件入门 From ChaMd5安全团队核心成员 MAX丶 鉴于硬件安全对于大多数新人是较少接触的,而这方面又非常吸...

3559
来自专栏杨建荣的学习笔记

无奈的硬件故障 (r7笔记第20天)

最近真是忙的厉害,感觉时间都不是自己的了,大周末的时间都排得满满当当,先是大半夜接到报警电话,接着碰到了让人无奈的硬件问题,一台服务器挂掉,结果上面有两个备库,...

2593
来自专栏腾讯云安全的专栏

你熟悉的Android Root 方式有哪些?|附演示视频

2055

扫码关注云+社区