samtools建立fasta索引

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了samtools建立fasta索引相关的知识,希望对你有一定的参考价值。

参考技术A 本文摘抄自: https://www.omicsclass.com/article/42
侵删!
fasta是常用的序列存储格式,很多软件(如GATK、IGV等)在导入序列以及进行快速查找时通常需要建立索引文件。下面就来介绍如何使用 samtools 便捷的建立fasta文件的索引以及快速进行序列提取。

1 建立索引

建立索引只需在Linux下输入命令:samtools faidx input.fa

这里序列文件为 input.fa,生成的索引文件以 .fai 结尾。需要注意的是,输入的fasta文件的每条序列除最后一行外,其余行的长度必须相同,否则会报错哦!最后生成的.fai文件如下, 共5列,以制表符分隔;

第一列 NAME : 序列的名称,只保留“>”后,第一个空白之前的内容;

第二列 LENGTH : 序列的长度,单位为bp;

第三列 OFFSET : 第一个碱基的偏移量,从0开始计数,换行符也统计进行;

第四列 LINEBASES : 除了最后一行外, 其他代表序列的行的碱基数, 单位为bp;

第五列 LINEWIDTH : 行宽, 除了最后一行外, 其他代表序列的行的长度,包括换行符,在windows系统中换行符为\r\n,要在序列长度的基础上加2。

2 提取序列

除建立索引外,还可以利用samtools方便的提取序列,例如:

samtools faidx input.fa chr2 > chr2.fa,会得到含chr2这条序列的fasta格式的文件,如果是多条序列,只需在文件后罗列需提取的序列ID即可,使用空格分隔,如 samtools faidx input.fa chr1 chr2 chr3 > chr.fa。

再如:samtools faidx input.fa chr2:1-1000 > chr2.fa,能得到chr2序列的第1到第1000个碱基的fasta格式的文件,同样可以提取多条序列。

samtools 安装

2samtools-faidx index

1、samtools faidx

1、samtools faidx 能够对fasta 序列建立一个后缀为.fai 的文件根据这个.fai 文件和原始的fasta文件, 能够快速的提取任意区域的序列

2、用法:samtools faidx genome.fa      #生成genome.fa.fai

3、例子

该命令对输入的fasta序列有一定要求:对于每条序列,除了最后一行外, 其他行的长度必须相同,  

>one 
ATGCATGCATGCATGCATGCATGCATGCAT 
GCATGCATGCATGCATGCATGCATGCATGC 
ATGCAT 
>two another chromosome 
ATGCATGCATGCAT 
GCATGCATGCATGC 

最后生成的.fai文件如下, 共5列,\t分隔;

one 66 5 30 31
two 28 98 14 15

第一列 NAME   :   序列的名称,只保留“>”后,第一个空白之前的内容;

第二列 LENGTH:   序列的长度, 单位为bp;

第三列 OFFSET :   第一个碱基的偏移量, 从0开始计数,换行符也统计进行;

第四列 LINEBASES : 除了最后一行外, 其他代表序列的行的碱基数, 单位为bp;

第五列 LINEWIDTH :  除了最后一行外, 其他代表序列的行的长度, 包括换行符, 在windows系统中换行符为\r\n, 要在序列长度的基础上加2;

提取序列:

samtools faidx input.fa               #生成索引input.fa.fai

samtools faidx input.fa chr1 > chr1.fa        #提取chr1序列

samtools faidx input.fa chr1:1-10000 > chr1.fa  #提取chr1序列上1-10000间的序列

2、samtools index

1、

  samtools index accepted_hits.bam            #生成索引accepted_hits.bam.bai

  samtools view accepted_hits.bam contig1         #提取比对到chr1序列reads

  samtools view accepted_hits.bam contig:1-10000    #提取比对到chr1序列上100-200区间的reads

2、

  samtools tview accepted_hits.bam ../genome.fa    #samtools tview运用要求要先对bam文件index索引

 

以上是关于samtools建立fasta索引的主要内容,如果未能解决你的问题,请参考以下文章

samtools的用法简介

2samtools-faidx index

区别samtools faid产生的.fai文件功能和bwa index 产生的四个文件的功能

利用samtools对参考基因组构建索引fai文件

sam文件转换为bam文件——SAMtools

tabix 操作VCF文件