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

相关文章

来自专栏测试开发架构之路

Filter Blue Light for Better Sleep(APP 推荐)

Filter Blue Light for Better Sleep By Carolyn Mohr 11 May, 2016 Many people lik...

33316
来自专栏知晓程序

用这个小程序,做独一无二的微信表情包!

本期,知晓程序(微信号 zxcx0101)推荐的「表情家园」小程序,让你从今以后能更愉快地聊天,更愉快地斗图。

402
来自专栏专知

【论文Slides分享】多媒体顶级会议ACM Multimedia 2017 China 论文宣讲会报告Slides资料分享

【导读】第25届ACM国际多媒体会议(ACM International Conference on Multimedia, 简称ACMMM)于2017年10月...

3016
来自专栏机器学习人工学weekly

机器学习人工学weekly-2018/4/15

注意下面很多链接需要科学上网,无奈国情如此 1. DeepMind的新工作,不用地图在城市里导航 Learning to navigate in cities ...

3308
来自专栏PaddlePaddle

【AI核心技术】课程十九:神经图灵机—寻址

UAI与PaddlePaddle联合推出的【AI核心技术掌握】系列课程持续更新中!

671
来自专栏企鹅号快讯

如何使用EndNote Match

貌似大家差不多都逃不出脱单、六级、雅思、SCI、考研、考博...的魔咒 那小博只好准备一篇 如何使用EndNote Match教程给你们祝早日出刊(笔芯~) 第...

1749
来自专栏Java技术栈

程序员最大的绝望(不是bug)

对于程序员来说,最大的烦恼可能并不是电脑里的bug,而更多的是坐在电脑前的人。 你知道PICNIC是什么意思吗?它并不是你想的那样是一次惬意的野餐。这是一个缩写...

3326
来自专栏量子位

AI续写的《冰与火之歌》怎么样?国外小哥录了一个评书版

允中 发自 凹非寺 量子位 报道 | 公众号 QbitAI ? 读……读粗来了…… 你知道《凛冬的寒风》么? 这个名字,归属于乔治啊啊马丁《冰与火之歌》系列小说...

2503
来自专栏Youngxj

超炫经典HTML5游戏 附游戏源码

2178
来自专栏CVer

糟了!TensorFlow 1.9.0 来了

关注 CVer 的同学应该知道,前不久Amusi整理了 TensorFlow相关的学习资料,并推出TensorFlow从入门到精通系列贴(最近事情比较多,更新的...

1053

扫描关注云+社区