pyfaidx：序列处理工具

pyfaidx是一个处理序列文件的软件，类似于samtools中的faidx。

参考地址：https://github.com/mdshw5/pyfaidx

安装：

pip install pyfaidx

基本方法：

读取文件：

参数简介：

split_char：分隔符

duplicate_action：指定对于重复的序列标识符该如何处理，默认值为 "stop"，即在遇到重复标识符时停止读取

key_function：传递一个函数作为参数来处理序列标识符

as_raw：是否以字符串形式打印

sequence_always_upper：是否总是以大写打印

```
from pyfaidx import Fasta
genes = Fasta('test.fasta')
```

打印所有的id：

```
genes.keys()
```

打印id和它对应的序列，可以用切片方法截取序列，也可以在截取的序列中继续截取：

```
genes['xxxxx'][100:200]
genes['xxxxx'][100:200][:10]
genes['xxxxx'][100:200][::-1]
```

打印序列：

```
genes['xxxxx'][100:200].seq
```

打印id：

```
genes['xxxxx'][100:200].name
```

打印起点：

```
genes['xxxxx'][100:200].start
```

打印终点：

```
genes['xxxxx'][100:200].end
```

打印id和起始点、终止点：

```
genes['xxxxx'][100:200].fancy_name
```

打印长度：

```
len(genes['xxxxx'])
```

打印互补序列：

```
genes['xxxxx'][100:200].complement
```

打印反向序列：

```
genes['xxxxx'][100:200].reverse
```

用get_seq方法来打印信息：

```
genes.get_seq('xxxxx', 201, 210)
```

打印反向互补序列：

```
genes.get_seq('xxxxx', 201, 210, rc=True)
```

用get_spliced_seq方法截取序列：

```
segments = [[1, 10], [20, 30]]
genes.get_spliced_seq('xxxxx', segments)
```

迭代读取，打印id和完整的信息

```
genes = Fasta('test.fasta')
for record in genes:
    print(record.name)
    print(record.long_name)
```

作为numpy数组有效地访问：

```
from pyfaidx import Fasta
import numpy as np
genes = Fasta('test.fasta')
np.asarray(genes['xxxxx'])
```

序列可以使用预读缓冲区在内存中缓冲，以实现快速顺序访问：

```
from timeit import timeit
fetch = "genes['xxxxx'][100:200]"
read_ahead = "import pyfaidx; genes = pyfaidx.Fasta('test.fasta', read_ahead=10000)"
no_read_ahead = "import pyfaidx; genes = pyfaidx.Fasta('test.fasta')"
string_slicing = "genes = {}; genes['xxxxx'] = 'N'*10000"

timeit(fetch, no_read_ahead, number=10000)
timeit(fetch, read_ahead, number=10000)
timeit(fetch, string_slicing, number=10000)
```

如果想要就地修改FASTA文件的内容，可以使用mutable参数。FastaRecord的任何部分都可以用相同长度的字符串替换，立即并永久地改变文件的内容：

```
genes = Fasta('test.fasta', mutable=True)
# 永久更改前10个为N
genes['xxxxx'][:10] = 'NNNNNNNNNN'
```

FastaVariant类提供了一种集成单核苷酸变体调用以生成一致序列的方法：

```
# het：杂合，hom：纯合，sample：样本名，call_filter='GT == "0/1"'：只保留杂合子
consensus = FastaVariant('chr22.fasta', 'chr22.vcf.gz', sample='xxxxx', het=True, hom=True, call_filter='GT == "0/1"')
# 打印索引变异位点
consensus['22'].variant_sites
# 打印位点对应序列
consensus['22'][16042790:16042800]

# 其他方法
Fasta('tests/data/chr22.fasta')['22'][16042790:16042800]
```

还可以使用pathlib指定路径：

```
from pyfaidx import Fasta
from pathlib import Path
genes = Fasta(Path('test.fasta'))
```

使用 fsspec 和 pyfaidx 从文件系统中读取 fasta 文件：

fsspec 是一个通用的文件系统规范抽象库，它为 Python 提供了一个简单且一致的 API，用于访问各种不同类型的文件和存储后端（如本地文件系统、S3、HDFS 等）。使用 fsspec 库，可以轻松地在不同类型的文件系统之间进行切换，同时还能够实现高效的文件系统操作和数据处理。

```
import fsspec
from pyfaidx import Fasta
of = fsspec.open("s3://broad-references/hg19/v0/Homo_sapiens_assembly19.fasta", anon=True)
genes = Fasta(of)
```

这个软件还提供了命令行软件：

```
# 可同时看多个id的id和序列aaa、bbb
faidx test.fasta aaa bbb

# 查看全名
faidx --full-names test.fasta aaa

# 只打印序列
faidx --no-names test.fasta aaa

# 打印互补序列
faidx --complement test.fasta aaa

# 打印反向序列
faidx --reverse test.fasta aaa

# 打印反向互补序列
faidx --reverse --complement test.fasta aaa

# --lazy 参数来对文件进行索引
faidx --lazy test.fasta aaa

# 根据上面索引不到则打印出替代的字符
faidx --lazy --default-seq='xxx' xxxxx>xxxxx

# 打印出id和长度（dp），第一列为id名，第二列为其长度
faidx --transform chromsizes test.fasta

# 打印出bed格式的文件
faidx --transform bed test.fasta 

# 打印出id名、起始位点、终止位点、A、T、C、G、N的统计信息
faidx --transform nucleotide test.fasta

# 对每条染色体/序列的碱基序列进行转置后输出
faidx --transform transposed test.fasta

# 将所有序列拆开保存为文件，并且以id名命名
faidx --split-files test.fasta

# --delimiter='_ ' ：使用下划线和空格作为序列名中的分隔符，最后的参数 000465.3 是你想要提取的染色体/序列的名称，即序列名中下划线和空格之后的部分。执行该命令后会输出指定序列的碱基序列
faidx --delimiter='_' test.fasta 000465.3

# --size-range 选项，表示只提取长度在 5500 到 6000 之间的染色体/序列，-i chromsizes 选项表示输出格式为 chrom.sizes 格式，每行第一列为符合条件的染色体/序列的名称，第二列为其长度（单位为 bp）
faidx --size-range 5500,6000 -i chromsizes test.fasta

# 根据指定的 BED 文件提取 FASTA 文件中对应的区域，并输出这些区域的碱基序列，-m 选项，表示输出时不会在序列间插入空行
faidx -m --bed regions.bed test.fasta

# -M 选项表示输出时插入的空行数
faidx -M --bed regions.bed test.fasta

# 去除扩展名后面的数字
faidx -e "lambda x: x.split('.')[0]" test.fasta -i bed
```