Python读取SAM文件
SAM（Sequence Alignment/Map）文件是用于存储高通量测序（HTS）数据的一种常见格式。在HTS数据分析中，我们经常会使用Python来处理和分析SAM文件。本文将介绍如何使用Python读取SAM文件，并展示一些实际应用的示例代码。
什么是SAM文件？
SAM文件是一种文本文件，用于存储序列对齐信息。它通常由测序仪生成的原始测序结果与参考基因组进行对比而得到。SAM文件中的每一行表示一个测序结果，包含了测序序列的名称、序列、质量分数和对齐位置等信息。
Python读取SAM文件的库
在Python中，我们可以使用pysam库来读取和处理SAM文件。pysam是一个常用的Python库，提供了一组函数和方法，用于处理SAM、BAM、CRAM和VCF等HTS文件格式。
首先，我们需要使用pip命令安装pysam库：
pip install pysam
安装完成后，我们就可以在Python脚本中导入pysam库并开始读取SAM文件了。
读取SAM文件
首先，我们需要使用pysam.AlignmentFile类来打开SAM文件。下面的示例代码展示了如何打开一个SAM文件并读取其中的一行数据：
import pysam
# 打开SAM文件
samfile = pysam.AlignmentFile("sample.sam", "r")
# 读取一行数据
alignment = samfile.readline()
# 打印对齐位置的参考序列名称和起始位置
print(alignment.reference_name, alignment.reference_start)
# 关闭SAM文件
samfile.close()
在上述示例中，我们首先使用pysam.AlignmentFile类打开了一个名为sample.sam的SAM文件，并指定模式为"r"，表示只读模式。
然后，我们使用readline方法读取了SAM文件中的一行数据，并将其存储在alignment变量中。通过访问alignment对象的属性，我们可以获取对齐位置的参考序列名称和起始位置。
最后，我们使用close方法关闭了SAM文件。
实际应用示例
接下来，让我们通过一个实际的示例来展示如何使用Python读取和处理SAM文件。
假设我们有一个SAM文件，其中存储了一系列测序结果。我们想要统计每个参考序列的对齐数量，并输出到一个文本文件中。
下面的示例代码展示了如何实现这一功能：
import pysam
# 打开SAM文件
samfile = pysam.AlignmentFile("sample.sam", "r")
# 统计每个参考序列的对齐数量
alignment_counts = {}
for alignment in samfile:
reference_name = alignment.reference_name
if reference_name in alignment_counts:
alignment_counts[reference_name] += 1
else:
alignment_counts[reference_name] = 1
# 输出结果到文本文件
with open("alignment_counts.txt", "w") as outfile:
for reference_name, count in alignment_counts.items():
outfile.write(f"{reference_name}\t{count}\n")
# 关闭SAM文件
samfile.close()
在上述示例中，我们首先创建了一个空字典alignment_counts，用于记录每个参考序列的对齐数量。
然后，我们使用for循环遍历SAM文件中的每一行数据。在循环中，我们首先获取当前对齐位置的参考序列名称，然后判断该参考序列是否已经存在于字典alignment_counts中。如果存在，则将对齐数量加一；如果不存在，则将对齐数量初始化为1。
接下来，我们使用with open语句打开一个名为alignment_counts.txt的文本文件，并将结果输出到该文件中。其中，outfile.write语句用于将参考序列名称和对齐数量写入文件。
最后，我们使用close方法关闭SAM文件。
通过上述示例代码，我们可以将每个参考序列的对齐数量统计结果输出到alignment_counts.txt文本文件中，供后续分析
