samtools-faidx用法

samtools faidx xxx.fa

能够对fasta 序列建立一个后缀为.fai 的文件

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

>one

ATGCATGCATGCATGCATGCATGCATGCAT

GCATGCATGCATGCATGCATGCATGCATGC

ATGCAT

>two another chromosome

ATGCATGCATGCAT

GCATGCATGCATGC

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

one 66 5 30 31

two28981415

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

第二列 LENGTH:   序列的长度， 单位为bp；

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

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

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

提取序列：

# 

对于UCSC的chr肯定是可以用的：

>chr1

>chr2

...

对于ensemble呢？可行

>1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF

>2 dna:chromosome chromosome:GRCh38:2:1:242193529:1 REF

...

# head -n  2  chr1.fa

>1

ATCG...

提取all：

对于提取序列来说，我们一般会选择用samtools、bedtools等，因为他们有非常完善的 *** 作 序列 的功能，比如取两个bed文件的交集，提取某段序列等等。但实际上，他们对于 区间 的 *** 作是不如R方便的，比如取基因body，取第一个外显子，取启动子区域等等。所以，我们有时候就会想能否在R里面，花式找了些区间之后，提取他们的序列，实际上，这也是可以的，需要分成几步。

首先我们需要先加载一些包

然后是提取基因的区间

然后根据基因区域来提取启动子序列

这里就出现了一个坑了，就是你在用Txdb包的时候，出现了warning，但用自建gtf构建的Txdb没有出现了warning。实际上，这种warning是因为你的区间取过界了。

而promoter_gtf没有报错是因为他不知道 你的染色体的长度范围在哪里 。

所以下面我的 *** 作都是基于promoter_gtf，这样还可以应用在 非模式物种 里面。首先你在构建Txdb的时候就需要给他传进去染色体长度了。但你怎么知道每条染色体是多长呢，这就需要你的fai文件了。fai文件是给基因组fa做索引的时候出现的，没有的话你可以直接用samtools faidx做下。

然后我们就可以在构建的时候手动输进对应的长度了，但这只是拟南芥这个组装完整的基因组，如果是一些非模式物种，scaffold级别的话，就会有数百条染色体长度了，根本不可能手动输。所以我们还需要读进R里面进行一些 *** 作。

这样我们提取启动子序列的时候就会报错了

然后我们只取前2个基因出来看看效果

本文摘抄自： 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 安装

---