假设我有一个 BAM 文件和几个位置，我想在此对齐中更仔细地检查它们。我的目标是找出这些位置是否在相同的读数上，以及每个读数在这两个位置上有哪些核苷酸。所以我的输出应该是在位置上出现的所有读数的表格，并列出相关的核苷酸。
输出应如下所示（前三个位置为 980、1350 和 1400；我组成了读取名称）：
| 读名字 | 980 | 1350 | 1400 | ... |
|---|---|---|---|---|
| ReAdN4m1L03 | 一个 | 吨 | ... | |
| ReAdN4m2L02 | 一个 | 吨 | G | ... |
| ReAdN4m3L01 | 吨 | G | ... | |
| ... | ... | ... | ... | ... |
最后，我想要一张通过读数连接核苷酸的表格。
到目前为止，我已尝试在所需位置读取 Bam 文件：
#Positions of interest
pos <- c(980, 1350, 1400)
#path to BAM file
path <- pathToBamFile
#name of reference used for mapping
refname <- "myreference"
#loading the bam file
bamFile <- BamFile(path)
gr <- GRanges(seqnames = refname, ranges = IRanges(start = pos, end = pos))
params <- ScanBamParam(which = gr, what = scanBamWhat())
aln <- scanBam(bamFile, param = params)
使用以下代码，我可以访问读取，但是如何在各个位置获取相应的核苷酸？
pos.no <- c(1:length(pos))
i <- 1
aln[[pos.no[i]]]$qname
我已经尝试使用aln[[i]]$cigar带有序列的 CIGAR 字符串aln[[i]]$seq来手动提取核苷酸，但这导致了越来越多的问题。
有没有什么简单的方法可以提取核苷酸？
我现在在互联网上搜索了很长时间，以寻找一个简单的解决方案，但没有成功。希望有人会在这里帮助我！非常感谢您提前抽出时间。