软件介绍之Hisat2

咱们《生信技能树》的B站有一个lncRNA数据分析实战,缺乏配套笔记,所以我们安排了100个lncRNA组装案例文献分享,以及这个流程会用到的100个软件的实战笔记教程

下面是100个lncRNA组装流程的软件的笔记教程

一、Hisat2介绍

Hisat是一种高效的RNA-seq实验比对工具。它使用了基于BWT和Ferragina-manzini (Fm) index 两种算法的索引框架。使用了两类索引去比对,一类是全基因组范围的FM索引来锚定每一个比对,另一类是大量的局部索引对这些比对做快速的扩展。比对原理可阅读文献原文:HISAT: a fast spliced aligner with low memory requirements.

二、Hisat2设计原则

1.优化了索引建立策略

hisat2应用了基于bowtie2的方法去处理很多低水平的用于构建和查询FM索引的操作。但是与其它比对器不同的是,该软件应用了两类不同的索引类型:代表全基因组的全局FM索引和大量的局部小索引,每个索引代表64000bp的区域。以人类基因组为例,创建了约48000个局部索引,每一个索引与其相邻索引重叠1024bp,最终可以覆盖这个30亿个碱基的基因组。这种存在交叉(overlap)的边界可以轻松的比对那些跨区域的read(可变剪切体)。尽管有很多索引,但是hisat会把他们使用合适方法压缩,最终只占4GB左右的内存。

2.采用了新的比对策略

RNA-seq产生的reads可能跨越长度比较大的内含子,在哺乳动物中甚至最长能达到1MB,同时外显子比较短,read也比较短。当使用100-bp reads时,会有很多read(模拟数据中大概34%)跨两个外显子的情况。为了更好的比对,将跨外显子的reads分成了三类:1)长锚定read,两个外显子中的每个都至少有16 bp;2)中间锚定read,一个外显子中具有8-15bp;3)短锚定read,仅与其中一个外显子比对上1-7bp。所以总的reads可以被划分为五类:1)不跨外显子的read 2)长锚定read 3)中间锚定read 4)短锚定read 5)跨两个外显子以上的read。在模拟的数据中,有25%左右的read是长锚定read,这种read在大多数情况下可以被唯一的定位到人的基因组上。5%为中间锚定read,对于这类,很多依赖于全局索引的算法就很难执行下去(需要比对很多次),而hisat,可以先将read中的长片段实现唯一比对,之后再使用局部索引对剩下的小片段进行比对(局部索引可以实现快速检索)。4.2%为短锚定read,因为这些序列特别短,因此只能通过在hisat比对其它read时发现的剪切位点或者用户自己提供的剪切位点来辅助比对。最后还有3%的是跨多个外显子的read,比对策略在hisat的online method中有介绍,文章中没有详解。比对过程中,中间锚定read、短锚定read、跨多个外显子read的比对占总比对时长的30%-60%,而且比对错误率很高!

三、软件安装

conda安装

conda install hisat2

四、hisat2 index建立

1.直接下载

直接在网站http://daehwankimlab.github.io/hisat2/download/#index下载已建好的index

wget -c https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz

2.自己构建index

hisat2可以自己建立index文件,如果没有现成的index的话,那就得自己去建立了,这个就比较麻烦而且耗内存和时间

# 其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites,但是先做的目的是为了提高运行效率
extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &
extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &
# 建立index, 必须选项是基因组所在文件路径和输出的前缀
hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19

五、Hisat2的用法

安装完成以后,可以使用hisat2 -h来查看软件的帮助文档。

  1. 软件用法:

2.常用参数:

-x :索引的前缀,也可以添加路径,例如 ~/genome/homo_GRCh38.dna.primary

-1 :双端测序的第一个文件,若有多组,则用逗号隔开,名字与2要匹配,如-1 flyA_1.fq,flyB_1.fq

-2 :类比上,名字与1要匹配,如-2 flyA_2.fq,flyB_2.fq

-U :单端数据文件

-S :保存到的sam文件,默认输出到标准输出

-p :线程数

--dta  :在下游使用stringtie组装的时候量身定做的参数。使用此选项,
       :HISAT2需要更长的锚长度才能从头发现拼接位点,这样可以减少与短锚的对齐,
       : 这有助于转录汇编程序显着提高计算和内存使用率。
       : 当然下游不使用stringtie也可以使用此参数减少短锚定read
       
--dta-cufflinks :下游使用cufflinks则需要添加此参数

-q :指定读取的测序文件是fastq文件(含有质量信息)

-f :指定读取的测序文件是fa文件

-t :打印加载索引文件和对齐读数所需的时钟时间

--phred33 :质量表示方法,如果使用fq文件,建议选择,同理有小众的--phred64

--min-intronlen :设置最小内含子的长度,默认值20

--max-intronlen :设置最大内含子长度,默认500000

--known-splicesite-infile :提供已知的剪接位点列表,HISAT2利用这些列表比对,可以用 hisat2_extract_splice_sites.py和gtf生成信息

--novel-splicesite-outfile :在结果中生成一个剪接位点的列表

--novel-splicesite-infile : 可以使用novel-splicesite-outfile所生成的剪接列表作为
                          : 新剪接列表,应该可以自己定义。为特定基因。

三、软件运行命令

hisat2 输出文件是sam格式,可通过管道符与Samtools工具连用,直接生成bam,并对bam文件进行sorted以方便后续数据处理

hisat2 -p 2 -x ${index} -1 $fq1 -2 $fq2 | samtools sort -@ 6 -o ./${sample}.bam -

批量运行脚本

mkdir 04.mapping
ls /home/data/lihe/lncRNA_project/03.trim/*_1.fq.gz > 1
ls /home/data/lihe/lncRNA_project/03.trim/*_2.fq.gz > 2
ls /home/data/lihe/lncRNA_project/03.trim/*_1.fq.gz | cut -d "/" -f 7 | cut -d "_" -f 1 > 0
paste 0 1 2 > config

cd 04.mapping

cat > 04.hg38_mapping.sh
index=/home/data/server/reference/index/hisat/hg38/genome
config=$1
number1=$2
number2=$3
cat $1 | while read id
do
    if((i%$number1==$number2))
    then
    arr=(${id})
    fq1=${arr[1]}
    fq2=${arr[2]}
    sample=${arr[0]}
    hisat2 -p 2 -x ${index} -1 $fq1 -2 $fq2 | samtools sort -@ 6 -o ./${sample}.bam -
    fi    ## end for number1
    i=$((i+1))
done

for i  in {0..6}
do 
(nohup bash 04.hg38_mapping.sh  config 7 $i 1>log.$i.txt 2>&1 & )
done 

四、结果解读

17968900 reads; of these:    #reads数目
  17968900 (100.00%) were paired; of these:   #能配对的reads,由于用的是trimmomatic后的paired clean data,故100%
    4434269 (24.68%) aligned concordantly 0 times    #不能比对到基因组上
    10748474 (59.82%) aligned concordantly exactly 1 time   #比对到基因组上,且唯一匹配,即unique map
    2786157 (15.51%) aligned concordantly >1 times  #比对到基因组上,但结果>1次,即不唯一匹配
    ----
    4434269 pairs aligned concordantly 0 times; of these:  #上述不能比对上的reads中
      269467 (6.08%) aligned discordantly 1 time #这些虽然比对上了,但是不合理,例如方向、R1 R2之间的距离等
    ----
    4164802 pairs aligned 0 times concordantly or discordantly; of these:   #剩余的不匹配reads中
      8329604 mates make up the pairs; of these:
        5311882 (63.77%) aligned 0 times
        2073153 (24.89%) aligned exactly 1 time   #作为单端来比对,这些可以唯一比对
        944569 (11.34%) aligned >1 times   #这些不唯一比对
85.22% overall alignment rate  #结果是唯一比对+不唯一比对+不能比对中可以单端比对的

五、导入IGV可视化

使用samtools 建立bam文件的index,然后可把bam文件导入igv可视化

ls /home/data/lihe/lncRNA_project/04.mapping/*.bam > 1
ls /home/data/lihe/lncRNA_project/04.mapping/*.bam | cut -d "/" -f 7 | cut -d "_" -f 1 | cut -d "." -f 1 > 0
paste 0 1 > config

cat > index.sh

config=$1
number1=$2
number2=$3
cat $1 | while read id
do
    if((i%$number1==$number2))
    then
    arr=(${id})
    input=${arr[1]}
    sample=${arr[0]}
    samtools index -@ 4 ${input} ./${sample}.bam.bai
    fi    ## end for number1
    i=$((i+1))
done

for i  in {0..4}
do 
(nohup bash index.sh  config 5 $i 1>log.$i.txt 2>&1 & )
done

文末友情推荐

与十万人一起学生信,你值得拥有下面的学习班:

(0)

相关推荐

  • NGS数据分析实践:06. 数据预处理 - 序列比对+PCR重复标记+Indel区域重比对+碱基质量重校正

    NGS数据分析实践:06. 数据预处理 - 序列比对+PCR重复标记+Indel区域重比对+碱基质量重校正 目录 1. 序列比对 1.1 参考基因组建索引 1.2 序列比对 2. 排序 3. PCR重 ...

  • 转录组学习一(软件安装)

    开篇:2017/10/7正式开始生信技能树论坛里的转录组入门从Linux下软件的安装 到 差异表达基因的功能注释及功能分析相关. 转录组相关软件的安装 任务 本地Windows电脑及服务器Linux系 ...

  • 比对软件-Bowtie2

    bowtie2 语法很重要!!!! Usage: bowtie2 [options]* -x <index> {-1 <m1> -2 <m2> | -U <r ...

  • Hisat2-Align 完成!-Win / Mac, 成功伴随失败。

    前两天推了下 Hisat2-build 插件,目前尚未对任何人开放.事实上,同一天,我也完成了 Hisat-Align,也就是,可以直接在 Windows 下跑Hisat2了.最近事情比较多,我抽空试 ...

  • 用了这么久的bowtie2,却不知道它的结果含义?

    今天是生信星球陪你的第529天 大神一句话,菜鸟跑半年.我不是大神,但我可以缩短你走弯路的半年~ 就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~ 这里有豆豆和花花的学习历程,从新手 ...

  • 转录组学习五(reads比对)

    转录组学习一(软件安装)转录组学习二(数据下载)转录组学习三(数据质控)转录组学习四(参考基因组及gtf注释探究)转录组学习五(reads的比对与samtools排序)转录组学习六(reads计数与标 ...

  • 多种组态软件介绍及优缺点

    多种组态软件介绍及优缺点

  • 学术干货 | Endnote 系列之软件介绍及基本功能演示(附带大招)

    文献管理还在用文件夹?那你就out了! 小编之前听同学讲过一个故事:某"老板"自从被学生教会了Endnote之后,每次实验室来了新学生,都要说:"嘿嘿嘿,Endnote好 ...

  • 音乐人练功 5个实用网站/软件介绍

    写歌遇到困难的话可以先缓一缓,说不定网络上有好用工具可以帮你解决难关 想让自己更加进步,需要反复的练习与改进 做音乐亦是如此,除了知识更需要实战练习! 或许你没办法长时间在工作桌前制作音乐,让自己累积 ...

  • 2009年:国内常见的PKM工具软件介绍

    PKM 2   PKM2 是基于内容的个人知识管理系统,它可以将您看到的所有文字.图片信息全部转储为 HTML 格式文档保存到数据库中.这些信息包括:你的笔记.网上的网页内容.本地机器里的文档内容.P ...

  • 【Affinity Photo基础教程】软件介绍

    欢迎来到Affinity Photo,如果你是刚开始使用Affinity Photo,这个教程来指导你通过一些基本的步骤,并给你一个想法,你可以开始用这个软件实现.当打开时Affinity Photo ...

  • 开源SlicerSALT软件介绍

    Slicer Shape AnaLysis Toolbox  (SlicerSALT) 是一款面向全球开源.免费的综合软件,其核心和Gui基于3DSlicer开发,可帮助生物医学科学家在其成像研究中精 ...

  • 通用脑-机接口研究软件介绍

    更多技术,第一时间送达 g.BCIsys-g.tec的脑机接口研究环境 g.BCIsys-g.tec基于Matlab/Simulink的系统 g.tec提供了完整的基于MATLAB的研发系统,包括数据 ...

  • 菱机妙用 | 伺服系统 Motorizer软件介绍及滚珠丝杆模块选型方法

    本 | 期 | 看 | 点 伺服系统 Motorizer软件整体介绍 伺服系统 Motorizer滚珠丝杠模块选型介绍 看点一 伺服系统 Motorizer软件整体介绍 点击视频快速掌握设置方法 灵活 ...

  • 华为MDC智能驾驶软件介绍

    平台软件零层逻辑架构如下图,可以看出,由基础层.功能层.应用层和服务层组成. 零层逻辑架构 从平台软件一层逻辑架构可以看出,MDC用了华为自研的越影操作系统.兼容Autosar标准的软件中间件,提供完 ...