基因组重测序的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软件也简单的给出了一个污染物种分布情况表,但是这个原理不一样。以后再讲
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 数据分析 常见技巧和经验总结
- Go by Example 中文版: Base64 编码
- Django3.0+supervisor+uvicorn+nginx进行线上部署
- 前端杂货铺上新
- 短视频系统源代码,实现前置摄像头水平翻转
- linux配置SOCK5代理
- 前端踩坑系列《五》
- linux上安装mitmproxy
- Jmeter(二十三) - 从入门到精通 - JMeter函数 - 上篇(详解教程)
- Elasticsearch学习笔记 -- 1
- Jmeter(二十四) - 从入门到精通 - JMeter函数 - 中篇(详解教程)
- [javascript] elementui和vue下复制粘贴上传图片
- SQL Server通过创建临时表遍历更新数据
- 对于 JavaScript 中循环之间的技术差异概述
- 初识 webpack 原理——自定义插件