Bash：根据读取长度从 BAM 文件中提取读取答案

【问题标题】：Bash: Extract reads from BAM files based on read lengthBash：根据读取长度从 BAM 文件中提取读取
【发布时间】：2022-06-10 18:46:41
【问题描述】：

关于堆栈溢出的第一个问题，希望您能帮助我。

假设一个 BAM 文件，我只想从中提取一定长度的读数（42 - 65 nt；第 10 列），但包含其余列的信息。示例 sn-p：

VH00693:3:AAANGKTM5:1:1507:7438:26974_AGTTATAGAC 256 ENST00000438504.2 352 0 32M * 0 0 CCTGCAGGAATATGGCTCCATCTTCATGGGCG CCCCCCCCCCC;CCCCCCCCCCCCCCCCCCCC NH:i:50 HI:i:4

VH00693:3:AAANGKTM5:1:1507:7438:26974_AGTTATAGAC 256 ENST00000438504.1 352 0 32M * 0 0 CCTGCAGGAATATGGCTCCATCTTCATGGGCGCCTGCAGGAATATGGCTCCATCTTCATGGGCG CCCCCCCCCCC;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC;CCCCCCCCCCCCCCCCCCCC NH:i:50 HI:i:4

我的尝试是使用 samtools 视图访问 BAM 文件 (Initial.bam)，并搜索适合所需读取大小的子字符串，然后将其解析为新的 BAM 文件 (Extract.bam）。

samtools view -h Initial.bam | \awk 'substr($0,1,1)=="@" || ($10>=42 && $10<=65)'| \samtools view -b > Extract.bam

但是，Extract.bam 仅包含 Initial.bam 的已提取标题部分（以“@”开头）。因此，标题提取工作，以及解析成一个新的 bam 文件。初始文件确实包含所需范围的读取，但那时我不知道如何调整我的代码 sn-p。你有什么建议吗？

【问题讨论】：

标签： bash extract bam

【解决方案1】：

如果您没有gawk，这应该适用于其他人

samtools view -h Initial.bam |

mawk '((_=length($10))%__-(_-__)%24)</^[@]/' __=42 | samtools view -b

【讨论】：

【解决方案2】：

找到了一个能够提供所需输出的改编版本。

samtools view -h Initial.bam | awk '{ if ( substr($0,1,1) == "@" || (length($10) >= 42 && length($10) <= 65)) {print $0} }' | samtools view -b > Extract.bam

【讨论】：