使用GenomicAlignments对BAM文件计数后无输出的原因排查

基于hg38.knownGene.gtf对双端RNA-seq样本完成修剪和比对后得到约2GB的.bam文件，使用GenomicAlignments包进行reads计数时输出文件为空。比对步骤使用了相同的参考文件，RStudio为最新版，计数代码如下：

input<- ("input.bam")
#total size of BAM file is around 2GB so, I used yieldSize to call .bam file.
data.FL <- BamFileList(input, yieldSize = 3000000) 
# reference file
gtf.file<-"hg38.knownGene.gtf"  #ucsc
TxDb<-makeTxDbFromGFF(gtf.file,format="gtf")

#GrangesList
exons<-exonsBy(TxDb,by="gene")
transcripts<-transcriptsBy(TxDb,by="gene") 
grl<-GRangesList(c(exons,transcripts))

#counting
data<-summarizeOverlaps(grl,data.FL,mode = "Union",singleEnd=F,ignore.strand=TRUE)

#writing the output in .csv:
vector <- as.vector(mcols(data)$data)
write.csv(vector,"vector.csv")

• GRangesList构建逻辑错误：直接将exons和transcripts用c()合并后转GRangesList，会导致区域结构混乱，且同时用外显子和转录本区域计数会造成区域重复，干扰reads匹配。应只选择一种区域（优先外显子）构建计数区域：

  # 仅保留外显子区域用于计数
  grl <- exonsBy(TxDb, by = "gene")
  

• 染色体命名不匹配：UCSC的hg38.gtf染色体名为chr1、chr2格式，但若比对工具输出的BAM文件使用无前缀的1、2命名，会导致注释区域与reads比对位置无法匹配。可分别检查两者的染色体命名：

  # 查看BAM文件的染色体信息
  seqinfo(BamFile(input))
  # 查看TxDb的染色体信息
  seqlevels(TxDb)
  

  若不一致，需统一命名（如给BAM文件染色体名添加chr前缀，或修改TxDb的seqlevels）。
• 计数结果提取方式错误：summarizeOverlaps返回的是SummarizedExperiment对象，正确提取计数矩阵的方式是assay(data)，而非mcols(data)$data——后者本身会返回空值：

  count_matrix <- assay(data)
  write.csv(count_matrix, "counts.csv")
  

• BAM索引文件缺失：确保BAM文件对应的.bai索引文件存在，且与BAM文件同名、放在同一目录下，GenomicAlignments读取BAM必须依赖索引文件快速定位区域。
• yieldSize分批读取异常：可先尝试不设置yieldSize，直接读取整个BAM文件测试是否能得到计数，排除分批读取导致的问题：

  data.FL <- BamFileList(input)
  

• reads与注释区域无重叠：读取少量reads检查比对位置是否落在注释区域内：

  gal <- readGAlignments(BamFile(input), n = 1000)
  hits <- findOverlaps(gal, grl)
  # 查看匹配到区域的reads数量
  length(hits)
  

  若结果为0，说明比对用的参考基因组与gtf文件不匹配（如混用NCBI和UCSC版本的hg38），需重新确认比对步骤的参考序列是否与gtf对应。

内容的提问来源于stack exchange，提问作者mehdi heidari