使用MAGeCK count分析CRISPR筛选数据时reads与sgRNA丢失问题求助

背景信息 #

• 处理5年前的CRISPR敲除筛选BAM格式测序数据，持有原始测序统计（含总reads数、比对reads数、未比对reads数）
• 使用以下命令完成BAM转FASTQ：

  samtools view -h -F 2048 filename.bam > tmp.bam
  bedtools bamtofastq -i tmp.bam -fq filename.fastq
  

• 核心问题：MAGeCK count输出的总reads数远低于原始统计（如仅5万条vs原始4400万条），77441条sgRNA中约61000条计数为0，重新下载转换数据后问题依旧

可能原因及排查方案 #

1. BAM转换时的过滤参数错误 #

你使用的samtools view -F 2048会过滤掉**补充比对（secondary alignments）**的reads，但原始测序统计的总reads包含所有测序产出序列，若BAM本身已做过过滤，叠加该参数会进一步丢失大量reads。

• 排查：
  • 用samtools view filename.bam | wc -l统计BAM内的比对reads数，和原始统计的比对reads数对比
  • 用samtools view -F 2048 filename.bam | wc -l统计过滤后的reads数，查看被剔除的比例
• 解决：若无需过滤补充比对reads，去掉-F 2048参数，直接执行samtools view -h filename.bam > tmp.bam

2. BAM仅存储了比对成功的reads #

原始测序统计的总reads包含未比对的reads，但如果当时的测序流程仅将比对到参考基因组的reads存入BAM，而CRISPR筛选的sgRNA可能因靶向区域特殊性未比对成功，这部分未比对reads会完全丢失。

• 排查：用samtools flagstat filename.bam查看BAM的总reads数，和原始统计的总reads数对比，若差距极大则说明BAM未保留未比对reads
• 解决：优先寻找原始FASTQ文件直接分析；若原始FASTQ丢失，联系实验室确认当年的数据处理流程是否剔除了未比对reads

3. MAGeCK使用的sgRNA参考库不匹配 #

实验所用的sgRNA库有77441条，但MAGeCK count指定的参考库序列与实际实验库不一致时，会导致大量reads无法匹配，最终计数为0。

• 排查：
  • 随机抽取转换后FASTQ的reads序列，手动与你的sgRNA库序列比对，确认匹配度
  • 检查MAGeCK count命令是否正确指定了实验对应的sgRNA库文件（--library参数）
• 解决：找回5年前实验使用的sgRNA库原始序列，生成匹配的参考库文件后重新运行count

4. BAM中的reads是截断/含接头的比对序列 #

CRISPR筛选测序通常仅靶向sgRNA核心区域（如20bp），但如果BAM是比对到参考基因组后的结果，reads可能被截断或残留接头序列，导致MAGeCK无法识别sgRNA。

• 排查：用samtools view filename.bam | head -10查看reads序列，与原始sgRNA序列对比完整性
• 解决：若存在截断/接头，使用MAGeCK的--trim-5/--trim-3参数修剪reads后再计数；或优先使用原始FASTQ重新处理

5. 旧测序数据的质量降解 #

5年前的旧数据可能存在质量降解，或当年测序流程存在缺陷，导致大量reads无法匹配sgRNA。

• 排查：用fastqc filename.fastq分析转换后FASTQ的质量分布，统计能匹配到sgRNA库的reads比例
• 解决：用fastp等工具对FASTQ做质量过滤（去除低质量reads、接头）后再运行MAGeCK count；若质量极差且样本仍保留，可考虑重新测序

内容的提问来源于stack exchange，提问作者Cassie Bishop