为什么 R 中的 pileup() 函数会过度计数读取？

Question

Pileup() 函数过度计算了变异位置的读取数。（我从某人那里继承了这段代码。）

我查看了 pileup 参数，但 none 似乎导致了多算。我在 IGV 中可视化变体，发现读取数量与 pileup 报告的不匹配。当 Pileup 报告读取次数较少（即 3 或 8）的变体的读取计数时，它是准确的。 Pileup 在报告具有大量读取（即 20、60）的变体的读取计数时不准确。

res <- pileup( bf, snp, 
        scanBamParam=ScanBamParam( flag = scanBamFlag( hasUnmappedMate=F, isProperPair=T, isDuplicate=F ), which = snp ), 
        pileupParam=PileupParam( distinguish_strands=F, min_base_quality=10, max_depth=1e4 ) )

其中 bf 是 bam 文件，snp 是列出找到的所有 snps 的 grange 对象。

> bf
class: BamFile 
path: Sample_4.split.bam
index: Sample_4.split.bai
isOpen: FALSE 
yieldSize: NA 
obeyQname: FALSE 
asMates: FALSE 
qnamePrefixEnd: NA 
qnameSuffixStart: NA 

> head(snp)
GRanges object with 6 ranges and 5 metadata columns:
              seqnames    ranges strand | paramRangeID            REF
                 <Rle> <IRanges>  <Rle> |     <factor> <DNAStringSet>
   rs13475703     chr1   5124338      + |         <NA>              A
    rs6408157     chr1   9547068      + |         <NA>              A

                             ALT      QUAL      FILTER
              <DNAStringSetList> <numeric> <character>
   rs13475703                  G    233.77        PASS
    rs6408157                  G     52.77        PASS

  -------
  seqinfo: 66 sequences from mm10 genome

> res[66,]
   seqnames      pos nucleotide count            which_label
66     chr1 87942764          G    61 chr1:87942764-87942764

对于上面的示例，对于 chr1:87942764 处的 G 变体，pileup 报告了 61 个读取。当我在 IGV 中打开 bam 文件时，我只计算了 41 个支持 G 变体的读取。

没有错误信息。只是报告的计数与我在 IGV 中看到的不匹配。

Answer 1

正如用户 merv 评论的那样：您是否将 IGV 配置为显示所有读数（查看 > 首选项... > 比对轨道选项 > 下采样读数）？ IGV 的默认设置是对读取进行下采样，因此可视化可能没有显示所有实际读取。

当我取消选中 "filter vendor failed reads" 时，读取计数与 pileup 报告的相匹配。谢谢 merv.

为什么 R 中的 pileup() 函数会过度计数读取？

Why is pileup() function in R over-counting reads?

r

function

bioinformatics

counting