查看原文
其他

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

2017-09-08 曾健明 生信技能树

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

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

这里主要参考这篇文章的图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

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

  1. ## Download and install Minia

  2. # http://minia.genouest.org/

  3. cd ~/biosoft

  4. mkdir Minia &&  cd Minia

  5. wget https://github.com/GATB/minia/releases/download/v2.0.7/minia-v2.0.7-bin-Linux.tar.gz

  6. tar -zxvf minia-v2.0.7-bin-Linux.tar.gz

  7. ~/biosoft/Minia/minia-v2.0.7-bin-Linux/bin/minia --help

  8. ## eg: ./minia -in reads.fa -kmer-size 31 -abundance-min 3 -out output_prefix

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

使用

step1:提取比对失败的reads

  1. samtools view -f4 jmzeng_recal.bam |perl -alne '{print "\@$F[0]\n$F[9]\n+\n$F[10]" }' >unmapped.fq

  2. perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fastq unmapped.fq -graph_data unmapped.gd -out_good null -out_bad null

  3. perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i unmapped.gd -png_all -o unmapped

  4. perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i unmapped.gd -html_all -o unmapped

  5. cd ~/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).

  1. ## http://kmergenie.bx.psu.edu/

  2. cd ~/biosoft

  3. mkdir KmerGenie &&  cd KmerGenie

  4. wget http://kmergenie.bx.psu.edu/kmergenie-1.7044.tar.gz

  5. tar zxvf kmergenie-1.7044.tar.gz

  6. cd kmergenie-1.7044

  7. make

  8. python setup.py install --user

  9. ~/.local/bin/kmergenie --help

  10. cd ~/data/project/myGenome/gatk/jmzeng/unmapped

  11. ~/.local/bin/kmergenie unmapped.fq

step3: 运行Minia

  1. cd ~/data/project/myGenome/gatk/jmzeng/unmapped

  2. ~/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

  1. cd ~/data/project/myGenome/gatk/jmzeng/unmapped

  2. perl ~/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

  3. perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i contigs.gd -png_all -o contigs

  4. perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i contigs.gd -html_all -o contigs

  5. perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fasta output_prefix.contigs.fa  -stats_assembly

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

  1. stats_assembly    N50 176

  2. stats_assembly    N75 113

  3. stats_assembly    N90 78

  4. stats_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软件也简单的给出了一个污染物种分布情况表,但是这个原理不一样。以后再讲


猜你喜欢

基因组 游记 | 工作资讯 

学习课程 | 好书分享 


菜鸟入门

Linux | Perl | R语言 | 可视化 

R包 | perl模块 | python模块


数据分析

ChIP-seq(上)ChIP-seq(下)RNA-seq | miRNA

WGS,WES,RNA-seq组与ChIP-seq之间的异同


编程实践

第0题 | 探索人类基因组序列 | 最后一题


直播基因组分析

我的基因组 | 解惑帖

一个标准的基因检测报告目录

生信技能树

直播我的基因组分析-目录-持续更新

【直播】我的基因组81:看看我的vcf文件的vaf分布情况

【直播】我的基因组82:如何对maf格式的突变文件统计vaf

批量IGV截图【直播】我的基因组83

从WGS测序得到的VCF文件里面提取位于外显子区域的【直播】我的基因组84

安装snpEFF工具并对VCF文件进行注释【直播】我的基因组85

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存