3
$\begingroup$

Sorry if this question seems basic. I am trying to run the MAGeCK pipeline to analyze CRISPR knockout screen data produced by someone in my lab around 5 years ago. I was given the data as bam files and also have the sequencing sample statistics from when the samples were run including the number of total reads, aligned reads, and unaligned reads. I converted the bam files back to fastq files using

samtools view -h -F 2048 filename.bam > tmp.bam

bedtools bamtofastq -i tmp.bam -fq filename.fastq

After converting the bam files back to fastq files, I ran the MAGeCK count function. I saw that none of the total reads matched between the MAGeCK count summary output file and the original sequencing statistics.

More problematically, a few samples had significantly fewer total reads (like 50,000) shown on the MAGeCK count summary output files relative to the original sequencing statistics (~44 million). As a result, many of the sgRNAs are not represented at all in the sample (there are ~61,000 zero counts relative to the total 77441 sgRNAs) which is impacting my analysis.

Can anyone help me understand why I might be losing sgRNAs using MAGeCK? Is it a problem with my source files? I would appreciate any help!

$\endgroup$
1
  • $\begingroup$ I know this is quite an old thread, but as I am experiencing something very similar, can I kindly ask you if you were eventually able to solve the issue and how? Thank you very much in advance! $\endgroup$
    – Marty
    Commented Jun 18 at 20:23

1 Answer 1

1
$\begingroup$

It's hard to say without looking more closely at the data, but you may be able to find a few guides that are drastically different in your counting, and see how prevalent they are within the original BAM file (before filtering). Perhaps the -F 2048 being used to filter out chimeric reads is overfiltering, and causing the discrepancy? Are unaligned reads included in the BAM - if not, does this explain the discrepancy?

Since the total reads doesn't match the original, I'd look into this. Are the extracted FASTQs consistent read length? Maybe some trimming was performed that removed the guides?

Do you know what reference the BAM files were mapped to? Perhaps this reference did not include the guide sequences, and were mapped to the genome instead? You may be able to extract this information from the BAM like this:

 samtools view -H filename.bam > filename.header

Ideally you'd have the original FASTQs to work with, but it sounds like those are unavailable...

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.