基因组重测序的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 条评论
登录 后参与评论

相关文章

来自专栏章鱼的慢慢技术路

用ARM实现音乐电子相册

1502
来自专栏后端云

compute node ha 主流开源实现

nova evacuate和热迁移很像。都是想实例从一个节点转移到另外一个节点。区别主要是热迁移在正常状态下进行的,疏散时在异常状态下进行的。用一个形象的比如就...

1093
来自专栏IT技术精选文摘

LVS集群的负载调度

1. 前言 本文将主要讲述在负载调度器上的负载调度策略和算法,如何将请求流调度到各台服务器,使得各台服务器尽可能地保持负载均衡。文章主要由两个部分组 成。第...

2049
来自专栏吉浦迅科技

DAY91:阅读Programming Model

我们正带领大家开始阅读英文的《CUDA C Programming Guide》,今天是第91天,我们正在讲解Unified Memory Programmin...

571
来自专栏字根中文校对软件

JcJc人工智能错别字校对系统API(网站错别字检查)

JcJc人工智能错别字校对系统API - 网站错别字检查 中文校对软件 http://cuobiezi.net/zigen_jiaodui.html App...

3955
来自专栏玉树芝兰

Tensorflow执行pip升级安装的坑

(由于微信公众号外部链接的限制,文中的部分链接可能无法正确打开。如有需要,请点击文末的“阅读原文”按钮,访问可以正常显示外链的版本。)

503
来自专栏ATYUN订阅号

深度学习图像识别项目(上):如何快速构建图像数据集

AiTechYun 编辑:yuxiangyu 在你还是个孩子时是否也是一个神奇宝贝迷?是否还记得里面的各种神奇宝贝,以及小智手中可以自动识别神奇宝贝的图鉴(Po...

1.1K6
来自专栏吉浦迅科技

在NVIDIA Jetson TX2上运行 openpose

CMU 开源的一个实时多人关键点检测库( real-time multi-person keypoint detection),包括人体关键点,手部关键点,脸部...

1852
来自专栏李蔚蓬的专栏

关于TensorFlow的安装和心得

随后在cmd处通过import tensorflow,查看本机的CUDA以及cuDNN的适机版本:

1193
来自专栏开源FPGA

基于FPGA的Uart接收图像数据至VGA显示

系统框图 ?   前面我们设计了基于FPGA的静态图片显示,并对一幅彩色图片提取了灰度,学习了RGB转Gray算法。这是基于一幅静态图片的,那么后面我们要怎么模...

2349

扫码关注云+社区