测序数据处理 —— 比对数据处理（Python）

最新推荐文章于 2024-10-13 21:12:05 发布

名本无名

最新推荐文章于 2024-10-13 21:12:05 发布

阅读量685

点赞数 5

分类专栏：测序数据处理文章标签： python linux 服务器

本文链接：https://blog.csdn.net/dxs18459111694/article/details/139209357

版权

测序数据处理专栏收录该内容

7 篇文章 0 订阅

订阅专栏

测序数据处理 —— 比对数据处理（Python）

文章目录

测序数据处理 —— 比对数据处理（Python）

介绍

前面详细介绍了比对数据的存储格式，下面我们介绍一下如何在 Python 中读取或修改比对结果。

pysam 是一个 Python 模块，它可以很容易地读取和操作存储在 SAM/BAM 文件中的比对结果。

它是 htslib C-API 的封装，提供对 SAM/BAM/VCF/BCF/BED/GFF/GTF/FASTA/FASTQ 文件的读写功能，以及访问 samtools 和 bcftools 软件包的命令行功能。

当前版本封装了 htslib-1.18、samtools-1.18 和 bcftools-1.18。使用 pip 进行安装

pip install pysam

SAM/BAM/CRAM

这三种文件的格式差不多，都是用来存储比对结果，其中 SAM 是文本类型文件，BAM 是二进制类型，CRAM 则是高度压缩格式，使用参考序列实现更好的压缩率。

我们一般会将结果保存为 BAM 文件。例如，读取 BAM 文件中的比对结果

import pysam

# 打开 BAM 文件
bamfile = pysam.AlignmentFile("example.bam", "rb")

# 遍历 BAM 文件中的所有 reads
for read in bamfile.fetch():
    print(read)

# 关闭 BAM 文件
bamfile.close()

注意

为了更快速地随机获取某些位置的比对结果，最好在使用之前为 BAM 文件创建索引。

AlignmentFile 对象提供了非常多的函数与方法，用于获取和修改比对结果信息。例如，检查是否为 BAM 文件构建了索引：

bamfile.check_index()  # True
# 或者
bamfile.has_index()    # True

如果没有找到索引文件（BAM 文件名称加 .bai 后缀）会抛出异常。

AlignmentHeader

BAM 文件的头部信息被封装为一个 AlignmentHeader 对象。例如

header = bamfile.header
header
# <pysam.libcalignmentfile.AlignmentHeader at 0x7f5ba306a450>

该对象包含很多属性和方法。例如

# 获取使用的参考基因组中染色体数目
header.nreferences
# 获取具体的染色体名称
','.join(header.references[:10])
# 'chr1,chr2,chr3,chr4,chr5,chr6,chr7,chrX,chr8,chr9'

# 获取染色体对应的 ID
header.get_tid('chr2')               # 1
# 获取 ID 对应的染色体名称
header.get_reference_name(1)         # 'chr2'
# 获取染色体的长度
header.get_reference_length('chr2')  # 243199373

AlignedSegment

具体的比对字段为 AlignedSegment 类型对象。该类存储了底层 C 中 samtools 比对信息结构的句柄。

对比对信息的访问会被转发到底层 C 结构，然后返回一个 python 对象。由于只转换了所需的数据，因此速度很快。

对于数据写入，底层的 C 结构会被原地更新。因为可变长度数据是连接在一起的，因此在更新字段时需要调整大小，所以这并不是建立比对结果的最有效方法。此外，BAM 文件中的比对结果可能会处于不一致的状态。

需要注意的一个问题是，序列应始终在质量分数之前设置。设置序列还会删除之前设置的任何质量分数。

对 AlignmentFile 对象迭代会返回 AlignedSegment 对象。例如

for read in bamfile:
    break
read.to_string()
# 'LH00190:192:22CNVYLT3:8:1145:40435:24441\t163\tchrM\t1\t2\t4M5I101M\t=\t1\t105\tCGATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTGGAGCC\tFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF5FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF\tAS:i:-45\tXS:i:-50\tXN:i:0\tXM:i:5\tXO:i:1\tXG:i:5\tNM:i:10\tMD:Z:0G0A0T0C68A32\tYS:i:-45\tYT:Z:CP'

可以获取比对信息。例如，获取比对位置

read.reference_name, read.reference_start, read.reference_end
# ('chrM', 0, 105)

获取 flag 值，以及几个根据 flag 值计算出的属性（is_xxxx 形式的属性）

print(read.flag)
print(read.is_read1, read.is_paired)
print(read.is_forward, read.is_reverse)
# 163
# False True
# True False

获取 CIGAR 信息

print(read.cigarstring)
print(read.cigartuples)
# 4M5I101M
# [(0, 4), (1, 5), (0, 101)]

其中 cigarstring 是由交替出现的整数和字符组成的字符串，表示操作的长度和类型，如果没有则返回 None。

cigartuples 结果以（操作、长度）元组列表的形式返回。对应表如下

`OP`	`Enum`	`ID`
`M`	`BAM_CMATCH`	`0`
`I`	`BAM_CINS`	`1`
`D`	`BAM_CDEL`	`2`
`N`	`BAM_CREF_SKIP`	`3`
`S`	`BAM_CSOFT_CLIP`	`4`
`H`	`BAM_CHARD_CLIP`	`5`
`P`	`BAM_CPAD`	`6`
`=`	`BAM_CEQUAL`	`7`
`X`	`BAM_CDIFF`	`8`
`B`	`BAM_CBACK`	`9`

获取所有 TAG

read.get_tags()
# [('AS', -77),
#  ('XN', 0),
#  ('XM', 4),
#  ('XO', 1),
#  ('XG', 18),
#  ('NM', 22),
#  ('MD', '0G1T0C68A54'),
#  ('YS', -15),
#  ('YT', 'CP')]

可以使用 get_tag 和 set_tag 获取和设置 TAG 的值。

print(read.get_tag('NM'))
print(read.set_tag('NM', 21))
print(read.get_tag('NM'))

比对到负链的 reads 的序列会以反向互补方式存储在 BAM 文件中。可以使用 get_forward_sequence 可以将反向互补的读数返回到其原始序列。

for read in bamfile:
    if read.is_reverse:
        break
read.seq
# CGATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTGGAGCC
read.get_forward_sequence()
# GGCTCCAGCGTCTCGCAATGCTATCGCGTGCACACCCCCCAGACGAAAATACCAAATGCATGGAGAGCTCCCGTGAGTGGTTAATAGGGTGATAGACCTGTGATCCATCG

计算比对到参考基因组上的长度

print(read.reference_length)
assert read.reference_length == read.reference_end - read.reference_start
# 105

计算查询序列的长度

read.query_length

该值与 BAM/SAM 文件中提供的序列长度一致，长度包括 softclip，等于 len(query_sequence)。

如果 BAM/SAM 文件中没有序列，则查询长度为 0。在这种情况下，可从使用 infer_query_length 方法从 CIGAR 比对信息中推断出读长。

还有很多方法和属性，可以自己探索一下。

fetch

fetch 函数用于获取覆盖基因组区域内比对上的 reads。

其中参数 contig，start，stop 以及 region 和 tid 用于指定基因组区域。 start 和 stop 表示以 0 为基础的半开放区间。为了向后兼容，reference 和 end 也可分别作为 contig 和 stop 的同义词。

如果缺少任何坐标，它们将被最小坐标（起始坐标）或最大坐标（停止坐标）取代。

请注意，region 字符串是以 1 为基础的，而 start 和 stop 表示以 0 为基础、半开放坐标的区间（如 BED 文件和 Python 切片）。

for read in bamfile.fetch(region='chrM', start=100, stop=105):
    print(read.reference_name)
    print(read.reference_start, read.reference_end)
    break
# chrM
# 0 150

如果 contig 或 region 设置为 '*' ，则返回 BAM 文件末尾的未比对上的 reads。设置为 '.' 则将从文件开头开始迭代。

for read in bamfile.fetch(region='*'):
    print(read.qname)
    break
# LH00190:192:22CNVYLT3:8:1101:33574:1753

如果没有设置 contig 或 region，将获取文件中的所有映射读数。读数将按参考序列顺序返回，但不一定是文件中的顺序。这种迭代模式仍然需要索引。如果没有索引，要设置参数 until_eof=True。

如果只设置了 contig，则将获取比对上 contig 的所有 reads。

set(read.reference_name for read in bamfile.fetch('chrM'))
# {'chrM'}

SAM 文件不允许随机访问。如果设置 region 或 contig 参数，则会引发异常。

前面提到的 contig 指代的可以是染色体或者一些特殊的片段。

pileup

pileup 方法是一种用于生成比对数据的堆积覆盖（pileup）表示的方法。可以表示在每个参考基因组位置上所有 reads 的累积信息。有助于分析序列比对的覆盖度、变异、和其他重要的信息。

类似于 fetch，可以指定基因组区域。例如

for pileupcolumn in bamfile.pileup("chrM", 100, 150):
    print(f"Coverage at base {pileupcolumn.reference_pos}: {pileupcolumn.nsegments}")

    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            print(f"\tbase in read {pileupread.alignment.query_name}: {pileupread.alignment.query_sequence[pileupread.query_position]}")
# Coverage at base 0: 8
# 	base in read LH00190:192:22CNVYLT3:8:2175:37883:19634: C
# 	base in read LH00190:192:22CNVYLT3:8:2181:3848:26860: A
# 	base in read LH00190:192:22CNVYLT3:8:2252:36385:4733: C
# 	base in read LH00190:192:22CNVYLT3:8:1145:40435:24441: C
# Coverage at base 1: 10
# 	base in read LH00190:192:22CNVYLT3:8:1134:41979:7217: A
# 	base in read LH00190:192:22CNVYLT3:8:1202:19141:11479: A
# 	base in read LH00190:192:22CNVYLT3:8:2175:37883:19634: A
# 	base in read LH00190:192:22CNVYLT3:8:2181:3848:26860: A
# 	base in read LH00190:192:22CNVYLT3:8:2252:36385:4733: A
# 	base in read LH00190:192:22CNVYLT3:8:1145:40435:24441: G
# 	base in read LH00190:192:22CNVYLT3:8:2119:47980:20948: A

pileup 方法生成一个 PileupColumn 类型的对象迭代器，每个对象代表一个参考基因组位置的覆盖情况。

其中，reference_pos 返回该位置所在基因组的位置；nsegments 表示有多少 reads 覆盖到该位置；pileups 返回比对到该位置的 reads 列表。

而每个 PileupColumn 包含多个 PileupRead 类型的对象，代表在该位置上的 reads。

其中 is_del 和 is_refskip 用于判断 reads 中该位置的碱基是否处于删失或 CIGAR 操作符为 N 的位置上。alignment 则返回 AlignedSegment 类型的对象。

还有 is_head 和 is_tail 用于判断该位置的碱基是不是 read 上的最左或最右侧的碱基。

count

count 用于计算基因组区域内的 reads 数量。例如

# 计算整个文件的 reads 数
total_reads = bamfile.count()
print(f"Total reads: {total_reads}")

# 计算特定染色体（例如 chr1）的 reads 数
chr1_reads = bamfile.count(contig="chr1")
print(f"Reads on chr1: {chr1_reads}")

# 计算特定区域（例如 chr1:10003-10023）的 reads 数
region_reads = bamfile.count(contig="chr1", start=10003, stop=10023)
print(f"Reads on chr1:10003-10023: {region_reads}")
# Total reads: 1661102
# Reads on chr1: 43494
# Reads on chr1:10003-10023: 3

count_coverage

count_coverage 方法用于计算每个碱基位置的覆盖度，并返回覆盖度的数组。例如

# 计算特定区域的覆盖度（例如 chr1:10003-10023）
coverage = bamfile.count_coverage('chr1', start=10003, stop=10023)

# 输出 A、C、G、T 四个碱基的覆盖度
a_coverage = coverage[0]
c_coverage = coverage[1]
g_coverage = coverage[2]
t_coverage = coverage[3]

print(f"A coverage: {a_coverage}")
print(f"C coverage: {c_coverage}")
print(f"G coverage: {g_coverage}")
print(f"T coverage: {t_coverage}")
# A coverage: array('L', [0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 3, 3, 0, 0, 0, 0, 2, 3, 0, 0])
# C coverage: array('L', [2, 2, 2, 0, 0, 0, 3, 3, 2, 0, 0, 0, 3, 2, 2, 0, 0, 0, 3, 2])
# G coverage: array('L', [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
# T coverage: array('L', [0, 0, 0, 2, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0])

get_index_statistics

get_index_statistics 用于返回索引中存储的每个染色体的比对/未比对上 reads 的统计信息，与 samtools idxstats 命令打印的统计信息类似。

# 获取索引统计信息
index_stats = bamfile.get_index_statistics()

# 输出索引统计信息
for stat in index_stats:
    print(f"Reference: {stat.contig}, Mapped reads: {stat.mapped}, Unmapped reads: {stat.unmapped}")
    break
# Reference: chr1, Mapped reads: 43494, Unmapped reads: 0

更多函数的使用方式，可以参考官方文档

pysam 是一个强大的库，不仅可以处理 BAM/SAM 文件，还可以用于读取和处理 FASTA 和 FASTQ 数据。

FASTA

pysam 提供了 FastaFile 类用于对 FASTA 格式的文件进行读写。例如

import pysam

# 打开 FASTA 文件
fasta_file = pysam.FastaFile("example.fasta")

# 获取序列名称列表
references = fasta_file.references
print(f"References: {references}")

# 读取特定序列
seq_name = references[0]  # 假设我们读取第一个序列
sequence = fasta_file.fetch(seq_name)
print(f"Sequence for {seq_name}: {sequence}")

# 读取特定范围的序列
start = 100
end = 200
sub_sequence = fasta_file.fetch(seq_name, start, end)
print(f"Sub-sequence for {seq_name} ({start}-{end}): {sub_sequence}")

# 关闭 FASTA 文件
fasta_file.close()

FASTQ

pysam 提供了 FastxFile 类用于对 FASTQ 和 FASTA 文件进行读写。例如

import pysam

# 打开 FASTQ 文件
fastq_file = pysam.FastxFile("example.fastq")

# 遍历 FASTQ 文件中的 reads
for entry in fastq_file:
    print(f"Name: {entry.name}")
    print(f"Sequence: {entry.sequence}")
    print(f"Quality: {entry.quality}")

# 关闭 FASTQ 文件
fastq_file.close()