浅析samtools mpileup命令的使用及其对应的Python使用接口（一）（小白友好）
来自： 更新日期：早些时候
~
作为一名新转行至生物信息学领域的初学者，我每天都在与各种文件格式及其对应的工具进行痛苦的斗争。最近，在阅读一些大神的文章时，我学到了一些有用的工具，其中之一就是samtools mpileup命令。由于之前没有使用过这个工具，我借助了samtools mpileup的帮助文档并结合一个小的数据集来学习它的用法。同时，作为坚定的Python拥护者，我还学习了Pysam包中这一功能的使用方法，这将在后续的文章中进行详细介绍。
本文将通过实际案例来直观展示samtools mpileup命令的运行结果，并尝试详细解释一些细节。首先，我选择了10x Genomics的1k PBMC数据集作为测试数据，你可以通过以下shell命令下载这个数据集：
为了使运行结果更加直观，我从下载到的BAM文件中提取了几个比对记录，并将它们整合到一个新的文件中。这个过程包括以下两步：
此外，我还单独下载了人的一号染色体的参考序列，以便于后续分析。现在，让我们来详细了解samtools mpileup命令的基本概念和使用方法。
pileup命令意为“堆叠”，它能够在给定基因组区间内的每个碱基位置（column）上，将映射到该位置上的所有reads堆叠起来，以展示输入BAM文件在指定基因组区间内的每个碱基位置的总体信息。
要开始使用samtools mpileup，首先需要设置基因组区间。可以通过-r/--region参数来指定。例如：
请注意，samtools mpileup默认会忽略一些alignment，包括UNMAP、SECONDARY、QCFAIL和DUP，这些alignment的FLAG值已被标记为256，即表示次要比对。如果你的示例BAM文件中包含这些类型的alignment，那么它们不会被输出。你可以通过重新设置参数--ff UNMAP,QCFAIL,DUP来显示这些被标记为次要比对的alignment。
输出结果中的字段分别代表以下内容：1）参考序列名称；2）基因组位置坐标；3）参考基因组上该位置的碱基；4）覆盖该基因组位置的reads数量；5）每个reads在该位置的碱基；6）碱基的质量分数。
同时，samtools mpileup命令能够正确识别并忽略reads开头被soft clip的一个碱基。此外，你还可以通过-l/--positions参数接BED文件来指定基因组区间，如果同时使用了-r和-l，那么最终生效的区间是两者交集。
在使用samtools mpileup时，还可以设置过滤条件以筛选alignment。例如，可以设置最低的Mapping Quality阈值（-q --min-MQ）和碱基质量分数阈值（-Q --min-BQ），或者排除特定的Read Group（-G --exclude-RG）。
为了更详细地分析结果，samtools mpileup命令在运行时提供了输出参考基因组序列的功能，通过-f / --fasta-ref参数来提供参考序列Fasta文件。这将改变输出结果中某些字段的表示方式。
在探索samtools mpileup的使用方法时，我尽量覆盖了这一命令的大部分用法。对于更深入的细节和使用技巧，我鼓励小伙伴们参考官方文档进行学习。同时，如果你在使用过程中发现任何不妥之处，欢迎提出宝贵的意见和建议。
本文内容为原创，可以自由转载，但请确保注明出处。感谢你的阅读！
浅析samtools mpileup命令的使用及其对应的Python使用接口（一）（小白友好）视频
相关评论：
本文将通过实际案例来直观展示samtools mpileup命令的运行结果，并尝试详细解释一些细节。首先，我选择了10x Genomics的1k PBMC数据集作为测试数据，你可以通过以下shell命令下载这个数据集：
为了使运行结果更加直观，我从下载到的BAM文件中提取了几个比对记录，并将它们整合到一个新的文件中。这个过程包括以下两步：
此外，我还单独下载了人的一号染色体的参考序列，以便于后续分析。现在，让我们来详细了解samtools mpileup命令的基本概念和使用方法。
pileup命令意为“堆叠”，它能够在给定基因组区间内的每个碱基位置（column）上，将映射到该位置上的所有reads堆叠起来，以展示输入BAM文件在指定基因组区间内的每个碱基位置的总体信息。
要开始使用samtools mpileup，首先需要设置基因组区间。可以通过-r/--region参数来指定。例如：
请注意，samtools mpileup默认会忽略一些alignment，包括UNMAP、SECONDARY、QCFAIL和DUP，这些alignment的FLAG值已被标记为256，即表示次要比对。如果你的示例BAM文件中包含这些类型的alignment，那么它们不会被输出。你可以通过重新设置参数--ff UNMAP,QCFAIL,DUP来显示这些被标记为次要比对的alignment。
输出结果中的字段分别代表以下内容：1）参考序列名称；2）基因组位置坐标；3）参考基因组上该位置的碱基；4）覆盖该基因组位置的reads数量；5）每个reads在该位置的碱基；6）碱基的质量分数。
同时，samtools mpileup命令能够正确识别并忽略reads开头被soft clip的一个碱基。此外，你还可以通过-l/--positions参数接BED文件来指定基因组区间，如果同时使用了-r和-l，那么最终生效的区间是两者交集。
在使用samtools mpileup时，还可以设置过滤条件以筛选alignment。例如，可以设置最低的Mapping Quality阈值（-q --min-MQ）和碱基质量分数阈值（-Q --min-BQ），或者排除特定的Read Group（-G --exclude-RG）。
为了更详细地分析结果，samtools mpileup命令在运行时提供了输出参考基因组序列的功能，通过-f / --fasta-ref参数来提供参考序列Fasta文件。这将改变输出结果中某些字段的表示方式。
在探索samtools mpileup的使用方法时，我尽量覆盖了这一命令的大部分用法。对于更深入的细节和使用技巧，我鼓励小伙伴们参考官方文档进行学习。同时，如果你在使用过程中发现任何不妥之处，欢迎提出宝贵的意见和建议。
本文内容为原创，可以自由转载，但请确保注明出处。感谢你的阅读！
浅析samtools mpileup命令的使用及其对应的Python使用接口（一）（小白友好）视频
相关评论：
