专栏首页生信技能树基因组重测序的unmapped reads assembly探究 【直播】我的基因组86

基因组重测序的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),作者:曾健明

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

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

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 100篇泛癌研究文献解读之snoRNAs

    为了分析不同类型、组织起源肿瘤的共性、差异以及新课题。TCGA于2012年10月26日-27日在圣克鲁兹,加州举行的会议中发起了泛癌计划。参考:https://...

    生信技能树
  • 时隔四年,无意中看到了双重高斯分布拟合

    如果我们有比较强的背景知识,或者看了如下分布的条形图,会下意识的猜想出是两个高斯分布的混合,但是想从数据的角度来探索,两个独立的高斯分布各自独立的均值和方差该如...

    生信技能树
  • 生信蓝领,一个不舍得分享的高通量数据分析框架

    安装bcbio框架 软件安装 配置参考基因组 流程配置 手动创建 脚本创建 简单实战 总结 当我跑完一些分析流程,比如说RNA-Seq,重测序分析以后,我就想到...

    生信技能树
  • 推理应用中内存体系结构的能量延迟精度的基本限制(computer science)

    本文获得了内存计算体系结构(IMC)的计算精度的基本限制。定义了一个IMC噪声模型和相关的SNR度量,并分析了它们之间的相互关系,以表明IMC的精度从根本上受到...

    用户8054058
  • Linux-批量替换文件内容

    瑞新
  • bug 回忆录(一)

    @author Ken @time 2020-09-24 22:00:55 @description 转载请备注出处,谢谢

    公众号---人生代码
  • ReplaceText

    在下面的列表中,必需属性的名称以粗体显示。任何其他属性(不是粗体)都被认为是可选的,并且指出属性默认值(如果有默认值),以及属性是否支持表达式语言。

    用户5475193
  • 使用c#对MongoDB进行查询(1)

        在MongoDB.Bson命名空间下存在一个BsonDocument类,它是MongoDB的文档对象,代表着MongoDB中不规则数据一条条实体模型。可...

    莫问今朝
  • leetcode栈之有效的括号

    这里借助栈,然后遍历每个char,针对 (、 [、 {分别push对应配对的char,其他的则判断stack是否为空或者pop出来的值是否与之相等,如果不等则返...

    codecraft
  • leetcode栈之有效的括号

    这里借助栈,然后遍历每个char,针对 (、 [、 {分别push对应配对的char,其他的则判断stack是否为空或者pop出来的值是否与之相等,如果不等则返...

    codecraft

扫码关注云+社区

领取腾讯云代金券