CRISPR Screen | 基因编辑筛选靶点 | MAGeCK | MAGeCKFlute

 

Twist Tier 4: 1001-2000 Oligos
除去10%的NTC，实际上的容量是300个基因（包括两个reporter）

打折前：1,645.00 USD

打折后：1,325.00 USD

 

2024年09月17日

正式开始small focus library screen

Library Mode

240 genes
6 gRNA per gene
160 NTC (240*6*10/9 * 10%)
Half # No-Site Controls + half # Intergenic Controls

 

add reporter gRNA after CRIPSR pick (here is mKate2)

check all 6 mKate2 gRNA in two plasmids

 

get results from CRISPRpick and copy gene name and gRNA to the new Excel file

add gRNA for reporter

go to twist to make the order

 

参考：/Users/zhixinli/Partners HealthCare Dropbox/Zhixin Li/LI-LAB-v0/MES-lab/CRISPR screen 

 

2024年09月06日

准备做一个focused library CRISPR screen做一个强有力的验证。

目前主要有三种sgRNA的来源：

  • https://www.addgene.org/pooled-library/liu-crispr-knockout-h3/
  • Brunello CP0043，https://www.addgene.org/pooled-library/broadgpp-human-knockout-brunello/
  • https://portals.broadinstitute.org/gppx/crispick/public （Human GRCh38(Ensembl v.112)、 CRISPRko、 SpyoCas9(NGG)、 Chen (2013) tracrRNA、 6 sgRNA per gene）

order the oligos from Twist Company

 

2024年08月27日

为了检查单细胞推理出来的基因表达相关性，与CRISPR的结果是否一致，我做了一个GSEA富集分析，结果correlation的top1-3显著富集，把我乐坏了，结果仔细一看符号方向不对。

这个问题困扰了我快一下午了，single和dual的结果也基本没有什么相关性，一些常见的hit也不符合常识。

后来在回家的路上我才慢慢意识到，sample size的重要性，这个已经是genome-wide的screen了，与GWAS异曲同工，GWAS是population的数据，缺点是sample的异质性太强，所以需要很大的sample size才能找到显著的hit。

回到CRISPR screen，有时我们只有1 vs 1的样本数，最多的时候就是3 vs 3，虽然样本本身没有异质性，但CRISPR这套系统的异质性也很强，Cas9是否表达，侵染了到底几个gRNA？细胞间是否存在相互作用？测序深度是否足够？

所有这些不确定因素导致CRISPR screen的结果非常dirty，SOX9和KRT20本身因为是直接的Target，所以富集也毫无意外，其他的间接的则没有那么明显。

 

所以，我的结论：目前该lab对CRISPR screen的理解和掌控其实是在非常肤浅的程度，数据基本无法利用。

如果是我，我会首先无脑加大样本量，去掉所有复杂的设计（这套系统准确度根本就不高），直接来个20-50组mKate High和Low的两组，这样我才敢说有足够的power来产生足够可信的hit。

实验方向的优化也是可以的，只是更加的费时费力，需要精通实验的experts。

 

这个lab其实非常菜，大部分时候仅仅停留在能懂的层次，偶尔进入会做的层次，但远远没达到精通的地步，这就导致太理论化、理想化，复杂fancy的设计，出来的全是一堆狗屁不同的结果。

举几个例子：

  1. 单细胞为了省钱，做index，人家GB上做得很好，我们做出来就是无法拆分，浪费钱；
  2. cDNA over expression，SOX9的variant就上了5个，结果RNA、cutrun的结果全是shit，OE系统不稳定；
  3. CRISPR screen做成功了一个epigenetic的小library，就直接上genome-wide，还搞出非常复杂fancy的设计，结果证明都是shit，系统灵敏度不行；
  4. cutrun直接上BAF各个subunit的抗体，抗体没有优化，这套体系的灵敏度远远不够，数据太dirty，根本无法回答某些具体的问题；

这种新lab的质控的优化基本没有，既没有什么前沿的课题和idea，也没有什么非常擅长的技术，老板把他过去“成功”的那一套体系全丢了，继续瞎玩吧，你还有多少时间可以浪费？

 

2024年08月22日

这次算是彻底把CRISPR Screen的原理和数据搞清楚了，但也开始意识到数据的局限性，以及实验设计的重要性。

核心就是让Cas9在cell line/organoid里表达，然后把sgRNA库导入进细胞，一般是处理3-7天，然后收集特定的细胞（可以直接是viability or flow sorting），靶向扩增sgRNA，测序分析。

目前主要有两种readout：

  1. viability，这个最简单，7天后，直接靶向扩增测序，相比day0，day7会死很多cell，通过对比sgRNA的丰度来确定dependency，这部分cell line已经被DepMap做完了；
  2. reporter，细胞内荧光标记，以及细胞膜表面的蛋白，就是两个细胞群High和Low，通过对比sgRNA的丰度来确定该基因的调控因子，这也就是实验设计个性化的地方；

 

建库时，一般为了省钱，都会pool在一起，然后拆分fastq。

然后用MAGeCK count来生成每个sgRNA的read count的table。

最后就是MAGeCK的 rra 或者 mle 来分析数据，处理的数据就是每个基因多个sgRNA，会综合考虑然后生成一个log2FC，然后生成Rank。

 

深度理解CRISPR Screen数据

在看原始数据的时候有很多tricky的地方，这个非常类似RNA-seq数据，连normalization都可以用CPM，不同的是一个gene有多个sgRNA。

稍微有点复杂的就是：

  • 每个gene的sgRNA的数量不一样，有的多达几十个，有的就4个；
  • sgRNA的丰度不一样，1 vs 10,100 vs 150，到底那个差异更大？
  • sgRNA的趋势不一样，同一个基因的不同sgRNA的趋势是相反的，怎么搞？

 

彻底搞懂之后，我发现CRISPR Screen很容易非常差，一些条件必须满足，否咋数据可信度很差：

  • CRISPR这个系统是否真的在cell里work，Cas9是否稳定表达
  • sgRNA的library如何？是否合理有效？
  • 细胞数据是否足够，处理时间是否合理？

 

有的分析取rank（P-value），有的取beta，有的就用log2FC。

我个人是觉得可信度不大，一个常识就是只有样本量足够大，RNA-seq的DEG P-value才有power。

这里的结果就只有一个样本，单样本DEG，你信吗？然后sgRNA数量多，数据参差不齐，出来的结果有个鬼的可信度。

 

有个做了3个replicate，还算比较合理。

 

关于SOX9和KRT20 dual reporter

SOX9 -》KRT20这个regulation本身是存在的，所以筛KRT20本身也会筛SOX9，反之则不成立。

 

数据分析参考：projects/SOX9_KRT20_CRISPR_Screen/whole_genome_CRISPR_screen.ipynb

分析工具：

  • https://sourceforge.net/p/mageck/wiki/Home/
  • https://sourceforge.net/p/mageck/wiki/advanced_tutorial/#tutorial-4-make-full-use-of-mageck-mle-for-more-complicated-experimental-design-eg-paired-samples-time-series

Paper参考：

可视化：

 

 

 

工具篇：MAGeCKFlute - Integrative analysis pipeline for pooled CRISPR functional genetic screens

安装折腾了几个小时，debug核心，根据关键提示Google，筛选找到答案即可。

error:

Error: package or namespace load failed for ‘depmap’:
.onLoad failed in loadNamespace() for 'depmap', details:
call: (function (cond)
error: error in evaluating the argument 'x' in selecting a method for function 'query': Failed to collect lazy table.
Caused by error in `db_collect()`:
! Arguments in `...` must be used.
✖ Problematic argument:
• ..1 = Inf
ℹ Did you misspell an argument name?
Error: loading failed
Execution halted
ERROR: loading failed


Solved:
downgrading dbplyr solves this issue

devtools::install_version("dbplyr", version = "2.3.4")

　　

 

 

 

系统性的学习一下CRISPR screen的数据分析

前篇：fastq demultiplexing - 确实得到了比较靠谱的fastq文件

 

这次比较麻烦，需要从Undetermined的fastq开始，问题就是index是RC，当时没有考虑到。

自己做有点麻烦，如果直接用最短的index，会出现很多假阳性，所以必须extend。

index在read的header里，因为是dual index，所以需要准备index fastq，然后用工具提取即可，一晚上能搞完80G+80G的paired fastqs。 

 

fastq非常dirty，需要先trim一下，fastp一下，80G只剩下65G。【完全可以不做】下面的方法花里胡哨，其实没啥用。

Find a tar file in the above location. It has 12 samples from HEK293T and an excel with sample information.

I trimmed the reads using fastp to remove adapters, poly-G tails and low-quality bases. I used PEAR to merge mate pairs to make a single read based on default overlap parameters.

The fastq.gz files you have now, treat them as single end read files. I believe you have the inzolia guide file for annotation, I shared the simple csv for MAGECX in the same folder. Will keep HT29 data in the same folder once I am done preprocessing it.

 【可以略过，没必要做merge，知道有这么个工具就行】

其次就是merge，因为测序长度可能超过了insert size，所以需要寻找overlap，然后merge，这个PEAR可以做。【bbmerge.sh也能做，但比较旧】

R1 --------------------------->
                                       <------------------------- R2

但PEAR会剩下很多unassembled的reads，这时候就用BBtools合并一下，reformat.sh，就是粗暴的拼接。【reformat就是简单的cat，并不是RC后paste】

reformat.sh in1=sample_R1.fq.gz in2=sample_R2.fq.gz out=sample.fq.gz

　　

所以强制改了PEAR的参数，使用所有的reads，这样就不会有unassembled的reads了。

pear -n 0 -e -p 1.0 -v 0 -j 20 -f merged.fq.unassembled.forward.fastq -r merged.fq.unassembled.reverse.fastq -o unassembled

　　

 

 

 

搞完之后，就可以跑MAGeCK的pipeline了。

实操分享-使用MAGeCK分析Bulk CRISPR Screen数据

conda create -c bioconda -n crispr_screen mageck
conda activate crispr_screen
conda update mageck

　　

可以直接用paired fastq

mageck count -l library.txt -n prefix --sample-label A,B,C,D --fastq A_1.fq.gz B_1.fq.gz C_1.fq.gz D_1.fq.gz --fastq-2 A_2.fq.gz B_2.fq.gz C_2.fq.gz D_2.fq.gz

 

首先要准备cas12_library_all.csv文件

这次的cas12 lib比较特殊，是四个20 bp的sgRNA被拼接成了一个sgRNA，我试了用bowtie2和hisat2去比对到全长sgRNA，基本是0比对。

conda install bioinfo::hisat
conda install bioconda::hisat2

hisat2-build ../library_ref.fasta gRNA_ref

hisat2 -p 15 -t --very-sensitive \
    -x /home/zz950/projects/demultiplexing/demultiplexed/sgRNA_ref/gRNA_ref \
    -1 fastq/D7_1uM_1_1.fq \
    -2 fastq/D7_1uM_1_2.fq \
    -S D7_1uM_1.sam 2> D7_1uM_1.GenomeMapStat.xls


bowtie2-build ../bowtie2_ref.fasta gRNA_ref

bowtie2 -x /home/zz950/projects/demultiplexing/demultiplexed/sgRNA_ref_bowtie2/gRNA_ref -1 fastq/D7_1uM_1_1.fq -2 fastq/D7_1uM_1_2.fq -p 15 -S D7_1uM_1.sam

　　

所以不得不拆分成四个sgRNA，然后继续用mageck来分析

mageck说是可以用fastq2，其实基本所有的信息都在fastq1里面，用了paired fastq反而会报错。

我们的样本

D7_1uM_1
D7_1uM_2
D7_DMSO_1
D7_DMSO_2
D2_DMSO_1
D2_DMSO_2
D2_0.1uM_1
D2_0.1uM_2
Plasmid

 

一行命令搞定：

mageck count -l ../cas12_library_all.csv -n results/all.cas12 --sample-label D7_1uM_1,D7_1uM_2,D7_DMSO_1,D7_DMSO_2,D2_DMSO_1,D2_DMSO_2,D2_0.1uM_1,D2_0.1uM_2,Plasmid --fastq D7_1uM_1_1.fq D7_1uM_2_1.fq D7_DMSO_1_1.fq D7_DMSO_2_1.fq D2_DMSO_1_1.fq D2_DMSO_2_1.fq D2_0.1uM_1_1.fq D2_0.1uM_2_1.fq Plasmid_1.fq

# 用了会报错
--fastq-2 D7_1uM_1_2.fq D7_1uM_2_2.fq D7_DMSO_1_2.fq D7_DMSO_2_2.fq D2_DMSO_1_2.fq D2_DMSO_2_2.fq D2_0.1uM_1_2.fq D2_0.1uM_2_2.fq Plasmid_2.fq

　　

对每个样本单独做也行

mageck count -l ../cas12_library_all.csv -n results/D7_1uM_1 --sample-label D7_1uM_1 --fastq D7_1uM_1_1.fq
mageck count -l ../cas12_library_all.csv -n results/D7_1uM_2 --sample-label D7_1uM_2 --fastq D7_1uM_2_1.fq --fastq-2 D7_1uM_2_2.fq 
mageck count -l ../cas12_library_all.csv -n results/D7_DMSO_1 --sample-label D7_DMSO_1 --fastq D7_DMSO_1_1.fq --fastq-2 D7_DMSO_1_2.fq 
mageck count -l ../cas12_library_all.csv -n results/D7_DMSO_2 --sample-label D7_DMSO_2 --fastq D7_DMSO_2_1.fq 
mageck count -l ../cas12_library_all.csv -n results/D2_DMSO_1 --sample-label D2_DMSO_1 --fastq D2_DMSO_1_1.fq --fastq-2 D2_DMSO_1_2.fq 
mageck count -l ../cas12_library_all.csv -n results/D2_DMSO_2 --sample-label D2_DMSO_2 --fastq D2_DMSO_2_1.fq
mageck count -l ../cas12_library_all.csv -n results/D2_0.1uM_1 --sample-label D2_0.1uM_1 --fastq D2_0.1uM_1_1.fq
mageck count -l ../cas12_library_all.csv -n results/D2_0.1uM_2 --sample-label D2_0.1uM_2 --fastq D2_0.1uM_2_1.fq 
mageck count -l ../cas12_library_all.csv -n results/Plasmid --sample-label Plasmid --fastq Plasmid_1.fq --fastq-2 Plasmid_2.fq

　　

OK，得到count matrix了，可以做下游分析了。

 

讨论之后，发现必须要用top2 gRNA作为construct，来计数。count也没有减少多少。

mageck count -l ../cas12_library_top2.csv -n results/all.top2.cas12 --sample-label D7_1uM_1,D7_1uM_2,D7_DMSO_1,D7_DMSO_2,D2_DMSO_1,D2_DMSO_2,D2_0.1uM_1,D2_0.1uM_2,Plasmid --fastq D7_1uM_1_1.fq D7_1uM_2_1.fq D7_DMSO_1_1.fq D7_DMSO_2_1.fq D2_DMSO_1_1.fq D2_DMSO_2_1.fq D2_0.1uM_1_1.fq D2_0.1uM_2_1.fq Plasmid_1.fq

　

mageck test -k all.top2.cas12.count.merge.txt -t D7_1uM -c D7_DMSO -n D7_comp
mageck test -k all.top2.cas12.count.merge.txt -t D2_0.1uM -c D2_DMSO -n D2_comp

　　

还有一种对比看时间效应：sample （D2 or D7） vs plasmid（Day 0）

但因为我们的plasmid数据不行，暂时就不对比了。

 

下游分析代码：projects/demultiplexing/demultiplexed/cas12_SOX9_dTAG_CRISPR_screen.ipynb 

 

从novogene得到了新的 demultiplexing 数据，对比了一下，数据量差不多，但是还是以他们为金标准，重新分析一遍吧，也挺快的。

 

 2057  zcat Inzolia_Plasmid/Inzolia_Plasmid_1.fq.gz Inzolia_Plasmid_Re/Inzolia_Plasmid_Re_1.fq.gz > Inzolia_Plasmid_1.fq
 2059  zcat Inzolia_Plasmid/Inzolia_Plasmid_2.fq.gz Inzolia_Plasmid_Re/Inzolia_Plasmid_Re_2.fq.gz > Inzolia_Plasmid_2.fq
 2061  zcat D7_DMSO_R2/D7_DMSO_R2_1.fq.gz D7_DMSO_R2_Re/D7_DMSO_R2_Re_1.fq.gz > D7_DMSO_R2_1.fq
 2062  zcat D7_DMSO_R2/D7_DMSO_R2_2.fq.gz D7_DMSO_R2_Re/D7_DMSO_R2_Re_2.fq.gz > D7_DMSO_R2_2.fq
 2064  zcat D2_DMSO_R1/D2_DMSO_R1_1.fq.gz D2_DMSO_R1_Re/D2_DMSO_R1_Re_1.fq.gz > D2_DMSO_R1_1.fq
 2065  zcat D2_DMSO_R1/D2_DMSO_R1_2.fq.gz D2_DMSO_R1_Re/D2_DMSO_R1_Re_2.fq.gz > D2_DMSO_R1_2.fq
 2066  zcat D7_1uM_R1/D7_1uM_R1_1.fq.gz D7_1uM_R1_Re/D7_1uM_R1_Re_1.fq.gz > D7_1uM_R1_1.fq
 2067  zcat D7_1uM_R1/D7_1uM_R1_2.fq.gz D7_1uM_R1_Re/D7_1uM_R1_Re_2.fq.gz > D7_1uM_R1_2.fq
 2068  zcat D2_100nM_R2/D2_100nM_R2_1.fq.gz D2_100nM_R2_Re/D2_100nM_R2_Re_1.fq.gz > D2_100nM_R2_1.fq
 2069  zcat D2_100nM_R2/D2_100nM_R2_2.fq.gz D2_100nM_R2_Re/D2_100nM_R2_Re_2.fq.gz > D2_100nM_R2_2.fq
 2070  zcat D7_DMSO_R1/D7_DMSO_R1_1.fq.gz D7_DMSO_R1_Re/D7_DMSO_R1_Re_1.fq.gz > D7_DMSO_R1_1.fq
 2071  zcat D7_DMSO_R1/D7_DMSO_R1_2.fq.gz D7_DMSO_R1_Re/D7_DMSO_R1_Re_2.fq.gz > D7_DMSO_R1_2.fq
 2072  zcat D2_DMSO_R2/D2_DMSO_R2_1.fq.gz D2_DMSO_R2_Re/D2_DMSO_R2_Re_1.fq.gz > D2_DMSO_R2_1.fq
 2073  zcat D2_DMSO_R2/D2_DMSO_R2_2.fq.gz D2_DMSO_R2_Re/D2_DMSO_R2_Re_2.fq.gz > D2_DMSO_R2_2.fq
 2074  zcat D7_1uM_R2/D7_1uM_R2_1.fq.gz D7_1uM_R2_Re/D7_1uM_R2_Re_1.fq.gz > D7_1uM_R2_1.fq
 2075  zcat D7_1uM_R2/D7_1uM_R2_2.fq.gz D7_1uM_R2_Re/D7_1uM_R2_Re_2.fq.gz > D7_1uM_R2_2.fq
 2076  zcat D2_100nM_R1/D2_100nM_R1_1.fq.gz D2_100nM_R1_Re/D2_100nM_R1_Re_1.fq.gz > D2_100nM_R1_1.fq
 2077  zcat D2_100nM_R1/D2_100nM_R1_2.fq.gz D2_100nM_R1_Re/D2_100nM_R1_Re_2.fq.gz > D2_100nM_R1_2.fq

　　

mageck count -l ../../cas12_library_top2.csv -n results/all.2gRNA.cas12 --sample-label D7_1uM_1,D7_1uM_2,D7_DMSO_1,D7_DMSO_2,D2_DMSO_1,D2_DMSO_2,D2_0.1uM_1,D2_0.1uM_2 --fastq D7_1uM_R1_1.fq D7_1uM_R2_1.fq D7_DMSO_R1_1.fq D7_DMSO_R2_1.fq D2_DMSO_R1_1.fq D2_DMSO_R2_1.fq D2_100nM_R1_1.fq D2_100nM_R2_1.fq 
# --fastq-2 D7_1uM_R1_2.fq D7_1uM_R2_2.fq D7_DMSO_R1_2.fq D7_DMSO_R2_2.fq D2_DMSO_R1_2.fq D2_DMSO_R2_2.fq D2_100nM_R1_2.fq D2_100nM_R2_2.fq

mageck test -k all.2gRNA.cas12.count.merge.txt -t D7_1uM -c D7_DMSO -n D7_comp
mageck test -k all.2gRNA.cas12.count.merge.txt -t D2_0.1uM -c D2_DMSO -n D2_comp

　　

mageck count -l ../../cas12_library_top2.csv -n results/all.2gRNA.cas12 --sample-label D7_1uM_1,D7_1uM_2,D7_DMSO_1,D7_DMSO_2,D2_DMSO_1,D2_DMSO_2,D2_0.1uM_1,D2_0.1uM_2,Plasmid --fastq D7_1uM_R1_1.fq D7_1uM_R2_1.fq D7_DMSO_R1_1.fq D7_DMSO_R2_1.fq D2_DMSO_R1_1.fq D2_DMSO_R2_1.fq D2_100nM_R1_1.fq D2_100nM_R2_1.fq Inzolia_Plasmid_1.fq
# --fastq-2 D7_1uM_R1_2.fq D7_1uM_R2_2.fq D7_DMSO_R1_2.fq D7_DMSO_R2_2.fq D2_DMSO_R1_2.fq D2_DMSO_R2_2.fq D2_100nM_R1_2.fq D2_100nM_R2_2.fq

　　

 

CRISRP共同量测序分为阳性筛选和阴性筛选。

阳性筛选指施加一定的筛选压力，经文库扰动后野生型细胞致死，仅有获得抗性的细胞存活。

与阳性筛选不同，阴性筛选是通过比较不同筛选时间点的sgRNA丰度差，获得缺失的sgRNA，分析其所对应的基因，从而找出可能影响细胞生长的候选基因。

MAGeCK is designed to identify positively and negatively selected sgRNAs and genes in genome-scale CRISPR/Cas9 knockout experiments.

Genes are ranked by the p.neg field (by default).

 

结果解读

Output file specification

https://sourceforge.net/p/mageck/wiki/output/#gene_summary_txt 

 

两种不同的compare方法

mageck test -k all.2gRNA.cas12.count.merge.txt -t D7_1uM -c D7_DMSO -n D7_comp
mageck test -k all.2gRNA.cas12.count.merge.txt -t D2_0.1uM -c D2_DMSO -n D2_comp


mageck mle -k all.2gRNA.cas12.count.merge.txt -d D2_designmatrix.txt -n D2_mle_comp
mageck mle -k all.2gRNA.cas12.count.merge.txt -d D7_designmatrix.txt -n D7_mle_comp

　　

第一种rra就是肯定有两种结果，一是positive，一是negative；

第二种就是mle，只有一个beta，出来的火山图比较奇怪，可以对比一下；

 

参考：

  • Crispr-cas9 screen analysis mageck
  • 可以用其他sgrna counter，但我还是喜欢用成熟的mageck
  • cas12 paper: Efficient gene knockout and genetic interactions: the IN4MER CRISPR/Cas12a multiplex knockout platform

 

posted @ 2024-01-31 08:09  Life·Intelligence  阅读(448)  评论(0)  编辑  收藏  举报
TOP